ACCEPTED FOR PUBLICATION IN IEEE TRANSACTIONS ON POWER SYSTEMS

1

Convex Relaxations of Chance Constrained
AC Optimal Power Flow

arXiv:1702.08372v3 [cs.SY] 4 Oct 2017

Andreas Venzke, Student Member, IEEE, Lejla Halilbasic, Student Member, IEEE, Uros Markovic, Student
Member, IEEE, Gabriela Hug, Senior Member, IEEE, and Spyros Chatzivasileiadis, Member, IEEE

Abstract—High penetration of renewable energy sources and
the increasing share of stochastic loads require the explicit representation of uncertainty in tools such as the optimal power flow
(OPF). Current approaches follow either a linearized approach or
an iterative approximation of non-linearities. This paper proposes
a semidefinite relaxation of a chance constrained AC-OPF which
is able to provide guarantees for global optimality. Using a
piecewise affine policy, we can ensure tractability, accurately
model large power deviations, and determine suitable corrective
control policies for active power, reactive power, and voltage.
We state a tractable formulation for two types of uncertainty
sets. Using a scenario-based approach and making no prior
assumptions about the probability distribution of the forecast
errors, we obtain a robust formulation for a rectangular uncertainty set. Alternatively, assuming a Gaussian distribution of the
forecast errors, we propose an analytical reformulation of the
chance constraints suitable for semidefinite programming. We
demonstrate the performance of our approach on the IEEE 24
and 118 bus system using realistic day-ahead forecast data and
obtain tight near-global optimality guarantees.
Index Terms—AC optimal power flow, convex optimization,
chance constraints, semidefinite programming, uncertainty.

I. I NTRODUCTION

P

OWER system operators have to deal with higher degrees
of uncertainty in operation and planning. If uncertainty is
not explicitly considered, increasing shares of unpredictable
renewable generation and stochastic loads, such as electric
vehicles, can lead to higher costs and jeopardize system
security. The scope of this work is to introduce a convex
AC optimal power flow (OPF) formulation which is able to
accurately model the effect of forecast errors on the power
flow, can define a-priori suitable corrective control policies
for active power, reactive power, and voltage, and can provide
near-global optimality guarantees.
Chance constraints are included in the OPF formulation
to account for uncertainty in power injections, defining a
maximum allowable probability of constraint violation. It is
generally agreed that the non-linear nature of the AC-OPF
along with the probabilistic constraints render the problem for
most instances intractable [1]. To ensure tractability of these
constraints, either a data-driven or scenario-based approach is
applied, or the assumption of specific uncertainty distributions
is required for an analytical reformulation of the chance
constraints. To deal with the higher complexity of chance
A. Venzke, L. Halilbasic and S. Chatzivasileiadis are with the Department
of Electrical Engineering, Technical University of Denmark, Kongens Lyngby,
Denmark.
U. Markovic and G. Hug are with the Power Systems Laboratory, ETH
Zurich, Zurich, Switzerland.

constrained OPF, existing approaches either assume a DCOPF [2]–[6], a linearized AC-OPF [7]–[10] or solve iteratively
linearized instances of the non-linear AC-OPF [11], [12].
Chance constrained DC-OPF results to a faster and more
scalable algorithm, but it is an approximation that neglects
losses, reactive power, and voltage constraints, and can exhibit
substantial approximation errors [13].
Refs. [2] and [3] formulate a chance constrained DCOPF assuming a Gaussian distribution of the forecast errors.
The work in [2] relies on a cutting-plane algorithm to solve
the resulting optimization problem, whereas the work in [3]
states a direct analytical reformulation of the same chance
constraints. This framework is further extended by the work
in [4] which assumes uncertainty sets for both the mean and
the variance of the underlying Gaussian distributions to obtain
a more distributionally robust formulation. The work in [5]
formulates a robust multi-period chance constrained DC-OPF
assuming interval bounds on uncertain wind infeeds. These
works [2]–[5] include corrective control of the generation units
to restore the active power system balance as a function of the
forecast errors. The work in [6] extends this corrective control
framework to include HVDC converter active power set-points
and phase shifting transformers in an N-1 security context.
Alternatively, the works in [7]–[10] use a linearization of the
AC power flow equations based on [14] to achieve a tractable
formulation of the chance constraints. As the operating point
is not known a-priori, the linearization is performed around
a flat start or no-load voltage, and not the actual operating
point. These works [7]–[10] focus on low-voltage distribution
systems with high share of photovoltaic (PV) production and
minimize PV curtailment subject to chance constraints on
voltage magnitudes. Scenario-based methods are applied to
achieve a tractable formulation. In this framework, line flow
limits and corrective control from conventional generation are
not considered. Furthermore, the utilized linearization in [7]–
[10] is designed for radial distribution grids and assumes no
voltage control capability of generation units.
In Ref. [11], an iterative back-mapping and linearization of
the full AC power flow equations is used to solve the chance
constrained AC-OPF. The recent work in [12] uses an iterative
procedure to calculate the full Jacobian, which is the exact AC
power flow linearization around the operating point. Assuming
a Gaussian distribution of the forecast errors, an analytical
reformulation of the chance constraints on voltage magnitude
and current line flow is proposed. Although this approach can
be shown to scale well, it is not convex and does not guarantee
convergence.

ACCEPTED FOR PUBLICATION IN IEEE TRANSACTIONS ON POWER SYSTEMS

(I) Non-convex
AC-OPF

(II) Non-convex chance
constrained AC-OPF
parametrize solution space

(IV) Convex chance
constrained AC-OPF
using affine policy

remove
rank-1

(III) Non-convex chance
constrained AC-OPF
using affine policy

relaxation gap

Fig. 1. We restrict the solution space of the non-linear chance constrained
AC-OPF to the parametrization by the affine policy. This problem is relaxed
by dropping the non-convex rank constraint. With relaxation gap we refer to
the gap between problems (IV) and (III).

In this work, we formulate convex relaxations of chance
constrained AC-OPF which allow us to provide guarantees
for the optimality of the solution, or otherwise upper-bound
the distance to the global optimum of the original non-linear
problem. Besides that, we include chance constraints for all
relevant state variables, namely active and reactive power,
voltage magnitudes and active and apparent branch flows. Two
tractable formulations of the chance constraints are proposed.
First, based on realistic forecast data and making no prior
assumptions about the probability distributions, we formulate
a rectangular uncertainty set and, subsequently, the associated
chance constraints. Second, assuming a Gaussian distribution
of the forecast errors, we provide an analytical reformulation
of the chance constraints.
A. Convex Relaxations and Relaxation Gap
In general, the AC-OPF is a non-convex, non-linear problem. As a result, identified solutions are not guaranteed to
be globally optimal and the distance to the global optimum
cannot be specified. Recent advancements in the area of
convex optimization with polynomials have achieved to relax
the non-linear, non-convex optimal power flow problem and
transform it to a convex semidefinite (SDP) or second-order
cone problem [15]–[17]. Formulating a convex optimization
problem results in tractable solution algorithms that can determine the global minimum. Within power systems, finding the
global minimum has two important implications. First, from
an economic point of view, it can result to substantial cost
savings [18]. Second, from a technical point of view, the global
optimum determines a lower or an upper bound of the required
control effort. The term relaxation gap denotes the difference
between the minimum obtained through the convex relaxation
and the global minimum of the original non-convex problem.
A relaxation is tight, if the relaxation gap is small. A relaxation
is exact, if the relaxation gap is zero, i.e. zero relaxation
gap is achieved when the minimum of the convex relaxation
coincides with the global minimum of the original non-convex,
non-linear problem. Since the work in [19] has shown cases in
which the semidefinite relaxation of [15] fails, it is necessary
to investigate the relaxation gap of the obtained solution,
and examine the conditions under which we can obtain zero
relaxation gap. In the work [20] a reactive power penalty is
introduced, which allows to upper bound the distance to global

2

optimum. In this work, we develop a penalized semidefinite
formulation for a chance constrained AC-OPF, which allows
us as well to determine an upper bound of the distance to
the global optimum. In Fig. 1 we illustrate the previously
explained concepts in the context of our work. With relaxation
gap, we refer to the gap between the semidefinite relaxation
and a non-linear chance constrained AC-OPF which uses the
affine policy to parametrize the solution space.
B. Main Contributions
In this work we propose a framework for a convex chance
constrained AC-OPF. The work in [21] makes a first step
towards such a formulation which takes into account security
constraints and uncertainty. The change of the system state is
described with an affine policy as an explicit function of the
forecast errors. A combination of the scenario approach and
robust optimization is used to ensure tractability of the chance
constraints [22]. For the convex relaxations we build upon the
SDP AC-OPF formulation proposed in [15]. The contributions
of our work are the following:
• To the best of our knowledge, this is the first paper that
proposes a convex formulation for the chance constrained
OPF that (a) is able to determine if it has found the global
minimum of the original non-convex problem1 , and (b)
if not, it is able to determine the distance to the global
minimum through the relaxation gap.
• In this paper, we introduce a penalty term on power
losses which allows us to obtain near-global optimality
guarantees and we investigate the conditions under which
we can obtain a zero relaxation gap.We show that this
penalty term is small in practice, leading to tight nearglobal optimality guarantees of the obtained solution.
• We formulate tractable chance constraints suitable for
semidefinite programming for two types of uncertainty
sets. First, using a piecewise affine policy, we state a
tractable formulation of the chance constrained AC-OPF
with convex relaxations that makes no prior assumptions
on the type of probability distribution. Using existing
data or scenarios, we determine a rectangular uncertainty
set; as the set and the chance constraints are affine or
convex, we can account for the whole set by enforcing
the chance constraints only at its vertices [22]. Second,
assuming Gaussian distributions, we formulate tractable
chance constraints for the optimal power flow equations
that are suitable for semidefinite programming. In that,
we also assume the correlation of different uncertain
variables. To the best of our knowledge, this is the first
paper that introduces a tractable reformulation of the
chance constrained AC-OPF with convex relaxations for
Gaussian distributions.
• The proposed framework includes corrective control policies related to active and reactive power, and voltage.
• Based on realistic forecast data and the IEEE 118 bus test
case, we compare our approach for both uncertainty sets
to the chance constrained DC-OPF formulation in [5],
1 As we will discuss later in this paper, in cases where a penalized SDP
formulation is necessary, this point corresponds to a near-global minimum.

ACCEPTED FOR PUBLICATION IN IEEE TRANSACTIONS ON POWER SYSTEMS

and the iterative AC-OPF in [12]. Compared to the DCOPF formulation, we find that the formulations proposed
in this paper are more accurate and significantly decrease
constraint violations. For the rectangular uncertainty set,
the affine policy complies with all considered chance
constraints and outperforms all other methods having
the lowest number of constraint violations. At the same
time, we obtain tight near-global optimality guarantees
which ensure that the distance to the global optimum is
smaller than 0.01% of the objective value. For a Gaussian
distribution, both the iterative AC-OPF and our approach
satisfy the constraint violation limit, with our approach
achieving slightly lower costs due to the corrective control capabilities. As the realistic forecast data we used
do not follow a Gaussian distribution, we also observed
that both approaches may exceed the constraint violation
limit at certain timesteps for that dataset.
The remainder of this work is structured as follows: In
Section II the convex relaxation of the chance constrained
AC-OPF problem is formulated. Section III introduces the
piecewise affine policy, defines corrective control policies and
states the tractable OPF formulation for both uncertainty sets.
Section IV states an alternative approach using a linearization
based on power transfer distribution factors (PTDFs). Section V investigates the relaxation gap for a IEEE 24 bus system
and presents numerical results for a IEEE 118 bus system
using realistic forecast data. Section VI concludes the paper.
The nomenclature is provided in Table I. An underline and
overline denote, respectively, the upper and lower bound of a
variable.

II. O PTIMAL P OWER F LOW F ORMULATION
A. Convex Relaxation of AC Optimal Power Flow
For completeness, we outline the semidefinite relaxation
of the AC-OPF problem as formulated in [15]. A power
grid consists of N buses and L lines. The set of generator
buses is denoted with G. The following auxiliary variables are
introduced for each bus k ∈ N and line (l, m) ∈ L:
Yk := ek eTk Y
Ylm :=
Yk :=
Ylm :=
Ȳk :=
Mk :=

(ȳlm + ylm )el eTl − (ylm )el eTm


1 <{Yk + YkT } ={YkT − Yk }
2 ={Yk − YkT } <{Yk + YkT }


T
T
1 <{Ylm + Ylm
} ={Ylm
− Ylm }
T
T
} <{Ylm + Ylm
}
2 ={Ylm − Ylm


−1 ={Yk + YkT } <{Yk − YkT }
2 <{YkT − Yk } ={Yk + YkT }
 T

ek ek
0
0
ek eTk
T

X := [<{V} ={V}]

(1)
(2)
(3)
(4)
(5)

3

TABLE I
N OMENCLATURE

Power grid
N , L, G
ck2 , ck1 , ck0
Y
ȳlm , ylm
nb
xlm
BAC
PTDFlm

Set of buses, lines and generators in the power network
Quadratic, linear and constant cost term of generator k
Admittance matrix
Shunt and series admittance of line (l, m)
Number of buses in the power network
Reactance of line (l, m)
Admittance matrix based on DC approximation
Power transfer distribution factor for line (l, m)
Optimal power flow

PGk , QGk
Vk
Plm , Slm
V
X
PDk , QDk
W
W (ζi )
W0
Bi
Biu , Bil
Wv
ρ
δopt

Active and reactive power generation at bus k
Voltage magnitude at bus k
Active and apparent branch flow on line (l, m)
Complex bus voltage vector
Real and imaginary bus voltage vector
Active and reactive power consumption at bus k
Matrix with product of voltages
Matrix W as a function of the forecast errors
Matrix W for forecasted system state
System change for forecast error i
System change for upper/lower limit on forecast error i
Matrix W for vertex v
Ratio of second to third eigenvalue of W
Near-global optimality measure
Uncertainty modeling

nW , W
PW
f
PW
ζ

dG
dW
γ
µ
cos(φ)
τ
Ns
β
v
nv , V
ζv
Λ
λ, η
κ

Number of wind farms and set of buses with wind farms
Wind infeeds
Forecasted wind infeeds
Wind forecast errors
Maximum violation probability of chance constraints
Generator participation factors
Wind deviation vector
Slack variable on generator participation factor
Weight for power loss penalty
Power factor of wind farms
Ratio of maximum reactive to active power
Number of scenarios
Confidence parameter
Vertices of rectangular uncertainty set
Number and set of vertices v
Forecast error for vertex v of uncertainty set
Covariance matrix
Eigenvalues and eigenvectors of covariance matrix
Limit on Gaussian forecast error

complex bus voltages. The non-linear AC-OPF problem can
be written using (1) – (7) as
X
min
{ck2 (Tr{Yk W } + PDk )2 +
W

k∈G

ck1 (Tr{Yk W } + PDk ) + ck0 }

(8)

subject to the following constraints for each bus k ∈ N and
line (l, m) ∈ L:
P Gk − PDk ≤ Tr{Yk W } ≤ P Gk − PDk

(9)

(6)

QG − QDk ≤ Tr{Ȳk W } ≤ QGk − QDk

(10)

(7)

2
k

(11)

−P lm ≤ Tr{Ylm W } ≤ P lm

(12)

2
S lm
T

(13)

Matrix Y denotes the bus admittance matrix of the power grid,
ek the k-th basis vector, ȳlm the shunt admittance and ylm the
series admittance of line (l, m) ∈ L, and V the vector of

k

V 2k ≤ Tr{Mk W } ≤ V
2

2

Tr{Ylm W } + Tr{Ȳlm W } ≤

W = XX

(14)

ACCEPTED FOR PUBLICATION IN IEEE TRANSACTIONS ON POWER SYSTEMS

The objective (8) minimizes generation cost, where ck2 , ck1
and ck0 are quadratic, linear and constant cost variables
associated with power production of generator k ∈ G.2 The
terms PDk and QDk denote the active and reactive power
consumption at bus k. Constraints (9) and (10) include the
nodal active and reactive power flow balances; P Gk , P Gk ,
QG and QGk are generator limits for minimum and maximum
k
active and reactive power, respectively. The bus voltages are
constrained by (11) with corresponding lower and upper limits
V k , V k . The active and apparent power branch flow Plm and
Slm on line (l, m) ∈ L are limited by P lm (12) and S lm
(13), respectively. To obtain an optimization problem linear
in W , the objective function is reformulated using Schur’s
complement:
X
min
αk
(15)
W, α

k∈G


c Tr{Yk W } + ak
√k1
ck2 Tr{Yk W } + bk

√


ck2 Tr{Yk W } + bk
 0 (16)
−1
√
where ak := −αk + ck0 + ck1 PDk and bk := ck2 PDk . In
addition, the apparent branch flow constraint (13) is rewritten:


Tr{Ylm W } Tr{Ȳlm W }
−(S lm )2
0
Tr{Ylm W }
(17)
−1
0
Tr{Ȳlm W }
0
−1
The non-convex constraint (14) can be expressed by:
W 0

(18)

rank(W ) = 1

(19)

The convex relaxation is introduced by dropping the rank constraint (19), relaxing the non-linear, non-convex AC-OPF to a
convex semidefinite program (SDP). The work in [15] proves
that if the rank of W obtained from the SDP relaxation is 1,
then W is the global optimum of the non-linear, non-convex
AC-OPF and the optimal voltage vector can be computed
following the procedure described in [23].
B. Inclusion of Chance Constraints
Renewable energy sources and stochastic loads introduce
uncertainty in power system operation. To account for uncertainty in bus power injections, we extend the presented OPF
formulation with chance constraints. A number of nW wind
farms are introduced in the power grid at buses k ∈ W and
modeled as
f
PW k = PW
+ ζk
(20)
k
f
where PW are the actual wind infeeds, PW
are the forecasted
values and ζ are the uncertain forecast errors. To simplify
notation, the resulting upper and lower bounds on net active
and reactive power injections are written in compact form as:
f
P k := P Gk − PDk + PW
+ ζk ,
k

(21)

f
P k := P Gk − PDk + PW
+ ζk
k

(22)

Qk := QGk − QDk

(23)

Qk := QG − QDk

(24)

k

2 In case renewable curtailment costs are assumed, this could introduce
negative linear costs, which may not result in a tight relaxation.

4

The convex chance constrained AC-OPF problem includes
chance constraints for each bus k ∈ N and line (l, m) ∈ L:
X
min
αk
(25)
W, α

k∈G

s.t. (9), (10), (11), (12), (17), (16), (18) for W = W0 (26)
n
(27)
P P k ≤ Tr{Yk W (ζ)} ≤ P k ,
Qk ≤ Tr{Ȳk W (ζ)} ≤ Qk ,
2
k

V ≤ Tr{Mk W (ζ)} ≤ V

(28)

2
k,

(29)

− P lm ≤ Tr{Ylm W (ζ)} ≤ P lm ,
"
2

Tr{Ylm W (ζ)} Tr{Ȳlm W (ζ)}
−(S lm )
Tr{Ylm W (ζ)}
−1
0
Tr{Ȳlm W (ζ)}
0
−1

W (ζ)  0

o

(30)
#
 0,

≥1−

(31)
(32)

The parameter  ∈ (0, 1) defines the upper bound on the
violation probability of the chance constraints (27) – (32).
The function W (ζ) denotes the system state as a function of
the forecast errors. The chance constrained AC-OPF problem
(25) – (32) is an infinite-dimensional problem optimizing over
W (ζ) which is a function of a continuous uncertain variable
ζ [21]. This renders the problem intractable, which makes it
necessary to identify a suitable approximation for W (ζ) [24].
In the following, an approximation of an explicit dependence
of W (ζ) on the forecast errors is presented.
III. P IECEWISE A FFINE P OLICY
We present a formulation of the chance constraints using a
piecewise affine policy, which approximates the system change
as a linear function of the forecast errors. This allows us
to include corrective control policies for active and reactive
power, and voltages. We propose a tractable formulation for
two types of uncertainty sets. First, using an approach based
on randomized and robust optimization, and making no prior
assumption on the underlying probability distributions, we
determine a rectangular uncertainty set. For that, it is sufficient
to enforce the chance constraints at its vertices. Second,
assuming a Gaussian distribution of the forecast errors, we
can provide an analytical reformulation of the linear chance
constraints and a suitable approximation of the semidefinite
chance constraints.
A. Formulation of Chance Constraints
The main idea is to describe the matrix W (ζ) as the sum
of the forecasted system operating state W0 and the change
of the system state Bi due to each forecast error. Similar to
[21], the matrix W (ζ) is approximated using the affine policy
W (ζ) = W0 +

nw
X

ζi Bi

(33)

i=1

where W0 and Bi are matrices modeled as decision variables.
Eq. (33) provides an affine parametrization of the solution
space for the product of real and imaginary part of bus voltages
described by W (ζ). The main advantages of the affine policy
are that it resembles affine corrective control policies and

ACCEPTED FOR PUBLICATION IN IEEE TRANSACTIONS ON POWER SYSTEMS

naturally allows to include these as well. Furthermore, as the
system change depends linearly on the forecast error, in case a
Gaussian distribution is assumed, an analytical reformulation
can be applied as we will show in Section III-E. Inserting (33)
in (27) – (32) yields:
nw
n
X
P P k ≤ Tr{Yk W0 } +
ζi Tr{Yk Bi } ≤ P k
(34)
Qk ≤ Tr{Ȳk W0 } +

i
nw
X

V 2k ≤ Tr{Mk W0 } +

ζi Tr{Ȳk Bi } ≤ Qk

i
nw
X

(35)
2

ζi Tr{Mk Bi } ≤ V k

i

− P lm ≤ Tr{Ylm W0 } +

nw
X

(36)

ζi Tr{Ylm Bi } ≤ P lm (37)

i

 2

Q
−S lm ΞP
lm Ξlm
 ΞP
−1
0 0
lm
Q
Ξlm
0
−1
nw
o
X
ζi Bi  0 ≥ 1 − 
W0 +

(38)

(39)

i

PnW
The terms ΞP
W } + i=1
ζi Tr{Ylm Bi } and
lm := Tr{Y
PnlmW 0
Q
Ξlm := Tr{Ȳlm W0 } + i=1 ζi Tr{Ȳlm Bi } denote the active
and reactive power flow on transmission line (l, m) ∈ L as a
function of the forecast errors. Note that the chance constraints
(34) – (39) are convex and can be classified in two groups: The
constraints (34) – (37) are linear scalar chance constraints and
the constraints (38) – (39) are semidefinite chance constraints.
B. Corrective Control Policies
The affine policy allows to include corrective control policies related to active power, reactive power, and voltage in the
AC-OPF formulation. In this work, the implemented policies
are generator active power control, generator voltage control,
and wind farm reactive power control.
Throughout the transmission system operation, generation
has to match demand and system losses. If an imbalance occurs, automatic generation control (AGC) restores the system
balance [25]. Hence, designated generators in the power grid
will respond to changes in wind power by adjusting their
output as part of secondary frequency control. The generator
participation factors are defined in the vector dG ∈ Rnb . The
term nb denotes the number of buses. The sum of the change
in generator active power set-points
P should compensate the
deviation in wind generation, i. e. k∈G dGk = 1. The wind
vector diW ∈ Rnb for each wind feed-in i in [1, nW ] has a
{−1} entry corresponding to the bus where the i-th wind farm
is located at. The other entries are zero. The line losses of
the AC power grid vary non-linearly with changes in wind
infeeds. To compensate for this change in system losses, we
add a slack variable γi to the generator set-points. This results
in the following constraints on each matrix Bi , bus k ∈ N
and wind feed-in i in [1, nW ]:
Tr{Yk Bi } = dGk (1 + γi ) + diWk

(40)

As a result of (40), it is ensured that each generator compensates the non-linear change in system losses according to its

5

participation factor. To constrain the magnitude of the slack
variable, a penalty term is added to the objective function (25),
where the term µ ≥ 0 is a penalty weight:
min

W, α, γ

X

αk + µ

nw
X

γi

(41)

i

k∈G

This penalty guides the optimization to a physically meaningful solution, i.e. it allows us to obtain rank-1 solution
matrices. The increase in losses due to deviations in wind infeeds is minimized. With this penalized semidefinite AC-OPF
formulation, near-global optimality guarantees can be derived
specifying the maximum distance to the global optimum [20].
The numerical results show that while this penalty is necessary
to obtain zero relaxation gap, in practice the deviation from
the global optimum is very small. This is investigated in detail
in Section V.
In power systems, automatic voltage regulators (AVR) are
installed as part of the control unit of generators. They keep
the voltages at the generator terminals to a value fixed by
the operator or a higher level controller [26]. The voltage setpoint at each generator k ∈ G is changed as a function of the
forecast errors [21] and can be retrieved using:
Vk (ζ)2 =Tr{Mk W0 } +

nw
X

ζi Tr{Mk Bi }

(42)

i=1

According to recent revisions in Grid Codes [27], renewable
generators such as wind farms have to be able to provide or
absorb reactive power up to a certain extent. This
q is2 often
P
specified in terms of a power factor cos φ :=
P 2 +Q2 . In
this paper, we include the reactive power capabilities of the
wind farms in the optimization. Note that these vary depending
on the magnitude of the actual wind infeed. For each k ∈ W
the constraints (23) and (24) are replaced by:
f
Qk := QGk − QDk + τ (PW
+ ζk )
k

(43)

f
τ (PW
k

Qk := QG − QDk −
+ ζk )
(44)
k
q
2
φ
where τ := 1−cos
cos2 φ . Using this procedure, active and reactive power set-points of FACTS devices and HVDC converter
can also be included in the optimization.
C. Piecewise Affine Policy
In Fig. 2 the affine policy for a wind infeed PWi is depicted.
By choosing an affine policy in the form of (33), the maximum
and minimum bounds of the uncertainty set are linearly
connected using the matrix Bi , and the OPF solution at the
bounds can be recovered. As the OPF is a non-linear problem,
the true system variation will likely not coincide with the
linearization. Hence, the affine policy of [21] is not exact at
the operating point W0 , but returns only an estimate W00 , i.e. a
non-physical higher rank solution. To obtain an exact solution
for W0 , i.e. a rank-1 solution, we introduce a modification to
the conventional affine policy by separating the linearization
between the maximum and minimum value into an upper part
Biu and a lower part Bil , and thereby introducing a piecewise
affine policy. Thus, we linearize between the operating point

ACCEPTED FOR PUBLICATION IN IEEE TRANSACTIONS ON POWER SYSTEMS

W00 + ζ i Bi

System change

W0 +

ζ i Biu

W00
Bi
W00 + ζ i Bi

Biu

Bil
W0

6

convex. Hence, it suffices to enforce the chance constraints at
the vertices v ∈ V of the uncertainty set.
The vector ζv ∈ RnW collects the forecast error bounds
for each vertex, i.e. the entries of this vector correspond to
the the deviation of each wind farm for a specific vertex v.
For each vertex, a corresponding slack variable γv is defined.
Based on our experience with the SDP solvers, we introduce
the following more numerically robust formulation:

W0 + ζ i Bil

Wv := W0 +

nW
X

ζ vi B i

(46)

i=1
f
PW
+ ζi
i

f
PW

i

f
PW
+ ζi

PWi

i

Fig. 2. Piecewise affine policy: The linearization between upper and lower
limit is split into two corresponding piecewise linearizations starting from
the exact operating point W0 . The red line indicates the true system behavior
and the dashed lines the approximation which is made with the corresponding
affine policy. This modification allows us to obtain the exact rank-1 solution
W0 , not the higher-rank approximation W00 .

and the maximum and minimum value of the uncertainty set,
respectively. We extend the work of [21], by ensuring that the
obtained solution is exact at the operating point. An additional
benefit of our approach is that we get a closer approximation
of the true system behavior, while the obtained control policies
are piecewise linear.
D. Tractable Formulation for Rectangular Uncertainty Set
In this section, we provide a tractable formulation of the
chance constraints for a rectangular uncertainty set. The proposed procedure is a combination of robust and randomized
optimization from [22] and which is applied to chance constrained AC-OPF in [21]. A scenario-based method, which
does not make any assumption on the underlying distribution
of the forecast errors, is used to compute the bounds of
the uncertainty set. Two parameters need to be specified,
 ∈ (0, 1) is the allowable violation probability of the chance
constraints and β ∈ (0, 1) a confidence parameter. Then, the
minimum volume hyper-rectangular set is computed, which
with probability 1 − β contains 1 −  of the probability mass.
According to [21], it is necessary to include at least the
following number of scenarios Ns to specify the uncertainty
set:
1
e
1
Ns ≥
(ln + 2nW − 1)
(45)
1−e−1 β
The term e is Euler’s number. The minimum and maximum
bounds on the forecast errors ζi ∈ [ζ i , ζ i ] are retrieved by
a simple sorting operation among the Ns scenarios and the
vertices, i.e. the corner points, of the rectangular uncertainty
set can be defined.
To obtain a tractable formulation of the chance constraints,
the following result from robust optimization is used: If
the constraint functions are linear, monotone or convex with
respect to the uncertain variables, then the system variables
will only take the maximum values at the vertices of the
uncertainty set [22]. The chance constraints (34) – (37) are
linear and the semidefinite chance constraints (38), (39) are

The matrix Wv denotes the power flow solution at the corresponding vertex v. The active and reactive power limits for
each bus k ∈ N and vertex v ∈ V can be written as:
v

f
Qk := QGk − QDk + τ (PW
+ ζvk )
k

(47)

f
QDk − τ (PW
+ ζvk )
k
f
PDk + PWk + ζvk
f
PDk + PW
+ ζvk
k

(48)

Qvk
v
Pk
P vk

:= QG −
k

:= P Gk −
:= P Gk −

(49)
(50)

We provide a tractable formulation of chance constraints (34)
– (39) for each vertex v ∈ V, bus k ∈ N and line (l, m) ∈ L:
v

P vk ≤ Tr{Yk Wv } ≤ P k
Qvk
V 2k

≤
≤

(51)

v
Tr{Ȳk Wv } ≤ Qk
2
Tr{Mk Wv } ≤ V k

− P lm ≤ Tr{Ylm Wv } ≤ P lm
"

−(S lm )2 Tr{Ylm Wv } Tr{Ȳlm Wv }
Tr{Ylm Wv }
−1
0
Tr{Ȳlm Wv }
0
−1

(52)
(53)
(54)
#
0

Wv  0
nW
X
Tr{Yk (Wv − W0 )} =
ζvi (dGk (1 + γv ) + diWk )

(55)
(56)
(57)

i=1

The constraint (57) links the forecasted system state to each
of the vertices. To enforce the semidefinite chance constraint
(39) for the uncertainty set, it suffices that Wv is positive
semidefinite at the vertices of the uncertainty set, i. e. (56)
is fulfilled. For illustrative purposes, in Fig. 3 a rectangular
uncertainty set is depicted for two uncertain wind infeeds PW1
and PW2 . The resulting optimization problem for a rectangular
uncertainty set of dimension nw minimizes the objective (41)
subject to constraints (26) and (51) – (57). Note that the
proposed formulation holds for an arbitrary high-dimensional
rectangular uncertainty set.
E. Tractable Formulation for Gaussian Uncertainty Set
In the following, it is assumed that the forecast errors ζ are
random variables following a Gaussian distribution with zero
mean and covariance matrix Λ. Assuming a Gaussian distribution can be helpful when there is insufficient amount of data at
hand, as it can provide a suitable approximation of the power
system operation under uncertainty. At the same time, through
the covariance matrix, geographical correlations between wind
farms, solar PV plants, or other types of uncertainty can be
captured. We give a direct tractable formulation of the chance

ACCEPTED FOR PUBLICATION IN IEEE TRANSACTIONS ON POWER SYSTEMS

f
PW
2

PW2
W4

7

PW2

W1

+ζ 2

W0 + κ2 B̃2u
(I)

f
PW
2

f
PW
2

W0

(II)

(III)
W0 − κ1 B̃1l
W3
f
PW
+ ζ1
1

W0 − κ2 B̃2l

W2
f
PW
1

f
PW
+ ζ1
1

PW1
f
PW
1

Fig. 3. Rectangular uncertainty set derived from a scenario-based method
displayed for two wind farms. It is sufficient to enforce the chance constraints
at the vertices of the uncertainty set. The vertices are denoted with circles.

constrained AC-OPF, as the work in [3] presented for the
chance constrained DC-OPF.
For a defined confidence interval 1 − , the uncertainty set
for a Gaussian distribution of the forecast errors is an ellipsoid.
First, the direction of linearization of the B matrices is rotated
to correspond to the ellipsoid axes which are described by the
eigenvectors ηi of the covariance matrix. The eigenvalues λi
describe the squared dimension of the ellipsoid in the direction
of its axes. Similar to the rectangular uncertainty set, we
introduce the following auxiliary variables for each ellipsoid
axis i in [1, nW ] and bus k ∈ W:
p
d˜G := dG ||ηi ||, d˜iWk := ηi , ζ˜i := λi
(58)
With B̃i we denote the matrices of the affine policy rotated in
the direction of the ellipsoid axes and (40) has to hold:
Tr{Yk B̃i } = d˜Gk (1 + γi ) + d˜iWk

(59)

Second, we use theoretical results on chance constraints from
the work in [28], which presents the theory for an analytical
reformulation of linear scalar chance constraints. To apply
the reformulation, we approximate the joint probability of
the chance constraint violation (34)–(39) with the violation
probability of each individual chance constraint, which is conservative. Applying the reformulation to the chance constraints
(34) – (37) yields for each bus k ∈ N and line (l, m) ∈ L:
v
u nw
uX
P k ≤ Tr{Yk W0 } ± t
(60)
κ2i Tr{Yk B̃i }2 ≤ P k
i

v
u nw
uX
Qk ≤ Tr{Ȳk W0 } ± t
κ2i Tr{Ȳk B̃i }2 ≤ Qk

(61)

i

v
u nw
uX
2
V 2k ≤ Tr{Mk W0 } ± t
κ2i Tr{Mk B̃i }2 ≤ V k

v
u nw
uX
≤ Tr{Ylm W0 } ± t
κ2i Tr{Ylm B̃i }2 ≤ P lm (63)
i

PW1

Fig. 4. Uncertainty set resulting from a Gaussian distribution of the forecast
errors considering correlation. The directions of approximation for the affine
policy are rotated corresponding to the eigenvectors of the covariance matrix.
The circles denote the points for which the definite chance constraint is
enforced. As a result, it holds for the whole dotted rectangular shape. The
indices (I) – (IV) denote the four quadrants of the uncertainty set for each of
which the complete set of constraints (55), (59) and (60) – (64) is included.

The term κi := Φ−1 (1 − )ζ˜i is introduced, where Φ−1
denotes the inverse Gaussian function. The chance constraint
(39) is P
a linear matrix inequality which ensures that the matrix
n
W0 + i w ζ̃i B̃i is positive semidefinite inside a confidence
interval 1 − . An analytical reformulation of this type of
constraint is not known [28].
Pn As a safe approximation, it
suffices to enforce that W0 + i w ζ̃i B̃i is positive semidefinite
at maximum corresponding deviations ±κi to ensure that (32)
is fulfilled. We include the following semidefinite constraints
for each ellipsoid axis i ∈ [1, nW ]:
W0 ± κi B̃i  0

(64)

This results in (39) holding for the outer rectangular approximation of the ellipsoid uncertainty set. The semidefinite
chance constraint on the apparent branch power flow can be
conservatively approximated by enforcing it for the smallest
rectangular set enclosing the ellipsoid, i.e. by including the
constraint (55) in the optimization.
The assumption of a multivariate Gaussian distribution of
the forecast errors leads to an uncertainty set which in two
dimensions can be described as an ellipse. For the case of
two wind farms with uncertain infeeds PW1 and PW2 this
configuration is depicted in Fig. 4. Incorporating the results on
the modification of the affine policy presented in Section III-C,
we add the constraints (60) – (63) not for B̃i but for both
B̃iu and B̃il and each of their combinations, splitting the
uncertainty set into four quadrants (I) – (IV) as depicted
in Fig. 4. The resulting optimization problem corresponds to
minimizing objective function (41) subject to constraints (26),
(55), (59) and (60) – (64) for each quadrant of the ellipsoid.

(62)

i

−P lm

(IV)
W0

f
PW
2

+ζ 2

W0 + κ1 B̃1u

IV. L INEARIZATION USING PTDF S
In the following, an alternative approach is presented which
is used as benchmark for comparison with the rest of the
approaches presented in this paper. To describe the system

ACCEPTED FOR PUBLICATION IN IEEE TRANSACTIONS ON POWER SYSTEMS

change as a function of the forecast errors, in this section we
introduce a linear approximation based on DC power flow.
This linear approximation uses the so-called power transfer
distribution factors (PTDFs) to estimate the change in line
loading due to a change in active power injections. This
approach has been used in the works in [3] and [29] in the
context of DC- and AC-OPF, respectively.
The PTDFs use the DC power flow representation, i. e.
assuming that the voltage magnitudes of all buses are equal
to 1 p.u. and the resistances of branches are neglected.
Hence, line losses are neglected and the generator participation
factors are defined without including the slack term γ. As we
assume constant voltage magnitudes, the semidefinite (32), the
voltage (29) and the reactive power (28) chance constraints
are dropped and the focus is on approximating the chance
constraints for the active power bus injection and active power
branch flow, Eqs. (27) and (30). The admittance matrix BDC is
constructed using only the line reactances xlm . The resulting
matrix is singular. Thus, one column and the corresponding
row are removed to obtain B̃DC . The vectors dG and diW
collect the generator participation factors and wind injections,
and d˜G and d˜iW denote the corresponding vectors with the first
entry removed. The PTDF for each line (l, m) ∈ L is defined
as follows:
−1
1
B̃DC
PTDFlm = (el − em )T xlm

(65)

The PTDFs provide an approximate linear relation between a
change in bus power injections and the change of the active
power flow over a transmission line. Assuming the maximum
and minimum bounds of the forecast errors are described by a
rectangular uncertainty set with vertices ζv from the previously
described scenario-based approach, we formulate a tractable
approximation of (27) and (30) for each bus k ∈ N , line
(l, m) ∈ L and vertex v ∈ V:
P vk

≤ Tr{Yk W0 } +

nW
X

ζvi (dGk +

diWk )

≤

v
Pk

(66)

i

−P lm ≤ Tr{Ylm W0 }+
nW
X
PTDFlm ζvi (d˜G + d˜iW ) ≤ P lm (67)
i

Assuming the forecast errors follow a Gaussian distribution
with zero mean and co-variance matrix Λ, we formulate a
tractable approximation of (27) and (30) for each bus k ∈ N
and line (l, m) ∈ L:
q
P k ≤ Tr{Yk W0 } ± Φ−1 (1 − ) d2Gk 1T Λ1 ≤ P k (68)
√
−P lm ≤ Tr{Ylm W0 } ± Φ−1 (1 − ) ΨT ΛΨ ≤ P lm (69)
The term 1 ∈ RnW denotes the vectors of ones. The vector
Ψ ∈ RnW contains for each wind feedin i ∈ [1, nW ] the
approximated change in line loading:
Ψi = PTDFlm (d˜G + d˜iW )

(70)

V. S IMULATION A ND R ESULTS
In this section, we first describe the simulation setup.
Subsequently, using the IEEE 24 bus test case, we investigate

8

the relaxation gap of the obtained solution matrices as a
function of the penalty weight. Detailed results on the IEEE
118 bus test case using realistic forecast data are provided
and our proposed approaches are compared to two alternative
approaches described in the literature.
A. Simulation Setup
The optimization problem is implemented in Julia using the
optimization toolbox JuMP [30] and the SDP solver MOSEK
8 [31]. A small resistance of 10−4 has to be added to each
transformer, which is a condition for obtaining zero relaxation
gap [15]. To investigate whether the relaxation gap of an
obtained solution matrix W is zero, the ratio ρ of the 2nd
to 3rd eigenvalue is computed, a measure proposed by [23].
This value should be around 105 or larger for zero relaxation
gap to hold, which means that the obtained solution matrix
is rank-2. The respective rank-1 solution can be retrieved by
following the procedure described in [23]. According to [15],
the obtained solution is then a feasible solution to the original
non-linear AC-OPF problem.
The work in [20] proposes the use of the following measure
to evaluate the degree of the near-global optimality of a
penalized SDP relaxation. Let f˜1 (x) be the generation cost
of the convex OPF without a penalty term and f˜2 (x) the
generation cost of the convex OPF with a penalty weight
sufficiently high to obtain rank-1 solution matrices. Then,
the near-global optimality can be assessed by computing the
˜
· 100%. The closer this parameter is to
parameter δopt := ff˜1 (x)
2 (x)
100%, the closer the solution is to the global optimum. Note
that this distance is an upper bound to the distance from global
optimality.
B. Investigating the Relaxation Gap
This section investigates the relaxation gap of the obtained
matrices. With relaxation gap, we refer to the gap between the
SDP relaxation and a non-linear chance constrained AC-OPF
which uses the affine policy to parametrize the solution space.
The IEEE 24 bus system with parameters specified in [32]
is used. The allowable violation probability is selected to be
 = 5%. Two wind farms with a forecasted infeed of 50 MW
and 150 MW and a maximum power of 150 MW and 400 MW
are introduced at buses 8 and 24, respectively. For illustrative
purposes, the forecast error for the rectangular uncertainty
is assumed to be bounded within ±50% of the forecasted
value with 95% probability. For the Gaussian uncertainty set,
a standard deviation of 25% of the forecasted value and
no correlation between both wind farms is assumed. Each
generator adjusts its active power proportional to its maximum
active power to react to deviations in wind power output.
For the rectangular uncertainty set, Fig. 5 shows the eigenvalue ratios ρ of the matrices W0 − W4 as a function of the
penalty weight µ. A certain minimum value for the weight
µ = 175 is necessary to obtain solution matrices with rank-1,
i.e. eigenvalue ratio ρ higher than 105 , at the operating state
W0 and the four vertices of the rectangular uncertainty set
W1 − W4 . The near-global optimality at µ = 175 for this test
case evaluates to a tight upper bound of 99.74%. If the penalty

105
0

100

200

300

400

4
5.9 ·10

Gen.
cost

5.85
0
200
0
-200

ρ(W0 )
ρ(W1 )
ρ(W2 )
ρ(W3 )
ρ(W4 )

100

200

300

1
0.8
0.6
0.4
0.2
0
1
0.8
0.6
0.4
0.2
0

1

2

1

2

400
Penalty
term

0

PW1 (p.u.)

1010

9

PW2 (p.u.)

Cost ($/h) Cost ($/h)

Eigenvalue ratio ρ

ACCEPTED FOR PUBLICATION IN IEEE TRANSACTIONS ON POWER SYSTEMS

100
200
300
400
Weight for power loss penalty µ (p.u.)

3

3
Time (h)

4

4

5

Bounds
[ζ 1 , ζ 1 ]
Forecast
f
PW
1

5

Bounds
[ζ 2 , ζ 2 ]
Forecast
f
PW
2

Fig. 7. Forecast data from hour 1 to hour 5. The bounds correspond to
the minimum and maximum values from the Ns sampled scenarios for the
rectangular uncertainty set.

Eigenvalue ratio ρ

Fig. 5. Eigenvalue ratios ρ, generation cost and penalty term as a function of
the power loss penalty weight µ for a IEEE 24 bus test case with two wind
farms and a rectangular uncertainty set.

Rectangular
uncertainty set
Gaussian
uncertainty set
Samples
Forecast
f
f
PW
, PW
1
2

0.8
0.6
0.4

1010

ρ(W0 )
ρ(W0 + κ1 B̃1u )

105

0.4

0.6
0.8
PW1 (p.u.)

1

ρ(W0 − κ1 B̃1l )
0

50

100

150

ρ(W0 + κ2 B̃2u )
ρ(W0 − κ2 B̃2l )

Cost ($/h) Cost ($/h)

PW2 (p.u.)

1

Fig. 8. Comparison of rectangular and Gaussian uncertainty set obtained from
the Ns scenarios sampled for hour 4.

58,227
58,222

Gen.
cost
0

50

100

150

60
30

Penalty
term
0
50
100
150
Weight for power loss penalty µ (p.u.)

Fig. 6. Eigenvalue ratios ρ, generation cost and penalty term as a function of
the power loss penalty weight µ for a IEEE 24 bus test case with two wind
farms and a Gaussian uncertainty set.

weight is increased beyond µ = 375 a higher rank solution is
obtained for the forecasted system state.
A similar observation can be made if a Gaussian distribution
is assumed for the forecast errors. Fig. 6 shows the eigenvalue
ratios ρ as a function of the penalty weight µ for the Gaussian
uncertainty set. A certain minimum value for the weight µ =
10 is necessary to obtain solution matrices with rank-1 at the
operating state W0 and the four end-point of the ellipsoid axes.
The generation cost is almost flat with respect to increasing
penalty weight and the near-global optimality at µ = 10 for
this test case evaluates to an upper bound larger than 99.99%.
As it is also observed, the necessary magnitude of the penalty
weight µ to obtain rank-1 solution matrices depends on the
test case and configuration.
C. IEEE 118 Bus Test Case
In this section, our proposed approaches using the affine policy and PTDFs are compared with two alternative approaches
described in the literature [5], [12]. We use the IEEE 118 bus

test case with realistic forecast data for the wind farms, and
Monte Carlo simulations to evaluate the constraint violations.
1) Simulation Setup: We use the IEEE 118 bus specifications from [33] with the following modifications: The bus
voltage limits are set to 0.94 p.u. and 1.06 p.u. As the upper
branch flow limits are specified in MW, the active line flow
limit is considered for branch flows. The line flow limits are
decreased by 30% and the load is increased by 30% to obtain a
more constrained system. Two wind farms with a rated power
of 300 MW and 600 MW are placed at buses 5 and 64. The
rated wind power corresponds to 24.1% of total load demand.
Realistic day-ahead wind forecast scenarios from [34] and [35]
are used for both wind farms. To create the scenarios, the
methodology described in [34] is used. The forecasts are based
on wind power measurements in the Western Denmark area
from 15 different control zones collected by the Danish transmission system operator Energinet. We select control zone 1
to correspond to the wind farm at bus 5 and zone 7 to the wind
farm at bus 64. We allow a constraint violation of  = 5% for
all considered approaches. In order to construct the rectangular
uncertainty set, the confidence parameter β = 10−3 is selected.
Then, a minimum of 314 scenarios are required according to
(45). The forecast is computed as mean value of the scenarios.
For the Gaussian uncertainty set, we compute the co-variance
matrix based on these 314 scenarios. Fig. 7 shows the forecast
data from hour 1 to hour 5 with the upper and lower bounds
specified by the maximum and minimum scenario values,
respectively. In Fig. 8 the rectangular and Gaussian uncertainty
set for hour 4 are shown.
In the following, the parameters for the corrective control

ACCEPTED FOR PUBLICATION IN IEEE TRANSACTIONS ON POWER SYSTEMS

policies are specified. A participation factor of 0.25 is
selected for the generators at buses i = {12, 26, 54, 61}, i.e.
dGi = 0.25. Wind farms have a reactive power capability of
0.95 inductive to 0.95 capacitive according to recent Grid
Codes [27]. The approaches using PTDFs assign a fixed
power factor cos φ to each wind farm. The affine policy
includes a generator voltage and wind farm reactive power
corrective control, assigning an updated set-point to generators
and wind farms based on the actual realization of the forecast
errors. To facilitate comparability, we use the same scenarios
for all approaches to compute the respective uncertainty sets.
We evaluate the constraint violations using Monte Carlo
simulations with 10’000 scenarios and MATPOWER AC
power flows [32]. We enable the enforcement of generator
reactive power limits in the power flow, i.e. P V buses
are converted to P Q buses once the limits are reached, as
otherwise high nonphysical overloading of the limits can
occur [36]. Furthermore, we distribute the loss mismatch
from the active generator set-points among the generators
according to their participation factors and rerun the power
flow to mimic the response of automatic generation control
(AGC).
2) Numerical Comparison to Alternative Approaches: In
the following, the main modeling assumptions of the respective
approaches and the type of chance constraints they include
are outlined. All approaches considering chance constraints
include corrective control of the active generator set-points.
•

•

Chance constrained DC-OPF [5] (DC-OPF): A robust formulation based on DC-OPF includes chance constraints
on active generator power and active branch flow. Interval
bounds on the forecast errors are assumed. Hence, we use
the scenarios to compute the interval bounds. A power
factor of 1 is assumed for wind farms.
Iterative chance constrained AC-OPF [12] (Iterative): At
each iteration the Jacobian is computed and the uncertainty margins resulting from the chance constraints
are updated until convergence is reached. The forecast
errors are assumed to follow a Gaussian distribution. The
covariance matrix constructed from the Ns scenarios is
used. Chance constraints on active and reactive generator
limits, voltage magnitudes and apparent line flows are
included in the formulation. A power factor of 1 is
assumed for wind farms, as no reactive power corrective
control is included in [12].

These two approaches are compared to the following approaches based on the formulations presented in this work:
•
•

•

AC-OPF with convex relaxations but without chance
constraints (AC-OPF) [15]
Chance constrained AC-OPF with convex relaxations,
using an affine policy for a Gaussian uncertainty set
(AP (Gauss)) including corrective control for wind farms,
generator voltages and generator active power.
Chance constrained AC-OPF with convex relaxations,
using an affine policy for a rectangular uncertainty set
(AP (Rect)) including corrective control for wind farms
and generator voltages and generator active power.

10

TABLE II
C OST OF UNCERTAINTY: G ENERATION COST IN RELATION TO AC-OPF
WITH CONVEX RELAXATIONS WITHOUT CONSIDERATION OF
UNCERTAINTY

Time step (h)

1

2

3

4

5

AP (Rect) (%)

0.774

0.740

0.755

0.921

1.562

PTDF (Rect) (%)

0.785

0.748

0.767

0.931

1.588

DC-OPF [5] (%)

-2.782

-2.826

-2.824

-2.673

-2.165

AP (Gauss) (%)

0.515

0.461

0.467

0.489

0.512

PTDF (Gauss) (%)

0.523

0.468

0.477

0.497

0.520

Iterative [12] (%)

0.519

0.465

0.473

0.494

0.516

•

•

Chance constrained AC-OPF with convex relaxations,
using PTDFs (PTDF (Gauss)) for a Gaussian uncertainty
set.
Chance constrained AC-OPF with convex relaxations,
using PTDFs (PTDF (Rect)) for a rectangular uncertainty
set.

In Table II the cost of uncertainty for the different approaches
and considered time steps are shown. The cost of uncertainty
represents the additional cost incurred by considering the
stochastic variables, and is defined as the difference between
the solution of the chance constrained and a baseline. In
this paper, the AC-OPF with convex relaxations but without
considering uncertainty is assumed as the baseline cost. From
Table II, we make the following observations. First, the DCOPF (with chance constraints) leads to a cost reduction, as no
losses are considered compared to the AC-OPF. Second, the
approaches stemming from robust optimization lead to a cost
increase of approximately 0.8% for time step 1 compared to
an increase of approximately 0.5% for the same time step
for the approaches assuming a Gaussian distribution. This
shows that the Gaussian uncertainty set is less conservative
as indicated in Fig. 8. For the rectangular uncertainty set,
the affine policy reduces the cost compared to the approach
using PTDFs. Comparing the approaches for the Gaussian
uncertainty set, again the affine policy results to the lowest
cost of uncertainty compared to the approach using PTDFs
and the iterative chance constrained AC-OPF. The reason for
that is that the affine policy includes corrective control for
voltages and both active and reactive power.
In Table III the violation probability of the chance constraints on active power, voltages, and active branch flows
are shown. Monte Carlo simulations using 10’000 scenarios
with MATPOWER AC power flows are conducted. A minimum violation limit of 10−3 p.u. for active generator limits
and 0.1% for voltage and line flow limits is considered to
exclude numerical errors. In all considered time steps, the ACOPF without consideration of uncertainty leads to insecure
instances and violates constraints on line and generator limits
on active power.
First, investigating the robust approaches using the rectangular uncertainty set the following observations can be made: The
robust DC-OPF formulation in [5] leads to insecure instances
for all time steps and violates both voltage and generator active

ACCEPTED FOR PUBLICATION IN IEEE TRANSACTIONS ON POWER SYSTEMS

TABLE III
V IOLATION PROBABILITY OF THE CHANCE CONSTRAINTS ON BUS
VOLTAGE , ACTIVE BRANCH FLOW AND ACTIVE GENERATOR POWER FOR
THE FORECAST DATA . M ONTE C ARLO SIMULATIONS USING 10’000
SCENARIOS WITH MATPOWER AC POWER FLOWS ARE CONDUCTED .
I NSECURE INSTANCES ARE MARKED IN BOLD .

11

TABLE IV
C OMPARISON OF VIOLATION PROBABILITY OF THE CHANCE
CONSTRAINTS ON BUS VOLTAGE , ACTIVE BRANCH FLOW AND ACTIVE
GENERATOR POWER FOR AFFINE POLICY AND ITERATIVE AC-OPF WITH
10’000 SAMPLES FROM A G AUSSIAN DISTRIBUTION .
Time step (h)

Time step (h)

1

2

3

4

1

2

3

4

5

5
Bus voltage

Bus voltage
AC-OPF (%)

0.0

0.1

0.2

1.0

4.3

DC-OPF [5] (%)

100.0

100.0

100.0

100.0

100.0

PTDF (Rect) (%)

19.5

20.9

14.6

13.0

12.9

AP (Rect) (%)

0.0

0.0

0.0

0.0

0.0

PTDF (Gauss) (%)

20.0

21.2

15.0

13.5

13.0

AP (Gauss) (%)

0.7

1.9

4.3

7.2

7.6

Iterative [12] (%)

0.0

0.0

0.0

0.0

0.0

Active power line limit

AP (Gauss) (%)

0.1

0.4

0.5

0.8

0.7

Iterative [12] (%)

0.0

0.0

0.0

0.0

0.0

Active power line limit
AP (Gauss) (%)

2.2

1.7

2.1

2.1

2.1

Iterative [12] (%)

1.5

1.7

1.7

2.1

2.1

Active generator limit
AP (Gauss) (%)

2.4

2.7

2.7

2.6

2.5

Iterative [12] (%)

2.7

2.6

2.3

2.4

2.1

AC-OPF (%)

17.7

18.8

14.9

32.5

46.5

DC-OPF [5] (%)

0.0

0.0

0.0

2.8

0.0

PTDF (Rect) (%)

0.0

0.0

0.0

0.0

0.0

TABLE V
P OWER LOSS PENALTY WEIGHT AND NEAR - GLOBAL OPTIMALITY
GUARANTEES FOR IEEE 118 BUS TEST CASE FOR ALL CONSIDERED TIME

AP (Rect) (%)

0.0

0.0

0.0

0.0

0.0

STEPS

PTDF (Gauss) (%)

4.6

11.1

13.1

9.3

7.8

AP (Gauss) (%)

4.6

3.7

0.9

2.6

5.8

Iterative [12] (%)

1.6

2.0

3.6

4.2

5.3

Active generator limit
AC-OPF (%)

46.4

48.8

45.9

45.5

DC-OPF [5] (%)

34.3

38.6

30.1

14.5

2.0

PTDF (Rect)

0.0

0.0

0.0

0.0

0.0

AP (Rect) (%)

0.0

0.0

0.0

0.0

0.0

PTDF (Gauss) (%)

2.6

3.6

2.8

3.1

5.5

AP (Gauss) (%)

0.0

0.2

0.4

2.3

0.7

Iterative [12] (%)

2.9

4.1

3.0

3.3

5.7

Penalty weight µ
(p.u.)

Near-global optimality
guarantee δopt (%)

AP (Rect)

100

≥ 99.99

AP (Gauss)

100

≥ 99.99

40.9

power constraints. The AC-OPF approach using PTDFs for
the chance constraints reduces the voltage violations but does
not comply with the 5% confidence interval. The AC-OPF
using the affine policy complies with the chance constraints
for all time steps while slightly decreasing the generation
cost compared to the approach using PTDFs. As the scenario
based method is conservative, there are nearly zero violations
occurring for the considered 10’000 samples for the approach
using the affine policy.
Second, we compare the different approaches which assume
a Gaussian distribution of the forecast errors. The affine policy
improves upon the approach using PTDFs and results to a
secure operation for time steps 1 to 3. For time steps 4 and 5
we observe a slight violation of the active power line and bus
voltage limit. This is due to the fact that we do not sample out
of a Gaussian distribution but out of a set of realistic forecast
scenarios, that apparently are not Gaussian distributed. The
iterative approach results to a secure operation for time steps
1 to 4 and slightly violates the active generator and branch
flow limit in time step 5.

In order to verify if these violations occur due to the mismatch between actual distribution and the assumed Gaussian
we repeat the 10’000 scenario evaluations for both affine
policy and the iterative chance constrained AC-OPF from [12].
We sample from the Gaussian distribution assumed for the
uncertainty set. The results are shown in Table IV. For all
5 time steps, both approaches comply with the 5% violation
probability. Hence, the occurring violations in Table III stem
from the mismatch between Gaussian distribution and actual
probability distribution. As shown in Table II, the affine policy
results in a slightly lower generation cost than the iterative ACOPF, as it includes corrective control policies. This leads us
to the following conclusions. First, that if the forecast errors
do follow a normal distribution both approaches demonstrate
good performance and do not exceed the violation limit. If the
data are not normally distributed, as is the case for the results
shown in Table III, none of the two methods can guarantee
that the violation probability will be below . The differences
in performance in that case are, as it would be expected, dataand system-specific. However, independent from the fact if
the underlying probability distribution is Gaussian or not, one
difference that remains is that the approach proposed in this
paper is more rigorous, since it provides guarantees regarding
the global optimality of the obtained solution and allows to
include corrective control policies related to reactive power
and voltage.
Table V lists the penalty weights and obtained near-global
optimality guarantees for the 5 time steps. Note that it is

ACCEPTED FOR PUBLICATION IN IEEE TRANSACTIONS ON POWER SYSTEMS

TABLE VI
S OLVING TIME FOR IEEE 118 BUS TEST CASE
AP (Rect)

AP (Gauss)

PTDF (Rect)

30 sec

10 min

15 sec

PTDF (Gauss)

DC-OPF [5]

Iterative [12]

15 sec

≤ 1 sec

4 sec

sufficient to define for both uncertainty sets a penalty weight of
µ = 100 p.u. to obtain zero relaxation gap, i.e. rank-1 solution
matrices and a near-global optimality guarantee of larger than
99.99%. This means that the maximum deviation from the
global optimum is smaller than 0.01% of the objective value.
Table VI lists the computational time of the different approaches. The optimization problems are solved on a desktop
computer with an Intel Xeon CPU E5-1650 v3 @ 3.5 GHz and
32 GB RAM. For all optimization problems except the iterative
approach, MOSEK V8 [31] is used. The iterative approach
utilizes the MATPOWER AC-OPF. The DC-OPF formulation
is the fastest, as the optimization problem is a linear program.
The computational time increases with increasing constraint
complexity. The SOC constraints in the formulation for the
Gaussian uncertainty set are computationally the most challenging. We observe that the iterative approach, despite the
need for computing a number of iterations, converges faster
than all approaches that utilize convex relaxations and an
SDP solver. Current trends expect the need of more rigorous
optimal power flow approaches in the future, that e.g. can
guarantee a global minimum. In that case the need for further
research to improve both the optimization solvers and the
convex formulations of the AC-OPF problem is apparent.
Possible directions to increase the computational speed of the
proposed approaches are the chordal decomposition technique,
outlined in [23], and distributed optimization techniques, e.g.
the alternating direction method of multipliers (ADMM) for
sparse semidefinite problems in [37]. The chordal decomposition technique can be applied to reduce the computational
burden of the semidefinite constraints on W0 and Bi (18). As
shown in [38] a speed-up by several orders of magnitude can
then be expected for large systems.
VI. C ONCLUSIONS
In this work, a convex formulation for a chance constrained
AC-OPF is presented which is able to provide near-global
optimality guarantees. The OPF formulation considers chance
constraints for all relevant variables, and has an explicit
representation of corrective control policies. Two tractable
formulations are proposed: First, a scenario-based method
is applied in combination with robust optimization. Second,
assuming a Gaussian distribution of forecast errors, we provide
an analytical reformulation of the chance constraints. Detailed
case studies on the IEEE 24 and 118 bus test systems are
presented. For the latter, we used realistic forecast data and
Monte Carlo simulations to evaluate constraint violations.
Compared to a chance constrained DC-OPF formulation,
we find that the formulations proposed in this paper are

12

more accurate and significantly decrease constraint violations.
Compared with iterative non-convex AC-OPF formulations,
both our piece-wise affine control policy and the iterative
AC-OPF do not exceed the constraint violation limit for
the Gaussian uncertainty set. Most importantly, our proposed
approach obtains tight near-global optimality guarantees which
ensure that the distance to the global optimum is smaller than
0.01% of the objective value. In our future work, besides
investigating chordal decomposition techniques, we include
security constraints in the proposed formulation by defining
a matrix W s (ζ) for each outage s of a generation unit or
transmission line.
ACKNOWLEDGMENT
The authors would like to thank Pierre Pinson for sharing
the forecast data, Line Roald for providing an updated version
of the code from [12], and Martin S. Andersen and Daniel K.
Molzahn for fruitful discussions.
R EFERENCES
[1] A. Nemirovski and A. Shapiro, “Convex approximations of chance
constrained programs,” SIAM Journal on Optimization, vol. 17, no. 4,
pp. 969–996, 2006.
[2] D. Bienstock, M. Chertkov, and S. Harnett, “Chance-constrained optimal
power flow: Risk-aware network control under uncertainty,” SIAM
Review, vol. 56, no. 3, pp. 461–495, 2014.
[3] L. Roald, F. Oldewurtel, T. Krause, and G. Andersson, “Analytical reformulation of security constrained optimal power flow with probabilistic
constraints,” in IEEE PowerTech, Grenoble, France, 2012.
[4] M. Lubin, Y. Dvorkin, and S. Backhaus, “A robust approach to chance
constrained optimal power flow with renewable generation,” IEEE
Transactions on Power Systems, vol. 31, no. 5, pp. 3840 – 3849, 2016.
[5] R. A. Jabr, S. Karaki, and J. A. Korbane, “Robust multi-period OPF with
storage and renewables,” IEEE Transactions on Power Systems, vol. 30,
no. 5, pp. 2790–2799, 2015.
[6] L. Roald, S. Misra, T. Krause, and G. Andersson, “Corrective control to
handle forecast uncertainty: A chance constrained optimal power flow,”
IEEE Transactions on Power Systems, vol. 32, no. 2, pp. 1626–1637,
2017.
[7] S. S. Guggilam, E. Dall’Anese, Y. C. Chen, S. V. Dhople, and G. B.
Giannakis, “Scalable optimization methods for distribution networks
with high PV integration,” IEEE Transactions on Smart Grid, vol. 7,
no. 4, pp. 2061 – 2070, 2016.
[8] K. Baker, E. Dall’Anese, and T. Summers, “Distribution-agnostic
stochastic optimal power flow for distribution grids,” in North American
Power Symposium (NAPS), Denver, US, 2016.
[9] T. Summers, J. Warrington, M. Morari, and J. Lygeros, “Stochastic optimal power flow based on convex approximations of chance constraints,”
in 2014 Power Systems Computation Conference, 2014.
[10] E. D. Anese, K. Baker, and T. Summers, “Chance-constrained AC
optimal power flow for distribution systems with renewables,” IEEE
Transactions on Power Systems, vol. PP, no. 99, pp. 1–1, 2017.
[11] H. Zhang and P. Li, “Chance constrained programming for optimal
power flow under uncertainty,” IEEE Transactions on Power Systems,
vol. 26, no. 4, pp. 2417–2424, 2011.
[12] J. Schmidli, L. Roald, S. Chatzivasileiadis, and G. Andersson, “Stochastic AC optimal power flow with approximate chance-constraints,” in
IEEE Power and Energy Society General Meeting, Boston, US, 2016.
[13] K. Dvijotham and D. K. Molzahn, “Error bounds on the DC power flow
approximation: A convex relaxation approach,” in IEEE 55th Conference
on Decision and Control (CDC), 2016, pp. 2411–2418.
[14] S. V. Dhople, S. S. Guggilam, and Y. C. Chen, “Linear approximations
to AC power flow in rectangular coordinates,” in 53rd Annual Allerton
Conference on Communication, Control, and Computing, 2015, pp. 211–
217.
[15] J. Lavaei and S. H. Low, “Zero duality gap in optimal power flow
problem,” IEEE Transactions on Power Systems, vol. 27, no. 1, pp.
92–107, 2012.

ACCEPTED FOR PUBLICATION IN IEEE TRANSACTIONS ON POWER SYSTEMS

[16] R. A. Jabr, “A conic quadratic format for the load flow equations of
meshed networks,” IEEE Transactions on Power Systems, vol. 22, no. 4,
pp. 2285–2286, 2007.
[17] X. Bai, H. Wei, K. Fujisawa, and Y. Wang, “Semidefinite programming
for optimal power flow problems,” International Journal of Electrical
Power & Energy Systems, vol. 30, no. 6, pp. 383–392, 2008.
[18] M. B. Cane, R. P. O’Neill, and A. Castillo, “History of optimal power
flow and formulations,” Federal Energy Regulatory Commission, 2012.
[19] B. C. Lesieutre, D. K. Molzahn, A. R. Borden, and C. L. DeMarco,
“Examining the limits of the application of semidefinite programming
to power flow problems,” in 49th Annual Allerton Conference on
Communication, Control, and Computing, 2011, pp. 1492–1499.
[20] R. Madani, S. Sojoudi, and J. Lavaei, “Convex relaxation for optimal
power flow problem: Mesh networks,” IEEE Transactions on Power
Systems, vol. 30, no. 1, pp. 199–211, 2015.
[21] M. Vrakopoulou, M. Katsampani, K. Margellos, J. Lygeros, and G. Andersson, “Probabilistic security-constrained AC optimal power flow,” in
IEEE PowerTech, Grenoble, France, 2012.
[22] K. Margellos, P. Goulart, and J. Lygeros, “On the road between robust
optimization and the scenario approach for chance constrained optimization problems,” IEEE Transactions on Automatic Control, vol. 59, no. 8,
pp. 2258–2263, 2014.
[23] D. K. Molzahn, J. T. Holzer, B. C. Lesieutre, and C. L. DeMarco,
“Implementation of a large-scale optimal power flow solver based
on semidefinite programming,” IEEE Transactions on Power Systems,
vol. 28, no. 4, pp. 3987–3998, 2013.
[24] A. Ben-Tal and A. Nemirovski, “Selected topics in robust convex
optimization,” Mathematical Programming, vol. 112, no. 1, pp. 125–
158, 2008.
[25] Ibraheem, P. Kumar, and D. P. Kothari, “Recent philosophies of automatic generation control strategies in power systems,” IEEE Transactions on Power Systems, vol. 20, no. 1, pp. 346–357, 2005.
[26] H. Vu, P. Pruvot, C. Launay, and Y. Harmand, “An improved voltage
control on large-scale power system,” IEEE Transactions on Power
Systems, vol. 11, no. 3, pp. 1295–1303, 1996.
[27] M. Tsili and S. Papathanassiou, “A review of grid code technical
requirements for wind farms,” IET Renewable Power Generation, vol. 3,
no. 3, pp. 308–332, 2009.
[28] A. Nemirovski, “On safe tractable approximations of chance constraints,” European Journal of Operational Research, vol. 219, no. 3,
pp. 707–718, 2012.
[29] S. Chatzivasileiadis, T. Krause, and G. Andersson, “Flexible AC transmission systems (FACTS) and power system security - A valuation
framework,” in IEEE Power and Energy Society General Meeting,
Detroit Michigan, US, 2011.
[30] M. Lubin and I. Dunning, “Computing in operations research using
Julia,” INFORMS Journal on Computing, vol. 27, no. 2, pp. 238–248,
2015.
[31] MOSEK ApS, MOSEK 8.0.0.37, 2016.
[32] R. D. Zimmerman, C. E. Murillo-Sánchez, and R. J. Thomas, “MATPOWER: Steady-state operations, planning, and analysis tools for power
systems research and education,” IEEE Transactions on Power Systems,
vol. 26, no. 1, pp. 12–19, 2011.
[33] “IEEE 118-bus, 54-unit, 24-hour system,” Electrical and Computer
Engineering Department, Illinois Institute of Technology, Tech. Rep.
[Online]. Available: http://motor.ece.iit.edu/data/JEAS IEEE118.doc
[34] P. Pinson, “Wind energy: Forecasting challenges for its operational
management,” Statistical Science, pp. 564–585, 2013.
[35] W. A. Bukhsh, C. Zhang, and P. Pinson, “An integrated multiperiod OPF
model with demand response and renewable generation uncertainty,”
IEEE Transactions on Smart Grid, vol. 7, no. 3, pp. 1495–1503, 2016.
[36] A. E. Efthymiadis and Y. H. Guo, “Generator reactive power limits and
voltage stability,” in Fourth International Conference on Power System
Control and Management (Conf. Publ. No. 421), 1996, pp. 196–199.
[37] R. Madani, A. Kalbat, and J. Lavaei, “ADMM for sparse semidefinite
programming with applications to optimal power flow problem,” in 54th
IEEE Conference on Decision and Control (CDC), 2015, pp. 5932–5939.
[38] R. A. Jabr, “Exploiting sparsity in SDP relaxations of the OPF problem,”
IEEE Transactions on Power Systems, vol. 27, no. 2, pp. 1138–1139,
2012.

13

Andreas Venzke (S’16) received the M.Sc. degree in
Energy Science and Technology from ETH Zurich,
Zurich, Switzerland in 2017. He is currently working towards the Ph.D. degree at the Department
of Electrical Engineering, Technical University of
Denmark (DTU), Kongens Lyngby, Denmark. His
research interests include power system operation
under uncertainty and convex relaxations of optimal
power flow.

Lejla Halilbasic (S’15) received the M.Sc. degree
in Electrical Engineering from the Technical University of Graz, Austria in 2015. She is currently
working towards the Ph.D. degree at the Department
of Electrical Engineering, Technical University of
Denmark (DTU), Denmark. Her research interests
include optimization of power system operation and
its applications to electricity markets.

Uros Markovic (S’16) received the M.Sc. degree in
Electrical Engineering and Information Technology
from ETH Zurich, Zurich, Switzerland in 2016. He
is currently working towards the Ph.D. degree at the
Power Systems Laboratory, ETH Zurich, where he
joined in March 2016. His research interests include
modeling, control and optimization of inverter-based
power system with low rotational inertia.

Gabriela Hug (S’05, M’08, SM’14) was born in
Baden, Switzerland. She received the M.Sc. degree
in electrical engineering in 2004 and the Ph.D.
degree in 2008, both from Swiss Federal Institute of
Technology (ETH), Zurich, Switzerland. After the
Ph.D. degree, she worked in the Special Studies
Group of Hydro One, Toronto, ON, Canada, and
from 2009 to 2015, she was an Assistant Professor
in Carnegie Mellon University, Pittsburgh, PA, USA.
She is currently an Associate Professor in the Power
Systems Laboratory, ETH Zurich. Her research is
dedicated to control and optimization of electric power systems.

Spyros Chatzivasileiadis (S’04, M’14) is an Assistant Professor at the Technical University of Denmark (DTU). Before that he was a postdoctoral
researcher at the Massachusetts Institute of Technology (MIT), USA and at Lawrence Berkeley National
Laboratory, USA. Spyros holds a PhD from ETH
Zurich, Switzerland (2013) and a Diploma in Electrical and Computer Engineering from the National
Technical University of Athens (NTUA), Greece
(2007). In March 2016 he joined the Center of
Electric Power and Energy at DTU. He is currently
working on power system optimization and control of AC and HVDC grids,
including semidefinite relaxations, distributed optimization, and data-driven
stability assessment.

