arXiv:1501.06613v1 [cs.CE] 26 Jan 2015

Adaptive Robust Transmission Network Expansion
Planning using Structural Reliability and
Decomposition Techniques
R. Mı́ngueza,∗, R. Garcı́a-Bertrandb , J. M. Arroyob
a

Independent consultant, C/Honduras 1, 13160, Torralba de Calatrava, Ciudad Real,
Spain
b
Department of Electrical Engineering, University of Castilla-La Mancha, Ciudad Real,
Spain

Abstract
Structural reliability and decomposition techniques have recently proved to
be appropriate tools for solving robust uncertain mixed-integer linear programs using ellipsoidal uncertainty sets. In fact, its computational performance makes this type of problem to be an alternative method in terms of
tractability with respect to robust problems based on cardinality constrained
uncertainty sets. This paper extends the use of these techniques for solving
an adaptive robust optimization (ARO) problem, i.e. the adaptive robust
solution of the transmission network expansion planning for energy systems.
The formulation of this type of problem materializes on a three-level mixedinteger optimization formulation, which based on structural reliability methods, can be solved using an ad-hoc decomposition technique. The method
allows the use of the correlation structure of the uncertain variables involved
by means of their variance-covariance matrix, and besides, it provides a new
interpretation of the robust problem based on quantile optimization. We also
compare results with respect to robust optimization methods that consider
cardinality constrained uncertainty sets. Numerical results on an illustrative example, the IEEE-24 and IEEE 118-bus test systems demonstrate that
the algorithm is comparable in terms of computational performance with
respect to existing robust methods with the additional advantage that the
correlation structure of the uncertain variables involved can be incorporated
∗

Corresponding author: rominsol@gmail.com, tlfn.: 00 34 926810046

Preprint submitted to European Journal of Operational Research

March 13, 2018

straightforwardly.
1. Introduction
Robust optimization techniques were developed to deal with the problem of designing solutions that are immune to data uncertainty by solving
equivalent deterministic problems, i.e. robust counterparts (RC). The first
published work about these techniques is proposed by Soyters [1], which
allowed uncertain parameters to take their worst possible values within a
given interval. However, although using this method the RC of any linear
mathematical programming problem remains linear, Soyters’s solutions are
excessively conservative, and [2, 3, 4] proposes the use of alternative uncertainty sets, the ellipsoidal uncertainty sets based on the Mahalanobis distance. This alternative method allows controlling the grade of conservatism,
and includes straightforwardly the correlation structure of the uncertain parameters involved by means of the variance-covariance matrix. In contrast,
the RC derived from any linear mathematical programming problem using
ellipsoidal uncertainty sets results in a conic quadratic problem. This feature,
which is not problematic for those cases where interior point algorithms are
applicable because they have been proved to be solvable in polynomial time,
is not particularly attractive for solving robust linear discrete optimization
models, as pointed out by [5]. It is precisely [5] the one that proposes the
use of an alternative uncertainty set that allows the RC problem to remain
as linear, even if binary variables are involved. Instead of allowing all uncertain parameters to take their worst possible values within the given interval
such as [1], it only allows a pre-established maximum number of parameters
Γ to do so. This alternative formulation also provides a robust solution in
terms of probability of infeasibility. Recently, [6] proposes an alternative
and decomposable solution technique inspired on structural reliability FirstOrder-Second-Moment (FOSM) methods, which allows reaching the optimal
solution of the RC problem using ellipsoidal uncertainty sets by solving two
kind of problems within an iterative scheme: one master problem of the
same size as the original RC problem, and one subproblem of considerable
lower size for each uncertain constraint. The method proposed by [6] has the
following advantages:
1. It improves tractability with respect to interior point algorithms, specially for large-scale problems and those including binary decision vari2

ables.
2. In addition and due to its relationship with reliability-based structural
techniques, it allows to calculate exact bounds for the probability of
constraint violation in case the probability distributions of uncertain
parameters are completely defined by using first and second moments
(mean and variance-covariance).
3. In terms of tractability, it constitutes an alternative formulation, specially for problems including binary variables, to that proposed by [5]
using cardinality constrained uncertainty sets.
Note that the formulation proposed by [6] was tested and proved to be
successful on mixed-integer linear mathematical programming problems with
hard constraints that must be satisfied for any possible realization of the uncertain data. This type of problems is relevant for optimal design, however,
problems with recourse [7, 8, 9] where the decision-maker must take some
decisions before discovering the actual value of the uncertain parameters and
having the opportunity to take further action once the uncertain parameters have been revealed, are more common on engineering problems and
can not take advantage of the features of the method proposed by [6] in its
present form. This framework is also known as Adaptive Robust Optimization (ARO).
A clear example of the practical importance of robust linear optimization with recourse (ARO) is the Transmission Network Expansion Planning
(TNEP) problem. It aims to resolve the issue of how to expand or reinforce
an existing electricity transmission network so that adequately serves system
loads over a given horizon. The main difficulty of this problem is to take
decisions under the great amount of uncertainty associated with i) demands,
ii) renewable generation, such as wind and solar power plants, and iii) equipment failure. There are several examples of the successful application of
robust optimization in TNEP problems based on cardinality constrained uncertainty sets proposed by [5], such as, [10, 11, 12, 13, 14, 15]. These works
use three-level formulations where first level minimizes the cost of expansion ([13] also minimizes the maximum regret), the decision variables for this
level are those associated with construction or expansion of lines, ii) second
level selects the realization of the uncertain parameters that maximizes the
system’s operating costs within the uncertainty set, the variables related to
this level are the uncertain parameters, i.e., generation capacity and demand,
and iii) in the third level the system operator selects the optimal decision
3

variables to minimize operating costs for given values of first and second level
variables.
This paper proposes a double iterative solution approach to the threelevel mixed-integer optimization problem associated with TNEP. The idea
is to replace cardinality constrained by ellipsoidal uncertainty sets, which
allows incorporating the correlation structure of the uncertain parameters
straightforwardly by means of the variance-covariance matrix. Inner level
problems (subproblems) are solved throughout an iterative method inspired
from structural reliability and decomposition techniques [6]. Note that, due
to the way the uncertainty sets are constructed, the subproblem does not
contain binary variables. Meanwhile, first level problem (master problem)
uses the constraint-and-column generation method (proposed by [16] ) initially used in this type of problems by [13] and [14], and also utilized by
[15].
The proposed method has the following advantages which might make it
attractive for practical use:
1. Computational performance is comparable with the method proposed
by [15] based on cardinality constrained uncertainty sets. Note that
this method proved to be the most efficient method to date for solving
this type of problem.
2. The subproblem can be solved using two nested linear optimization
problems: i) the first one takes the optimal operating decisions for
given realizations of the uncertain parameters involved, and ii) the second one updates the worst possible realization of the uncertain parameters in terms of operating costs. This second problem, under certain
circumstances, have analytical solutions and otherwise can be solved
efficiently by using interior point algorithms.
3. The approach proposed in this paper requires the definition of the protection level, which establishes approximately the quantile of the uncertain objective function to be optimized. Besides, the proposed method
would not require the definition of bounds for the uncertain parameters,
which might be of interest for certain cases.
The remainder of the paper is structured as follows. Section 2 describes
the adaptive robust formulation of the TNEP problem in compact matrix
form, and discusses the solution approach proposed by [15] inspired in the
works by [12, 13, 14]. In Section 3 the definition of the uncertainty set and
the description of the solution approach are presented. Numerical results
4

for different networks are given in Section 4 and compared with those obtained using the method proposed by [15]. Finally, the paper is concluded in
Section 5.
2. Transmission network expansion planning problem: ARO compact formulation
Note that for clarity in the exposition, we adopt the same nomenclature
used by [15]. Thus, the robust TNEP problem can be written in compact
matrix form as the following three-level mathematical programming problem:
!
Minimize
x

cT x + Maximum
d∈D

Minimum bT y
y ∈ Ω(x, d)

(1)

subject to
cT x ≤ Π
x ∈ {0, 1},

(2)
(3)

where x is the vector of first stage binary variables representing the investment vs no investment in reinforcing or building new lines, c is the investment
cost vector, d is the vector of second stage continuous variables representing
the random or uncertain parameters, i.e. generation capacities and demands,
D is the uncertainty set, b is the vector including operating costs, and y
is the vector of third stage continuous variables referring to the operating
variables. These operating variables include powers consumed, power flows,
power produced by generating units, load shed by demands, and voltage angles at buses (see reference [15] for a detail description of the formulation).
Π represents the maximum budget for investment in transmission expansion.
Finally, Ω(x, d) defines the feasibility region for the operating variables y, as
a function of investment decisions x and given realizations of the uncertain
parameters d, as follows:

Ax + By = E : λ



F x + Gy ≤ K : µ
Ω(x, d) =
(4)
I eq y = d : α



I ineq y ≤ d : ϕ,

where A, B, E, F , G and K are matrices of constant parameters dependent on the network configuration and element characteristics, I eq selects
5

the components of y that are equal to the uncertain parameters (demands),
and I ineq selects the components of y that are limited by the uncertain parameters (i.e. maximum power generation and maximum load shedding).
The first set of equality constraints correspond to constraints enforcing the
power balance at every bus, the power flow through each line, and fixing the
voltage angle of the reference bus. The first set of inequality constraints are
associated with line flow limits, and limits on the voltage angles at every bus.
Note that λ, µ, α and ϕ are the dual variable vectors associated with each
set of constraints, respectively.
For a detailed physical interpretation of the mathematical formulation
(1)-(3) we recommend the paper by [14].
2.1. Bi-level approach proposed by [15]
[15], analogously to [12], proposed dealing with problem (1)-(3) by decomposing and iteratively solving a subproblem and a master problem. The master variables correspond to x, i.e. the vector of first stage binary variables.
Thus, for given values of these master variables, the subproblem corresponds
to:
Maximum
Minimum bT y .
(5)
d∈D
y ∈ Ω(x, d)
Considering the dual problem associated with the third-level, (5) results
in the following single level maximization problem:
Maximize f dual = (E − Ax)T λ − (K − F x)T µ + dT (α − ϕ)
d, λ, µ, α, ϕ

(6)

subject to
B T λ − GT µ + I Teq α − I Tineq ϕ = b
µ ≥ 0
ϕ ≥ 0,

(7)
(8)
(9)

d ∈ D.

(10)

Regarding the uncertainty set D, [15] proposed the use of the following

6

set:

X

G
ˆG G
dG
= d¯G
i
i − di zi ; ∀i ∈ Ω
D
ˆD D
dD
= d¯D
i
i + di zi ; ∀i ∈ Ω

(11)
(12)

ziG ≤ ΓG
r ; ∀r

(13)

ziD ≤ ΓD
r ; ∀r

(14)

ziG ∈ {0, 1}; ∀i ∈ ΩG
ziD ∈ {0, 1}; ∀i ∈ ΩD ,

(15)
(16)

i∈(ΩG ∩ Ωr )

X

i∈(ΩD ∩

Ωr )

where dG
i is the uncertain generation limit related to the ith variable within
ˆG
vector d, d¯G
i is the corresponding nominal value, di is the maximum positive
distance from the nominal value that can take the random parameter, and
¯D ˆD
ΩG is the set of indices of the generating units. Analogously, dD
i , di , di
and ΩD correspond to the same values but for demands. The only limitation
D
introduced by this simplified set is that uncertainty budgets ΓG
r and Γr
must be integer values, which in our opinion does not dwarf the benefits of
robust optimization from the practical perspective. Alternatively, instead of
uncertainty budgets associated with generation, demand and regions ΓG
r and
D
Γr , a unique uncertainty for each region Γr , or for the system Γ could be
used instead.
Note that this set takes advantage of the simplification proposed by [17]
because the uncertain parameters might only take either the nominal, or
upper or lower limits of their uncertainty range, and the fact noticed by [14]
that the worst realization of “nature” try to make generation capacities to
be as lower as possible and demand loads as higher as possible. This set
contains a reduced number of binary variables with respect to [12], which
reduce the number of binary variables by half, increasing the percentage of
reduction with respect to [13] and [14].
Subproblem (6)-(10) is a bilinear mathematical programming problem,
which can be linearized and transformed into a mixed-integer linear mathematical programming problem at the expense of introducing binary variables
associated with the uncertainty set (for more details, see the transformation
of dT (α − ϕ) introduced by [15]).
The optimal solution of subproblem (6)-(10) provides the uncertain parameter values d(i) to construct primal cuts for the master problem, which
7

corresponds to the following optimization problem at iteration k:
Minimize
cT x + γ
x, y (i) ; ∀i = 1, . . . , k − 1

(17)

γ ≥ bT y (i) ; ∀i = 1, . . . , k − 1
γ ≥ 0
Ax + By (i) = E; ∀i = 1, . . . , k − 1
F x + Gy (i) ≤ K; ∀i = 1, . . . , k − 1
I eq y (i) = d(i) ; ∀i = 1, . . . , k − 1
I ineq y (i) ≤ d(i) ; ∀i = 1, . . . , k − 1,

(18)
(19)
(20)
(21)
(22)
(23)

subject to

cT x ≤ Π
x ∈ {0, 1}.

(24)
(25)

Note that the master problem, besides variable γ, includes one operating
variable vector y (i) for each realization of the uncertain parameters obtained
from subproblem (6)-(10).
3. Proposed solution technique
This section justifies and explain in detail the proposed solution technique. However, before focusing on the mathematical description of the solution method, we justify the type of uncertainty set used in this paper.
3.1. Ellipsoidal uncertainty set
Analogously to [3] and [6], we propose the following uncertainty set. Let
us assume that uncertain parameter vector d = (dG , dD )T has nominal or
expected values d̄ and variance-covariance matrix Σ, respectively. According
to [3], the ellipsoidal uncertainty set can be written using the Mahalanobis
distance as follows:
n
o
T

D(β) = d − d̄ (Σ)−1 d − d̄ ≤ β 2 ,
(26)

which constitutes the uncertainty set in (10). In contrast with the method
proposed by [15], one unique parameter β controls the size and protection
level of the ellipsoidal set.
8

Note that only first and second order moments of the random parameters
without lower and upper bounds are initially considered in this paper for the
ellipsoidal uncertainty set. There are two reasons, firstly, the simplification
of the process as it will be shown in next subsections, and secondly,
because


it is straightforward to include also the interval limitations d̄ − d̂, d̄ + d̂
of the random parameters within the proposed solution method.
In order to simplify the problem, the uncertainty set (26) is transformed
using an affine mapping into a ball of radius β, resulting in the following set
of alternative constraints:
d = d̄ + Lz; kzk ≤ β,

(27)

where L is the mapping matrix which can be obtained from Cholesky decomposition of variance-covariance matrix Σ = LLT , z represents a perturbation
vector and k · k stands for Euclidean norm. Note that using also the fact noticed by [14] that the worst realization of “nature” try to make generation
capacities to be as lower as possible and demand loads as higher as possible,
we could reduce the uncertainty set by half including the following additional
constraints:
G

(28)

D

(29)

dG ≤ d̄

dD ≥ d̄ ,

which force maximum power productions and demand loads to take, respectively, values below and above their corresponding nominal values.
3.2. Relationship of subproblem with First Order Second Moment methods
Structural reliability techniques have been applied successfully to solve
certain class of stochastic programming problems [18], and also robust optimization design problems [6]. Trying to establish an analogy between subproblem (5) and the problem treated in [18], which attempts to optimize a
given quantile of the random objective function, we could express the former
as follows:
f S (d) = Minimum (bT y),
(30)
y ∈ Ω(x, d)
however, no information about the worst realization of nature is included
yet, for that reason it is necessary to use the optimal value dML of the uncertain parameters within the uncertainty set (27)-(29) which provide the
9

worst-minimum operational costs. It is obtained from solving the following
problem:
qβ = maximum [f S (d)] .
(31)
d ∈ D(β)
Let assume that dML is known, in that case it is straightforward to calculate the optimal solution of subproblem (5) and values of the operating
variables y ∗ as follows:
qβ =

Minimum bT y,
y ∈ Ω(x, d)

(32)

d = dML : η

(33)

subject to
where η is the dual variable vector associated with constraint (32), which provides information about how much the optimal objective function changes
when uncertain parameters change around dML . For those values of the optimal operating variables y ∗ , we could try to find a new value of the uncertain
parameter vector in the vicinity of dML but providing an even worst operating
cost as follows:
∂qβ T
(d − dML) = qβ + η T (d − dML )
Maximize qβ +
∂d
d

(34)

subject to
dG ≤ d̄
D

G

(35)

D

(36)
(37)

d ≥ d̄ ,
d = d̄ + Lz,
ztz

(1/2)

≤ β.

(38)

Since the ellipsoidal uncertainty set defines a convex region, the optimal
solution dML must be located at the boundary of the ellipsoidal uncertainty
set, i.e. constraint (38) must hold as an equality constraint. The optimal
solution of the random parameters from problem (34)-(38) corresponds to
dML, otherwise it would mean that it is possible to obtain another value of
the random parameters inside the ellipsoidal set providing a worst value of the
operational cost, so that dML would not correspond to the optimal solution,
contradicting our initial hypothesis. This result might seem rather obvious,
10

however, this formulation (32)-(38) allows the proposal of a decomposition
technique to look for the optimal solution of subproblem (5).
Before describing the algorithm, let transform the formulation (34)-(38)
into the following equivalent problem at the optimum (see reference [18] for
the formal proof):

(1/2)
β = Minimum
(z t z)
,




d


subject to
(39)
d = d̄ + Lz,





∂qβ T
ML
qβ +
(d − d ) ≥ qβ . 
∂d
Note that the last constraint in (39) is the linear approximation of the
equation:
g(x, y ∗ , d) = f S (d) ≥ qβ
(40)

which avoids using operating variables and solving (30) to get (40). This
problem structure is analogous to that given by [18], however in that work
the objective function depends explicitly on the uncertain parameters and
there is no need to use a first order approximation.
Problem (39) corresponds to the formulation associated with the structural reliability method FOSM [see 19, 20, 21, 22, 23, 24, 25] used to calculate the probability of obtaining an operational cost higher than qβ , i.e.
Prob(f S (d) > qβ ). The random variables involved d belong to an n-dimensional
space, which for given values of the decision x and operating y ∗ variables
can be divided into two regions with respect to the limit-state equation
g(x, y ∗ , d) = f S (d) − qβ , the undesirable (U) and desirable (D) regions. That
is,
U ≡ {d|g(x, y ∗ , d) = f S (d) − qβ > 0},
(41)
D ≡ {d||g(x, y ∗ , d) = f S (d) − qβ ≤ 0}.
Note that once the qβ associated with the worst operational costs within
the uncertainty set is known, obtaining a greater value of that cost associated
with another realization of the uncertain parameters outside the uncertainty
set is an undesirable situation.
In this context, the dependent normally distributed random variables d
are transformed into independent standard normal random variables z that
can be integrated straightforwardly using the linear Hasofer and Lind transformation [20], L is the Cholesky factorization of the variance-covariance
11

matrix, i.e. Σ = LLT . This factorization is always possible for positive definite matrices, and variance-covariance matrices must be positive definite by
definition because their eigenvalues represent variances which must be always
positive. This transformation is equivalent to the mapping transformation in
(27). The optimal solution z ∗ corresponds to the closest point to the origin
located on the limit-state equation in the standard and independent normal
random space, β is the minimum distance so-called reliability index in the
structural reliability scientific community, and d∗ = dML is the point of maximum likelihood, i.e. the actual values of the uncertain parameters located on
the limit-state equation where the probability is higher, and it represents the
most likely values of the random parameters that produce a non desirable
operational cost.
The final probability of obtaining an undesirable operational cost is related to the reliability index by the relation (see [6] or [18] for more details):

Prob f S (d) > qβ = Φ(−β),
(42)
and the probability of a desirable operational cost is

Prob f S (d) ≤ qβ = 1 − Φ(−β) = Φ(β),

(43)

where Φ(·) is the cumulative distribution function of the standard normal
random variables. This method provides the exact probability if the limitstate equation is linear in the standard normal random space, i.e., if the
resulting limit-state distribution is normally distributed, which is the case if
uncertain parameters are normally distributed. Nevertheless, it is an approximation in case parameters do not follow a multivariate normal distribution,
which is usually the case becasue we only use expectations d̄ and variancecovariance matrix Σ. According to expressions (42) and (43), qβ corresponds
to a upper quantile of the uncertain minimum operational cost random variable, or an approximation of the quantile, depending on whether the random
parameters involved are normally distributed or not.
From this perspective, subproblem (5) can be interpreted as finding the
quantile of the uncertain minimum operating costs given by expression (43).
According to [18] the solution method proposed in this paper attempts to
minimize the Value-at-Risk of total transmission costs.
Figure 1 shows a graphical illustration of the interpretation of subproblem
5 under FOSM method. Panel (a) shows the circumferences associated with
12

different transformed joint probability density function contours for the bidimensional case (two random variables involved), including that associated
with the ellipsoidal uncertainty set (gray dashed line) whose distance from the
origin is equal to β. It is also shown the operating cost contour equal to the
quantile qβ , which is tangent with respect the ellipsoidal set circumference at
the point of maximum likelihood. Note that combination of random variables
whose operating costs are higher than qβ are in the undesirable region. If
random variables are normally distributed, those operating cost contours
are linear in the transformed space and the multidimensional problem can
be reduced to one single dimension related to the standard normal random
variable (see panel (b)). The probability of the operating costs to be inside
the undesirable region, which requires a multidimensional integration can be
calculated as Φ(−β) and it corresponds to the gray shadow area in panel (b).
This probability can also be represented in the space of the operating cost
random variable (see panel (b)), where Φ(−β) is equal to the shadow area
associated with operating costs exceeding qβ .

Design or Maximum
likelihood point

Uncertainty set
fZ (z1 ,z2)

z2

Φ(− β )

Joint probability density
function contours
β

Increasing operating
cost contours

N(0,12)

0
(b)

β

Undesirable
region

z1

Φ(− β )

Desirable region

qβ contour
qβ
(a)

Oper. cost
(c)

Figure 1: 3-D graphical interpretation of the FOSM method.

13

3.3. Proposed decomposition method
Once the uncertainty set is properly defined and the relationship between
subproblem (5) and FOSM methods has been established, we focus on our
methodological proposal for solving problem (1)-(3). It is based on iteratively
solving the following problems:
Subproblem: Consisting on formulation (32)-(38), which is the equivalent
with respect to (5) including the uncertainty set definition (27)-(29).
Master problem: It consist of the minimization of the objective function
(17) subject to constraints (18)-(25).
This iterative scheme is known as “outer” because it tries to solve the
complete TNEP problem. This outer iterative scheme have proved to be
successful for solving TNEP problems using cardinality constrained uncertainty sets, however, the key for the success of the proposed methodology is in
the resolution of the min-max subproblem. We propose a block coordinate
descent algorithm [26] that takes advantage of sub-subproblems (32)-(33)
and (34)-(38) to solve each of them optimizing the objective function over
a subset of variables while the remaining variables are fixed to given values. This process is iteratively solved until no further improvement of the
objective function is achieved. This iterative scheme is known as “inner”.
Convergence of this algorithm is only guaranteed for convex problems where
each optimization sub-subproblem has a unique optimizer (strict convexity
satisfies this requirement), which is the case for both sub-subproblems in this
specific application.
Before proceeding to explain how the subproblem can be solved, the proposed outer iterative scheme is described step by step in the following algorithm, which is analogous to that proposed by [14] and in [15]:
Algorithm 3.1. (Outer iterative method).
Input: Selection of uncertainty factor β and the tolerance of the process ε.
These data are selected by the decision maker.
Step 1: Initialization. Initialize the iteration counter to ν = 1, and upper
and lower bounds of the objective function z (up) = ∞ and z (lo) = −∞.

14

Step 2: Solving the master problem at iteration ν. Solve the master
problem (17) subject to constraints (18)-(25). The result provides values of the decision variables x(ν) and γ(ν) . Update the optimal objective function lower bound z (lo) = cT x(ν) + γ(ν) . Note that at the first
iteration the optimal solution corresponds to the no investment case,
alternatively, we could start with any other vector of decision variables.
Step 3: Solving subproblem at iteration ν. For given values of the decision variables x(ν) , we calculate the worst operating costs within the
(ν)
uncertainty set or quantile qβ , obtaining also the corresponding uncertain parameters dML
(ν) . This is achieved by solving sub-subproblems
(32)-(33) and (34)-(38). Update the optimal objective function upper
(ν)
bound z (up) = cT x(ν) + qβ .
Step 4: Convergence checking. If (z (up) − z (lo) )/z (up) ≤ ε go to Step 5,
else update the iteration counter ν → ν + 1 and continue in Step 2.
Step 5: Output. The solution for a given tolerance corresponds to x∗ =
x(ν) .

Next step is the description of the method to solve the subproblem, as
pointed out previously it consists on a block coordinate descent algorithm
that takes advantage of formulations (32)-(33) and (34)-(38). The key is to
chose an initial value of the uncertain parameter vector at iteration κ = 1, as
candidate for maximum likelihood point dML
(κ) . We recommend the selection of
low values of generation capacities and high values of demand loads because
we know that the optimal maximum likelihood point is around this area of
the uncertainty set. Once this selection is done, sub-subproblem (32)-(33) is
solved replacing constraint (33) by:
d = dML
(κ) : η (κ) .
(κ)

(44)

The solution provides qβ , the optimal values of operating variables y ∗(κ)
for that particular realization of the random parameters and the derivatives
of the corresponding worst operational cost with respect to the uncertain
parameters, i.e. η (κ) . Next, we try update the point of maximum likelihood
in a attempt to find either a point inside the ellipsoidal uncertainty set, this
15

case occurs if the initial point dML
(κ) is outside the ellipsoidal set, or a point in
the boundary of the ellipsoidal set with an even worst operational cost. This
is done solving sub-subproblem (34)-(38) replacing the objective function
(34) by:
(κ)
Maximize qβ + η T(κ) (d − dML
(45)
(κ) ) .
d
The advantage of using this sub-subproblem is that [6] provides the analytical solution, which for this particular case is as follows:
dML
(κ+1) = d̄ + β q

Ση (κ)
η T(κ) Ση (κ)

.

(46)

This inner” iterative process continues until no further improvement in
the worst operational cost is achieved, and/or the maximum likelihood point
does not change between consecutive iterations. The corresponding algorithm is summarized as follows:
Algorithm 3.2. (Inner iterative method).
Input: The same input arguments than the outer algorithm are used and
also the actual values of the decision variables from the outer problem
x(ν) , which are constants during this iterative process.
Step 1: Initialization. Initialize the iteration counter to κ = 1, and the
initial values of the random parameters dML
(κ) .
Step 2: Determining the minimum operational cost at iteration κ.
For given values of the random parameters dML
(κ) , obtain the operating
variables providing the lowest operational costs for the actual iteration
by solving problem (32)-(33) but replacing constraint (33) by (44). The
(κ)
result provides qβ , the optimal values of operating variables y ∗(κ) for
that particular realization of the random parameters and the derivatives of the corresponding worst operational cost with respect to the
uncertain parameters, i.e. η (κ) .
Step 3: Convergence checking. If the iteration counter is equal to (κ =
1) go to Step 4, otherwise, proceed to check convergence, i.e., if ||dML
(κ) −
ML
d(κ−1) || ≤ ε go to Step 5, else continue in Step 4.
16

Step 4: Evaluating the new points of maximum likelihood. For given
values of the operating variables y ∗(κ) , we update the point of maximum
likelihood using the analytical expression (46). Thus, updated representative values of the random parameters for the next iteration dML
(κ+1)
are obtained. The iteration counter must then be updated κ → κ + 1
before proceeding to Step 2.
(κ)

Step 5: Output. The solution for a given tolerance corresponds to qβ ,
y ∗(κ) and dML
(κ) .

3.4. Bounded random variables
So far, the method developed in this paper deals with uncertainty coefficients by using the first and second probability moments associated with
their probability distributions, and without assuming any particular type of
distribution. This feature makes the method to be an approximation unless
the random parameters truly follows a normal distribution. However, most
of robust uncertainty sets proposed in the literature deal with random variables within a given interval. This information could be easily incorporated
in the proposed method by including the following constraint set into the
subsubproblem (34)-(38):
−d̂ ≤ d − d̄ ≤ d̂.

(47)

These constraints constitute bounded convex polytopes both in the original space associated with random model parameters d and also in the transformed space z. To solve this modified subproblem, the best strategy is as
follows:
1. To use the analytical solution (46) and check afterwards whether constraints (47) hold, which means that the optimal solution has been
achieved.
2. Otherwise, it means that the solution is not optimal according to those
bounds, then the optimal solution is found by solving problem (34)(38), including constraints (47) and removing the quadratic constraint
(38). This results in a linear mathematical programming formulation.
3. Finally, in case the solution from item 2 is also out of the ellipsoidal
uncertainty set, the complete problem (34)-(38) including constraints
17

(47) must be solved. This is a linear mathematical programming problem with one quadratic constraint that can be solve efficiently using
interior point algorithms.
Alternatively, after checking the analytical solution from item 1, the conic
problem in item 3 could be used directly instead of the linear problem in item
2.
4. Numerical case studies
In this section, we present numerical experiments of our model proposal
and its comparison with the method proposed by [15], which proved to be
the most computationally efficient method to solve the TNEP problem using
cardinality constraint uncertainty sets. We use the same examples provided
by [15], i.e. the Garver system illustrative example [27], and two realistic
case studies: the IEEE 24-bus system [28] and the IEEE 118-bus test system
[29]. To avoid providing redundant information numerical data will not be
displayed in this manuscript, those readers interested in that information can
find it in [15] and references therein.
All examples have been implemented and solved using GAMS [30, 31]
and CPLEX 12, respectively, on a PC with four processors clocking at 2.39
GHz and 3.2 GB of RAM memory. Note that computing times reported
correspond to the sum of running times used by the solvers in order to solve
the masters and subproblems until the final solution is attained. We use the
GAMS command model.resusd. We use the same tolerance for all problems
and equal to ǫ = 10−6 .
4.1. Illustrative example. Garver system
The proposed model is illustrated with the Garver 6-bus system, depicted
in Figure 2. This system comprises 6 buses, 3 generators, 5 inelastic demands
and 6 lines. Nominal values of generation capacities and demands and their
offering and bidding prices can be found in [15]. The load-shedding cost is
equal to the bid price of each demand. It is considered that a maximum of
three lines can be installed between each pair of buses. Line data are obtained
from Table I of [32], including construction costs, and the maximum available
investment budget is e40 million.
The investment return period of each line is considered to be 25 years, and
the discount rate is 10%, resulting in an annual amortization rate of 10%.
18

Existing lines
New lines null correlation
New lines negative correlation
240MW
80MW
bus 5

150MW
bus 1

350MW
bus 3
40MW

bus 2

240MW

600MW
bus 6

bus 4
160MW

Figure 2: Garver’s 6-bus test system.

Since investment cost is annualized, the weighted factor σ is equal to the
number of hours of a year, i.e. 8760, obtaining an annualized load-shedding
and power generation cost.
Regarding the uncertainty set, the standard deviations of power generation capacity and demand load are equal to 50% and 20% of their nominal values, respectively, divided by 2.3263. Note that we have selected this
standard deviation so that in case the random parameters are normally distributed, the probability of those random parameters to be within the range
defined by [15] is equal to 0.99. Initially the random parameters are assumed
to be independent.
We have solved the robust TNEP problem for different reliability indexes
β. These values have been selected in an attempt to find similar objective
functions with respect to those obtained in [15] for different combinations of
the uncertainty budgets ΓG and ΓD . In order to obtain statistically sound
conclusions about computing times, we have solved the problem 100 times
and obtained mean and standard deviation executions times. Results from
this numerical experiments are given in Table 1.

19

Table 1: Computational results for Garver’s 6-bus example
[15] method
Proposed method
Optimal
Optimal
♯ iter. Mean (s) Std. (s) β
♯ iter. µt (s) σt (s)
sol. (Me)
sol. (Me)
ΓG
ΓD
ΓG
ΓD
ΓG
ΓD
ΓG
ΓD

=3
=5
=2
=3
=1
=2
=0
=0

35505.31

3

0.644

0.054 6.58 35505.31

2

0.16

0.02

25832.22

4

0.936

0.094

4.3

25841.82

3

0.95

0.05

5861.92

4

1.011

0.072

2.6

5315.16

3

0.44

0.06

440.07

2

0.428

0.065

0

440.07

2

0.17

0.06

From this table the following observations are pertinent:
1. The proposed method is comparable in average computational time
with respect to the method proposed by [15] for all cases, from the first
and worst possible case corresponding to Soyster’s solution [1] up to the
deterministic case when uncertainty budgets are null. Regarding the
standard deviation of computing times the values are also comparable.
2. Both methods require an small number of outer iterations.
3. The worst possible solution associated with β = 6.58 is the same than
the analogous proposed by [15]. The reason is that all random parameters take their minimum and maximum possible
qP values associated with
8
2
zi -values, i.e. equal to 2.3263. Note that
i=1 2.3263 = 6.58.
4. The deterministic solution is also the same because it corresponds to
the case where all uncertain parameters are equal to their nominal
values.
5. For intermediate situations similar solutions have been found, which
according to the reliability indexes β = 4.3 and β = 2.6 used and assuming normally distributed parameters, they correspond to quantiles
Φ(β) of 0.99999146 and 0.99533881.
If we focus on the robust solution associated with the reliability index
β = 4.3, the lines constructed are shown in Figure 2 in broken lines. At
the point of maximum likelihood, only generation capacities of generators
20

at nodes 3 and 6, i.e those with higher capacities, reach their minimum
values. The expansion investment is equal to 21.239Me. However, let assume
that those power generators correspond to wind farm fields and their power
productions are negatively correlated with coefficient −0.8, this makes the
possibility of both generators at nodes 3 and 6 to be at their minimum
capacities simultaneously remote. If we solve the robust problem under this
new circumstance, the annualized optimal solution is reduced to 7739.721Me
with the same level of uncertainty, however, the cost of investment is higher
38.616Me because besides the additional lines constructed in the previous
case, three additional lines joining buses 4 and 6 are constructed (shown in
gray continuous lines in Figure 2). In this case, at the point of maximum
likelihood, only generation capacities of generators at nodes 1 and 3 reach
their minimum values. This result proves that the correlation structure of
generation capacities associated with renewable and intermittent sources is
crucial for the appropriate expansion of the network.
Finally, we perform an additional simulation experiment using a reliability index β = Φ−1 (0.9) = 1.28155 to check whether the optimal solution of
the subproblem corresponds to the 0.9-quantile. Note that we have selected
this level of protection because the maximum likelihood point at the optimal
solution does not contain any limit of the random variables, in addition we
divide the standard deviations of all random variables by 10 to make sure
that we avoid sampling negative values of generations capacities and demand
loads. This is important if we are interested in checking whether or not the
optimal qβ∗ = 440.227Me for this example corresponds to the 0.9-quantile.
We also consider the negative correlation between generation capacities at
nodes 3 and 6. The procedure is simple, we sample random values of generation capacities and demand loads using the following expression:
dsim = d̄ + Lz sim

(48)

where z sim are randomly generated samples from an independent standard
normal distribution. For each realization of the uncertain parameters, problem (32)-(33) replacing (33) by the expression d = dsim is solved, which
allows calculating the minimum operating costs qβsim for those uncertain parameters. Repeating the process 1000 times, the probability of not exceeding
the optimal value qβ∗ = 440.227 is equal to the number of times the cost is
lower or equal than qβ∗ divided by the total number of samples. The probability obtained from Monte Carlo simulation is 0.907 which confirms that
the optimal cost corresponds to the 0.9-quantile.
21

Figure 3 shows the histogram associated with operational cost simulation results, the normal distribution fitted and also the optimized quantile
qβ∗ = 440.227. Note that the cost distribution is normal, because uncertain
parameters were supposed to be normally distributed. For this reason the
quantile minimization is exact.
Garver simulation experiment
0.2
0.18

Sample data
Normal fit normal
qβ

0.16
0.14

Φ( β )
Density

0.12
0.1
0.08

Φ(− β )
0.06
0.04
0.02
0
430

432

434

436
438
Operational cost (Meuro)

440

442

444

Figure 3: Simulation experiment for Garver’s 6-bus test system.

Although this result is satisfactory, let us remind readers that this exact probabilities have been obtained under very specific circumstances which
barely happen in practice, for the rest of cases the solution is an approximation of the quantile.
4.2. IEEE 24-bus Reliability Test System
The following case study described is based on the IEEE 24-bus Reliability
Test System (RTS) [33], depicted in Fig. 4. The system comprises 24 buses,
34 existing corridors which can accept a maximum of three equal lines and
7 new corridors, 10 generating units and 17 loads. Data for lines in existing
corridors are taken from [33], while line data for new corridors are obtained
from Table I in [34]. Investment costs are e20 million, and analogously to
22

the Garver illustrative example, the investment return period of each line
is considered to be 25 years, and the discount rate is 10%, resulting in an
annual amortization rate of 10%. Data related to location of generators
230 kV
22

18
21
17

23

16

19
14

Synch.
Cond.

15

24

20

11

3

13

12

10

9

6
4
5

8

138 kV

1

2

7

Existing lines
New lines null correlation
New lines negative correlation

Figure 4: IEEE 24-bus reliability test system (RTS)

and demands in the network, maximum power offered, prices for generating
units, power bid and the bid prices for demands can be found in [15]. The
load-shedding cost is equal to 100 times the bid price of each demand.
Regarding the uncertainty set, power generation capacity can increase
or decrease a maximum of 50% of their nominal values, while demand are
allowed to change a maximum of 20%. Regarding the uncertainty set and
analogously to the previous case, the standard deviations of power generation
capacity and demand load are equal to 50% and 20% of their nominal values,
respectively, divided by 2.3263. Note that we have selected this standard
23

deviation so that in case the random parameters are normally distributed,
the probability of those random parameters to be within the range defined
by [15] is equal to 0.99. Initially the random parameters are assumed to be
independent.
We have solved the robust TNEP problem for different reliability indexes
β. These values have been selected in an attempt to find similar objective
functions with respect to those obtained in [15] for different combinations of
the uncertainty budgets ΓG and ΓD . In order to obtain statistically sound
conclusions about computing times, we have solved the problem 100 times
and obtained mean and standard deviation executions times. Results from
this numerical experiments are given in Table 2.
Table 2: Computational results for IEEE 24 RTS case study without correlation

[15] method
Optimal
♯ iter. Mean (s) Std. (s)
sol. (Me)
ΓG = 10
ΓD = 17
ΓG = 7
ΓD = 12
ΓG = 3
ΓD = 5
ΓG = 0
ΓD = 0

β

500752.14

3

1.51

0.30

454917.47

3

1.80

0.39

9

348941.40

3

1.12

0.11

219161.52

2

0.414

0.08

Proposed method
Optimal
♯ iter. µt (s) σt (s)
sol. (Me)

12.09 500752.14

2

1.34

0.19

455488.75

3

2.35

0.71

4.5

345282.84

3

0.78

0.10

0

219161.52

2

0.28

0.06

From Table 2 we can extract the same conclusions than in the previous
illustrative example:
1. The proposed method is comparable in average computational time
with respect to the method proposed by [15] for all cases, from the first
and worst possible case corresponding to Soyster’s solution [1] up to the
deterministic case when uncertainty budgets are null. Regarding the
standard deviation of computing times the values are also comparable.
2. Both methods require an small number if outer iterations.
3. The worst possible solution associated with β = 12.09 is the same than
the analogous proposed by [15]. The reason is that all random parame24

ters take their minimum and maximum possible
qP values associated with
27
2
zi -values, i.e. equal to 2.3263. Note that
i=1 2.3263 = 12.09.
4. The deterministic solution is also the same because it corresponds to
the case where all uncertain parameters are equal to their nominal
values.
5. For intermediate situations similar solutions have been found. For
instance, the optimal solution associated with the reliability index
β = 4.5, and assuming normally distributed parameters, corresponds
to quantile Φ(β) of 0.999996602326875.
If we focus on the robust solution associated with the reliability index
β = 9, the lines constructed are shown in Figure 4 in broken lines. Only three
new lines are selected, from bus 3 to 9, from 11 to 14 and from 14 to 16. The
investment cost is 14.307Me. At the point of maximum likelihood, there
are several generation capacities reaching their minimum possible values, for
instance generators 3 and 4 at buses 7 and 13, respectively. If we assume
that those power generators correspond to wind farm fields and their power
productions are negatively correlated with coefficient −0.8, this makes the
possibility of both generators at nodes 3 and 4 to be at their minimum capacities simultaneously remote. If we solve the robust problem under this new
circumstance, the annualized optimal solution is reduced to 432351.604Me
with the same level of uncertainty, however, the cost of investment is higher
16.606Me and instead of constructing the lines from the previous case, only
two lines from bus 3 to 24 and from 15 to 24 are selected as new lines
(shown in gray continuous lines in Figure 4). This result also confirms that
the correlation structure of generation capacities associated with renewable
and intermittent sources is important for the appropriate expansion of the
network.
Finally, we perform an additional simulation experiment using a reliability
index β = Φ−1 (0.9) = 1.28155 to check whether the optimal solution of the
subproblem corresponds to the 0.9-quantile. Note that we have selected this
level of protection because the maximum likelihood point at the optimal
solution does not contain any limit of the random variables, which is the
case with the standard deviation used in the example.
This is important if we are interested in checking whether or not the optimal qβ∗ = 251719.780Me for this example corresponds to the 0.9-quantile.
We also consider the negative correlation between generation capacities at
nodes 7 and 13. The procedure is analogous to that used in the previous
25

illustrative example. Repeating the process 1000 times, the probability of
not exceeding the optimal value qβ∗ = 251719.780 is equal to the number of
times the cost is lower or equal than qβ∗ divided by the total number of samples. The probability obtained from Monte Carlo simulation is 0.889 which
confirms that the optimal cost corresponds to the 0.9-quantile.
Figure 5 shows the histogram associated with operational cost simulation
results, the normal distribution fitted and also the optimized quantile qβ∗ =
251719.780. Note that the cost distribution is normal, because uncertain
parameters were supposed to be normally distributed. For this reason the
quantile maximization is exact.
IEEE 24 RTS simulation experiment

5

1.5

x 10

Sample data
Normal fit normal
qβ

Φ( β )

Density

1

Φ(− β )

0.5

Φ(− β )

0

1

1.5

2
2.5
3
Operational cost (Meuro)

3.5

4
5

x 10

Figure 5: Simulation experiment for IEEE 24 RTS test system.

Although this result is satisfactory, let us remind readers that this exact probabilities have been obtained under very specific circumstances which
barely happen in practice, for the rest of cases the solution is an approximation of the quantile.

26

4.3. IEEE 118-bus test system
Finally, we run additional computational tests using the IEEE 118-bus
test system [29]. The system comprises 118 buses, 186 existing lines, 54
generating units and 91 loads. In addition, it is possible to construct up to
61 additional lines to duplicate each one of the following existing lines: 8,
12, 23, 32, 38, 41, 51, 68, 78, 96, 104, 118, 119, 121, 125, 129, 134, 159, 7,
9, 36, 117, 71, 131, 133, 147, 103, 65, 144, 168, 4, 13, 132, 69, 66, 67, 5,
89, 29, 167, 145, 70, 42, 90, 16, 174, 98, 99, 185, 93, 94, 128, 164, 97, 153,
146, 116, 163, 31, 92, 130. Data for lines in existing corridors are taken from
[29]. Investment costs are e100 million, and analogously to both previous
examples, the investment return period of each line is considered to be 25
years, and the discount rate is 10%, resulting in an annual amortization rate
of 10%.
Data for generation capacities and demand loads are given in [15]. The
load-shedding cost is equal to 1.2 times the bid price of each demand.
Regarding the uncertainty set, the standard deviations of power generation capacity and demand load are equal to 50% and 50% of their nominal values, respectively, divided by 2.3263. Note that we have selected this
standard deviation so that in case the random parameters are normally distributed, the probability of those random parameters to be within the range
defined by [15] is equal to 0.99. Random parameters are assumed to be
independent.
We have solved the robust TNEP problem for different reliability indexes
β. These values have been selected in an attempt to find similar objective
functions with respect to those obtained in [15] for different combinations of
the uncertainty budgets ΓG and ΓD . In order to obtain statistically sound
conclusions about computing times, we have solved the problem 100 times
and obtained mean and standard deviation executions times. Results from
this numerical experiments are given in Table 3.
From Table 3 the same observations than in previous examples are obtained. However, in this
qPcase, the worst possible solution is associated with
1
2
β = 28.02 because
i=1 452.3263 = 28.02. Note that computational
times are below 10 seconds, which prove the method to be applicable in real
world energy transmission networks.

27

Table 3: Computational results for IEEE 118-bus test system

[15] method
Optimal
♯ iter. Mean (s) Std. (s)
sol. (Me)
ΓG = 54
ΓD = 91
ΓG = 35
ΓD = 60
ΓG = 15
ΓD = 20
ΓG = 0
ΓD = 0

β

Proposed method
Optimal
♯ iter. µt (s) σt (s)
sol. (Me)

31994.36

4

7.73

0.80

28.02 31994.36

2

0.42

0.25

30032.93

5

53.46

60.36

21.0

29827.30

4

1.84

0.49

23352.97

4

24.26

26.05

9.5

21828.82

5

4.60

0.23

13929.23

2

0.56

0.07

0

13929.23

3

2.65

0.34

5. Conclusions
This paper extends the applicability of method [6], based on first-order
second-moment methods from structural reliability and decomposition techniques, into the adaptive robust optimization problem associated with transmission network expansion planning. The method is specially suitable for
problems where first and second order moments of the probability distributions of the uncertain parameters involved are available.
The proposed method has the following advantages:
1. The feasibility region defined by (27) is convex and only depends on
variables d.
2. The resulting sub-subproblem (32)-(33) for fixed values of the uncertain
parameters is linear and it does not content binary variables, in contrast with respect to TNEP problems based on cardinality constrained
uncertainty sets.
3. The sub-subproblem (34)-(38) resulting from fixing the operating variables to given values has analytical solution. In case additional bounds
on uncertain parameters are included, the worst possible scenario requires solving a linear programming problem with one quadratic constraint, which can be solved efficiently using interior point algorithms.
4. The outer iterative method between subproblem and master problem
provides upper and lower bounds of the objective function, which al28

lows checking if the problem is being solved appropriately during the
iterative process. Besides, it converges in a small number of iterations.
5. This method allows including the correlation structure of the random
variables involved easily, this is very important specially for wind power
farms at different sites and demand loads, which might have significant
correlations affecting the optimal expansion strategy.
6. The interpretation based on structural reliability allows selecting the
level of conservationism in terms of approximate probabilities of exceeding the design operational costs, because this parameter allows
selecting the quantile of the cost function to be optimized. This strategy is equivalent to the value-at-risk strategy widely used in financial
and stochastic applications.
In summary, this paper proposes a new decomposition algorithm to attain the exact solution of the TNEP problem derived from using a two-stage
adaptive robust strategy with ellipsoidal uncertainty sets. In addition, this
procedure set the basis for alternative and more sophisticated robust uncertainty sets based on reliability methods. This would allow incorporating not
only expected values and the correlation structure of uncertain parameters
but appropriate marginal distributions of those parameters. This is a subject
for further research.
Acknowledgements
This work has been partially funded by projects by Junta de Comunidades
de Castilla-La Mancha, under project POII-2014-012-P; the Ministry of Science of Spain, under CICYT Project ENE2012-30679; and by the European
Commission, under Grant Agreement Number 309048.
References
[1] A. L. Soyster, Convex programming with set-inclusive constraints and
applications to inexact linear programming, Oper. Res. 21 (5) (1973)
1154–1157. doi:10.1287/opre.21.5.1154.
[2] A. Ben-Tal, A. Nemirovski, Robust convex optimization, Mathematics of
Operations Research 23 (4) (1998) 769–805. doi:10.1287/moor.23.4.769.

29

[3] A. Ben-Tal, A. Nemirovski, Robust solutions of uncertain linear programs, Operations Research Letters 25 (1) (1999) 1 – 13.
doi:http://dx.doi.org/10.1016/S0167-6377(99)00016-4.
[4] A. Ben-Tal, A. Nemirovski, Robust solutions of linear programming
problems contaminated with uncertain data, Mathematical Programming 88 (3) (2000) 411–424. doi:10.1007/PL00011380.
[5] D. Bertsimas, M. Sim, The price of robustness, Oper. Res. 52 (1) (2004)
35–53. doi:10.1287/opre.1030.0065.
[6] R. Mı́nguez, Robust solutions of uncertain mixed-integer linear programs using structural-reliability and decomposition techniques, CoRR
abs/1409.8593.
URL http://arxiv.org/abs/1409.8593
[7] G. B. Dantzig, Linear programming under uncertainty, Management Sci.
1 (1955) 197–206.
[8] A. Thiele, T. Terry, M. Epelman, Robust linear optimization with recourse, Ioe technical report tr09-01, University of Michigan (2010).
[9] D. Bertsimas, E. Litvinov, A. Sun, J. Zhao, T. Zheng, Adaptative robust optimization for the security constrained unit commitment problem, IEEE Trans. Power Syst. 28 (1) (2013) 52–63.
[10] P. Wu, H. Cheng, J. Xing, The interval minimum load cutting problem
in the process of transmission network expansion planning considering
uncertainty in demand, IEEE Trans. Power Syst. 23 (3) (2008) 1497–
1506.
[11] H. Yu, Y. Chung, K. P. Wong, Robust transmission network expansion
planning method with taguchi’s orthogonal array testing, IEEE Trans.
Power Syst. 26 (3) (2011) 1573–1580.
[12] R. A. Jabr, Robust transmission network expansion planning with uncertain renewable generation and loads, IEEE Trans. Power Syst. 28 (4)
(2013) 4558–4567.
[13] B. Cheng, J. Wang, L. Wang, Y. He, Z. Wang, Robust optimization
of transmission expansion planning: Minimax cost vs. minimax regret,
IEEE Trans. Power Syst. (2014) 1–9doi:10.1109/TPWRS.2014.2313841.
30

[14] C. Ruiz, A. J. Conejo, Robust transmission expansion planning, European Journal of Operational Research (2014) 1–
12doi:http://dx.doi.org/j.ejor.2014.10.030.
[15] R. Mı́nguez, R. Garcı́a-Bertrand, Robust transmission network expansion planning in energy systems: Improving computational performance,
CoRR abs/1501.05480.
URL http://arxiv.org/abs/1501.05480
[16] B. Zeng, L. Zhao, Solving two-stage robust optimization problems using a column-and-constraint generation method, Operations Research
Letters 41 (5) (2013) 457–461.
[17] R. A. Jabr, Robust self-scheduling under price uncertainty using conditional value-at-risk, IEEE Trans. Power Syst. 20 (4) (2005) 1852–1858.
[18] R. Mı́nguez, A. J. Conejo, R. Garcı́a-Bertrand, Reliability and decomposition techniques to solve certain class of stochastic programming
problems, Reliability Engineering & System Safety (96) (2011) 314–323.
doi:10.1016/j.rees.2010.09.011.
[19] A. M. Freudenthal, Safety and the probability of structural failure,
Transactions, ASCE 121 (1956) 1337–1397.
[20] A. M. Hasofer, N. C. Lind, Exact and invariant second moment code
format, J. Engrg. Mech. 100 (EM1) (1974) 111–121.
[21] R. Rackwitz, B. Fiessler, Structural reliability under combined load sequences, Comput. Struct. 9 (1978) 489–494.
[22] O. Ditlevsen, Principle of normal tail approximation, J. Engineering
Mechanics Div., ASCE 107 (6) (1981) 1191–1208.
[23] M. Hohenbichler, R. Rackwitz, Non-normal dependent vectors in structural safety, J. Engineering Mechanics Div., ASCE 107 (6) (1981) 1227–
1238.
[24] R. E. Melchers, Structural reliability analysis and prediction, 2nd Edition, John Wiley & Sons, New York, 1999.

31

[25] R. Mı́nguez, M. Castillo, Reliability-based optimization in engineering
using decomposition techniques and FORMS, Structural Safety 31 (3)
(2009) 214–223.
[26] A. J. Conejo, E. Castillo, R. Mı́nguez, R. Garcı́a-Bertrand, Decomposition techniques in mathematical programming. Engineering and science
applications, Springer-Verlag Berlin Heidelberg, New York, 2006.
[27] L. L. Garver, Transmission network estimation using linear programming, IEEE Trans. Power Apparat. Syst. 89 (7) (1970) 1688–1697.
[28] R. T. S. T. Force, The IEEE Reliability Test System-1996, IEEE Trans.
Power Syst. 14 (3) (1979) 1010–1020.
[29] IEEE 118-bus system, Available: http://motor.ece.iit.edu/ data/data
118 bus.pdf.
[30] http://www.gams.com.
[31] A. Brooke, D. Kendrick, A. Meeraus, R. Raman, GAMS: A user’s guide,
GAMS Development Corporation, Washington, 1998.
[32] L. P. Garcés, A. J. Conejo, R. Garcı́a-Bertrand, R. Romero, A bi-level
approach to transmission expansion planning within a market environment, IEEE Trans. Power Syst. 24 (3) (2009) 1513–1522.
[33] Reliability Test System Task Force, The IEEE Reliability Test System
1996, IEEE Transactions on Power Systems 14 (3) (1999) 1010–1020.
[34] R. Fang, D. J. Hill, A new strategy for transmission expansion in competitive electricity markets, IEEE Trans. Power Syst. 18 (1) (2003) 374–
380.

32

