Robust Transmission Network Expansion Planning in Energy Systems: Improving Computational Performance

# Robust Transmission Network Expansion Planning in Energy Systems: Improving Computational Performance

## Abstract

In recent advances in solving the problem of transmission network expansion planning, the use of robust optimization techniques has been put forward, as an alternative to stochastic mathematical programming methods, to make the problem tractable in realistic systems. Different sources of uncertainty have been considered, mainly related to the capacity and availability of generation facilities and demand, and making use of adaptive robust optimization models. The mathematical formulations for these models give rise to three-level mixed-integer optimization problems, which are solved using different strategies. Although it is true that these robust methods are more efficient than their stochastic counterparts, it is also correct that solution times for mixed-integer linear programming problems increase exponentially with respect to the size of the problem. Because of this, practitioners and system operators need to use computationally efficient methods when solving this type of problem. In this paper the issue of improving computational performance by taking different features from existing algorithms is addressed. In particular, we replace the lower-level problem with a dual one, and solve the resulting bi-level problem using a primal cutting plane algorithm within a decomposition scheme. By using this alternative and simple approach, the computing time for solving transmission expansion planning problems has been reduced drastically. Numerical results in an illustrative example, the IEEE-24 and IEEE 118-bus test systems demonstrate that the algorithm is superior in terms of computational performance with respect to existing methods.

1

## 1 Introduction

The aim of Transmission Network Expansion Planning (TNEP) is to resolve the issue of how to expand or enhance an existing electricity transmission network in order to adequately satisfy demand for a given horizon. The main difficulty of this problem is decision-taking with the great amount of uncertainty associated with i) demand, ii) generation with renewable energy sources, such as wind and solar power plants, and iii) equipment failure. Two main frameworks have been used to tackle these uncertainties, stochastic programming (BirgeL:97, ) and robust optimization (Soyster:73, ; El-GhaouiL:97, ; El-GhaouiOL:98, ; Ben-TalN:98, ; Ben-TalN:99, ; Ben-TalN:00, ; BertsimasS:04, ).

In stochastic programming problems, uncertain data is assumed to follow a given probability distribution and is usually dealt with by using scenarios or finite sampling from the joint probability density function. Examples of the successful application of stochastic programming in TNEP problems are given by (CarrionAA:07, ; LopezPQ:07, ; GarcesCGR:09, ; RohSW:09, ) among others. However, the number of scenarios needed to represent the actual stochastic processes can be very large, which may result in intractable problems. Conversely, robust optimization can avoid the intractability issue related to stochastic programming approaches. Note that tractability is not the only difficulty stemming from stochastic programming, but it is the most important aspect from the point of view given in this paper. Examples of the successful application of robust optimization in TNEP problems are given by (WuCX:08, ; YuCW:11, ). However, in this paper we focus on the seminal works proposed by (Jabr:13, ), (ChenWWHW:14, ) and (RuizC:14, ).

The three works on which our proposal is based formulate the problem using an Adaptive Robust Optimization (ARO) framework (ThieleTE:10, ; BertsimasLSZZ:13, ). Therefore, all of them use three-level formulations: i) the first level minimizes the cost of expansion ((ChenWWHW:14, ) also minimizes the maximum regret), the decision variables for this level are those associated with construction or expansion of lines, ii) the 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 variables to minimize operating costs for given values of first and second level variables. Note that in this paper we do not consider contingencies explicitly, although, they could be taken into account within robust formulation following the guidelines set out by (StreetOA:11, ).

Jabr (Jabr:13, ) merges the second and third levels into one single-level maximization problem (so called subproblem in this paper) by using the dual of the third level. It takes advantage of the fact that all uncertain parameters need to be equal to their upper or lower limits in the worst case scenario, which allows the use of binary variables to deal with these uncertain parameters. The limitation of this simplification is that the uncertainty budget, typically used in robust optimization, must be an integer. However, in our opinion, this does not detract from the benefits of robust optimization ,while it does simplify the resolution of the problem substantially. In contrast, in references (ChenWWHW:14, ; RuizC:14, ) the second and third levels are merged into one single level maximization problem by using the Karush-Kuhn-Tucker (KKT) conditions for the third level. In this particular case, the number of constraints, continuous and binary variables of the subproblem increase with respect to the other approach.

Regarding the resolution of the resulting bi-level problem, Jabr (Jabr:13, ) sets out a Benders (ConejoCMG:06, ) approach where the dual information from subproblems is used to make additional Benders cuts. The main drawback of this approach is the slow convergence typical of this type of decomposition algorithms, which forced the author to include additional linear constraints in order to improve convergence. On the other hand, Cheng et al. and Ruiz and Conejo (ChenWWHW:14, ; RuizC:14, ) applied a constraint-and-column (ZengZ:13, ) generation method solely based on primal cuts. This method is computationally advantageous with respect to Benders decomposition and converges in a small number of iterations (authors reported solution results with only 3-4 iterations required). Note that the methodologies proposed by (ChenWWHW:14, ) and RuizC:14 () are basically the same from a mathematical point of view, although, they analyse different situations and Ruiz and Conejo RuizC:14 () set out a different way to construct the uncertainty sets.

In this paper an alternative approach to solving the three-level mixed-integer optimization problem associated with TNEP is set out. The concept elaborated on in this paper is to combine the way the inner level problems (subproblems) are solved by (Jabr:13, ) and the constraint-and-column generation method used by (ChenWWHW:14, ) and RuizC:14 () for the first level problem (master problem). Note that although the master problem solved by means of primal cuts is very efficient, computing times associated with the approaches adopted by (ChenWWHW:14, ) and RuizC:14 () increase exponentially with respect to the size of the problem. This is due to the subproblem, which consumes most of the resolution time due to the use of the KKT optimality conditions for the third level in the second level problem. The resulting model simultaneously calculates the third level primal and dual variables, and as the primal variables are not required to construct the primal cuts of the master, this alternative consisting in use of the KKT conditions is inefficient.

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 approaches for finding solutions set out in references (Jabr:13, ), (ChenWWHW:14, ) and RuizC:14 (). The detailed formulation is only given in the A. In Section 3 the definition of the uncertainty set and the description of the approach for finding a solution are set out. Numerical results for different networks are given in Section 4 and compared with those obtained using the method proposed by (ChenWWHW:14, ) and RuizC:14 (). Finally, the paper is concluded in Section 5.

## 2 Transmission network expansion planning problem: ARO compact formulation

According to references (Jabr:13, ), (ChenWWHW:14, ) and RuizC:14 (), the robust TNEP problem can be written in compact matrix form as the following three-level mathematical programming problem:

 Minimize\boldmathx\unboldmath⎛⎜⎝% \boldmathc\unboldmathT\boldmathx\unboldmath+cMaximum\boldmathd\unboldmath∈\boldmathD\unboldmathMinimum\boldmathy\unboldmath∈Ω(\boldmathx\unboldmath,\boldmathd\unboldmath)\boldmath% b\unboldmathT\boldmathy\unboldmath⎞⎟⎠ (1)

subject to

 \boldmathc\unboldmathT\boldmathx% \unboldmath ≤ Π (2) x ∈ {0,1}, (3)

where is the vector of first stage binary variables representing the investment vs no investment in reinforcing or building new lines, is the investment cost vector, is the second stage continuous variables vector representing the random or uncertain parameters, i.e. generation capacities and level of demand, is the uncertainty set, is the vector that includes operating costs, and is the third stage continuous variables vector referring to operating variables. These operating variables include power consumed, power flows, power produced by generating units, load shed by demand and voltage angles at buses (see A for a detailed description of the formulation). represents the maximum budget for investment in transmission expansion. Finally, defines the feasibility region for the operating variables , as a function of investment decisions and given realizations of the uncertain parameters , as follows:

 Ω(\boldmathx\unboldmath,\boldmathd\unboldmath)=⎧⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪⎩\boldmathA\unboldmath\boldmathx% \unboldmath+\boldmathB\unboldmath\boldmathy\unboldmath=\boldmathE\unboldmath:\boldmathλ\unboldmath\boldmathF\unboldmath\boldmathx\unboldmath+% \boldmathG\unboldmath\boldmathy\unboldmath≤% \boldmathK\unboldmath:\boldmathμ\unboldmath\boldmathI\unboldmatheq\boldmathy\unboldmath=\boldmathd\unboldmath:\boldmathα\unboldmath\boldmathI\unboldmathineq\boldmathy\unboldmath≤\boldmathd\unboldmath:\boldmathφ\unboldmath, (4)

where , , , , and are matrices with constant parameters dependent on the network configuration and element characteristics, selects the components of that are equal to the uncertain parameters (demand), and selects the components of that are limited by the uncertain parameters (i.e. maximum power generation and maximum load shedding). The first set of equality constraints correspond to those enforcing the power balance at every bus, the power flow through each line, and fixing the voltage angle of the reference bus. The second 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 reference RuizC:14 ().

### 2.1 The bi-level approach put forward by Jabr

Jabr (Jabr:13, ) proposes tackling the problem (1)-(3) by decomposing and iteratively solving a subproblem and a master problem. The master variables correspond to , i.e. the first stage binary variables vector . Thus, for given values of these master variables, the subproblem corresponds to:

 Maximum\boldmathd\unboldmath∈\boldmathD\unboldmathMinimum\boldmathy\unboldmath∈Ω(\boldmathx\unboldmath,\boldmathd\unboldmath)\boldmath% b\unboldmathT\boldmathy\unboldmath. (5)

Considering that the dual problem associated with the third-level is equal to:

 Maximize\boldmathλ\unboldmath,\boldmathμ\unboldmath,\boldmathα\unboldmath,\boldmathφ\unboldmath(\boldmathE\unboldmath−% \boldmathA\unboldmath\boldmathx\unboldmath)T\boldmath% λ\unboldmath−(\boldmathK\unboldmath−\boldmathF% \unboldmath\boldmathx\unboldmath)T\boldmathμ% \unboldmath+\boldmathd\unboldmathT(\boldmathα% \unboldmath−\boldmathφ\unboldmath) (6)

subject to

 \boldmathB\unboldmathT\boldmathλ% \unboldmath−\boldmathG\unboldmathT\boldmathμ% \unboldmath+\boldmathI\unboldmathTeq\boldmathα\unboldmath−\boldmathI\unboldmathTineq% \boldmathφ\unboldmath = b (7) μ ≥ 0 (8) φ ≥ 0, (9)

it can be substituted into (5), which results in the following single level maximization problem:

 Maximize\boldmathd\unboldmath,\boldmathλ\unboldmath,%\boldmath$μ$\unboldmath,\boldmathα\unboldmath,% \boldmathφ\unboldmathfdual=(\boldmathE\unboldmath−\boldmathA\unboldmath% \boldmathx\unboldmath)T\boldmathλ\unboldmath−(% \boldmathK\unboldmath−\boldmathF\unboldmath\boldmathx% \unboldmath)T\boldmathμ\unboldmath+\boldmathd% \unboldmathT(\boldmathα\unboldmath−\boldmathφ\unboldmath) (10)

subject to (7), (8), (9), and

 \boldmathd\unboldmath∈\boldmathD\unboldmath. (11)

Subproblem (10)-(11) 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 (see (Jabr:13, ) for more details). This simplification consists in assuming that the uncertain parameters are either at the nominal, upper or lower limits of their uncertainty range, as shown by (WuCX:08, ) and (Jabr:13, ).

The optimal solution for the subproblem (10)-(11) provides the dual variable information and the uncertain parameter values to construct Benders cuts for the master problem. The main problem of this type of decomposition algorithm is the slow convergence property, which forced the author to include additional linear constraints obtained from a number of previously computed realizations of the uncertainty vector. These additional cuts accelerate convergence.

### 2.2 Bi-level approach proposed by Chen et al. and Ruiz and Conejo

Chen et al. (ChenWWHW:14, ) and Ruiz and Conejo RuizC:14 () also propose dealing with problem (1)-(3) by decomposing and iteratively solving a subproblem and a master problem. The master variables correspond to , i.e. the vector of first stage binary variables. However, these authors set out a different solution strategy. In particular, the third-level minimization problem is substituted into the subproblem (5) because of its KKT conditions, with which the following single level maximization problem is obtained:

 (12)

subject to

 \boldmathA\unboldmath\boldmathx\unboldmath+\boldmathB\unboldmath\boldmathy\unboldmath = E (13) \boldmathI\unboldmatheq\boldmathy% \unboldmath = d (14) 0 = \boldmathb\unboldmath+\boldmathB\unboldmathT\boldmathλ\unboldmath−\boldmathG\unboldmathT\boldmathμ\unboldmath+\boldmathI\unboldmathTeq\boldmathα\unboldmath−\boldmathI\unboldmathTineq\boldmathφ\unboldmath (15) 0 ≤ \boldmathK\unboldmath−\boldmathF\unboldmath\boldmathx\unboldmath−\boldmathG\unboldmath% \boldmathy\unboldmath⊥\boldmathμ\unboldmath≥0 (16) 0 ≤ \boldmathd\unboldmath−\boldmathI\unboldmathineq\boldmathy\unboldmath⊥\boldmathφ% \unboldmath≥0 (17) d ∈ \boldmathD\unboldmath, (18)

where constraint (15) results from differentiating the Lagrangian of the third-level problem with respect to third-level variables , and constraints (16) and (17) represent the complementary conditions associated with inequality constraints from (4).

Subproblem (12)-(18) is a nonlinear mathematical programming problem, which can be linearized and transformed into a mixed-integer linear mathematical programming problem at the expense of using the Fortuny-Amat (FortunyAmatM:81, ) transformation, which requires including a binary variable for each inequality constraint. For more details about Fortuny-Amat transformation related to this specific application, see reference RuizC:14 ().

The optimal solution of subproblem (12)-(18) provides the uncertain parameter values to construct primal cuts for the master problem at iteration , which corresponds to the following optimization problem:

 Minimize\boldmathx\unboldmath,\boldmathy\unboldmath(i);∀i=1,…,k−1\boldmathc% \unboldmathT\boldmathx\unboldmath+γ (19)

subject to

 γ ≥ \boldmathb\unboldmathT\boldmathy% \unboldmath(i);∀i=1,…,k−1 (20) γ ≥ 0 (21) \boldmathA\unboldmath\boldmathx\unboldmath+\boldmathB\unboldmath\boldmathy\unboldmath(i) = \boldmathE\unboldmath;∀i=1,…,k−1 (22) \boldmathF\unboldmath\boldmathx\unboldmath+\boldmathG\unboldmath\boldmathy\unboldmath(i) ≤ \boldmathK\unboldmath;∀i=1,…,k−1 (23) \boldmathI\unboldmatheq\boldmathy% \unboldmath(i) = \boldmathd\unboldmath(i);∀i=1,…,k−1 (24) \boldmathI\unboldmathineq\boldmathy% \unboldmath(i) ≤ \boldmathd\unboldmath(i);∀i=1,…,k−1 (25) \boldmathc\unboldmathT\boldmathx% \unboldmath ≤ Π (26) x ∈ {0,1}. (27)

This master problem, besides variable , includes one operating variable vector for each realization of the uncertain parameters obtained from subproblem (12)-(18). Unlike the Benders master problem set out in Jabr:13 (), this algorithm converges in a few iterations without needing to include additional primal cuts besides those included by the iterative process itself. Note that the additional constraints included by Jabr to improve convergence of the Benders algorithm, obtained from a prespecified number of previously computed realizations of the uncertainty vector, is a similar strategy with respect to the column-and constrained generation method. However, in the latter case, those realizations are obtained by the process itself, which speeds up convergence.

## 3 Proposed solution technique

In this section the solution technique is set out and justified in detail. However, before focusing on the mathematical description of the solution method, we will justify the uncertainty set type used in this paper.

### 3.1 Uncertainty set

To deal with uncertainty we use cardinality constrained uncertainty sets (BertsimasS:04, ). In particular, we use the same definition of uncertainty set as that given in (Jabr:13, ) but distributed by type of random parameter (generation capacity or demand ) and per region :

 \boldmathD\unboldmath=⎧⎨⎩dGi=¯dGi+^dGizGi;∀i∈ΩGdDi=¯dDi+^dDizDi;∀i∈ΩD (28)

where is the uncertain generation limit related to the th variable within vector , is the corresponding nominal value, is the maximum positive distance from the nominal value that can take the random parameter, and is the set of indices of the generating units. Analogously, , , and correspond to the same values but for demand. Additionally, the following constraints associated with budget uncertainty must hold:

 ∑i∈(ΩG∩Ωr)|zGi| ≤ ΓGr;∀r (29) ∑i∈(ΩD∩Ωr)|zDi| ≤ ΓDr;∀r, (30)

where is the set of indices for region , and are auxiliary continuous variables with the following characteristics and , and and are the maximum number of random parameters for generation capacity and demand, respectively, which may reach their lower or upper limits within region . We include the discrimination by region according to that set out by (RuizC:14, ).

It must be stressed that the auxiliary variables and are initially assumed to be continuous. However, in order to solve the robust TNEP such as by (Jabr:13, ), these variables will be considered as binary. The only limitation introduced by this simplification is that uncertainty budgets and must be integer values, which in our opinion does not detract from the benefits of robust optimization from a practical perspective.

In addition to the simplification proposed by Jabr (Jabr:13, ), we also consider the discovery of Ruiz and Conejo (RuizC:14, ) consisting of the worst outcome in “nature” is for generation capacities to be as low as possible whilst demand loads is as high as possible. According to this, our uncertainty sets are finally defined as follows:

 dGi = ¯dGi−^dGizGi;∀i∈ΩG (31) dDi = ¯dDi+^dDizDi;∀i∈ΩD (32) ∑i∈(ΩG∩Ωr)zGi ≤ ΓGr;∀r (33) ∑i∈(ΩD∩Ωr)zDi ≤ ΓDr;∀r (34) zGi ∈ {0,1};∀i∈ΩG (35) zDi ∈ {0,1};∀i∈ΩD. (36)

Alternatively, instead of uncertainty budgets associated with generation, demand and regions and , a unique uncertainty for each region , or for the system could be used instead.

Note that we do not use the uncertainty set defined by Ruiz and Conejo (RuizC:14, ) because in its current form it does not provide appropriate interpretation with respect to robustness. The definition of the uncertainty set given by (RuizC:14, ) is:

 ∑i∈Ωr(dGmaxi−dGi)∑i∈Ωr(dGmaxi−dGmini) ≤ ΓGr;∀r (37) ∑i∈Ωr(dDi−dDmini)∑i∈Ωr(dDmaxi−dDmini) ≤ ΓDr;∀r. (38)

This formulation presents the following interpretation problem. If or , the uncertainty variables associated with generation take the minimum possible values and those random variables related to load demand take the maximum possible values . This scenario corresponds to the maximum level of uncertainty (Soyster:73, ) as pointed out by (RuizC:14, ). However, if or , the uncertainty variables associated with generation take the maximum values possible and those random variables related to load demand take the minimum possible values . This scenario does not correspond to the no uncertainty case as pointed out in (RuizC:14, ), but rather to the most favourable one from the system operation perspective, and capacity expansion planning using this uncertainty budget would probably result in null investments and the resulting network could be congested for all possible values of the random parameters within the uncertainty set. Because of this it is difficult to interpret intermediate situations and .

### 3.2 Proposed decomposition method

Once the uncertainty set is properly defined, we focus on our methodological proposal for solving problem (1)-(3). It is based on the following observations:

1. The subproblem defined by Jabr (Jabr:13, ) is much more efficient than that proposed by Cheng et al. (ChenWWHW:14, ) and Ruiz and Coenjo RuizC:14 (). Note that the number of binary variables is equal to the number of uncertain parameters (cardinality of plus cardinality of , i.e. ), which are the only variables involved including the dual variables of the third-level problem and the uncertain parameters themselves. Conversely, the subproblem defined using the KKT conditions of the third level problem includes the uncertain parameters, the primal and dual variables of the third-level problem. The number of binary variables associated with the maximum generation capacities and upper bounds of load-shedding are just equal to the number of uncertain parameters, and we have to include two additional binary variables per limit associated with the following variables: line flow, voltage angles, power generation and load-shedding. This results in a much bigger and complex subproblem to solve.

2. Furthermore, the only variables required for the master problem set out by (ChenWWHW:14, ) and RuizC:14 () are the uncertain parameters. Therefore, calculation of the primal variables within the subproblem is useless and inefficient.

3. The master problem proposed by (ChenWWHW:14, ) and RuizC:14 () based on primal cuts is at an advantage in computing terms with respect to the Benders master problem based on dual cuts set out by (Jabr:13, ). It converges in a few iterations and does not require including additional cuts artificially.

For the reasons given above, we use as a subproblem a slight modification of that proposed by Jabr (Jabr:13, ) and as a master problem that set out by (ChenWWHW:14, ) and RuizC:14 (), as follows:

Subproblem:

The subproblem is aimed at obtaining for given values of the investment variables, i.e. for a given network configuration, the worst outcome for uncertain parameters producing the maximum optimum operational cost. In terms of equations, it consists of the maximization of the objective function (10) subject to constraints (7), (8), (9) and the uncertainty set definition (31)-(36).

Master problem:

Given the worst outcome for the uncertain parameters from subproblems, the master problem is aimed at selecting the investment, i.e. the network configuration, that minimizes both the investment and the least desirable (maximum) operational costs associated with each set of values of the uncertain parameters. In terms of equations, it consists of minimizing the objective function (19) subject to constraints (20)-(27).

The only additional detail required to properly define the method is the linearization of the bilinear term included in the objective function of the subproblem, i.e. . Considering equations (31)-(32), this bilinear term becomes:

 \boldmathd\unboldmathT(\boldmathα\unboldmath−\boldmathφ\unboldmath)=∑∀i∈ΩGdGiφGi+∑∀i∈ΩDeidDiφDi+∑∀i∈ΩDdDiαDi=∑∀i∈ΩG(¯dGi−^dGizGi)φGi+∑∀i∈ΩDei(¯dDi+^dDizDi)φDi+∑∀i∈ΩD(¯dDi+^dDizDi)αDi=∑∀i∈ΩG(¯dGiφGi−^dGizGiφGi)+∑∀i∈ΩDei(¯dDiφDi+^dDizDiφDi)+∑∀i∈ΩD(¯dDiαDi+^dDizDiαDi). (39)

Note that for questions of consistency, the negative sign on the compact expression for is implicitly included in its components and , which are defined as negative in the detailed formulation given in the A. is the percentage of load shed by demand .

The terms to be linearized are , , and , which can be replaced by the following continuous variables , , and , respectively, iff the following set of constraints is included:

 φG,minizGi ≤ zGφi≤φG,maxizGi;∀i∈ΩG (40) φG,mini(1−zGi) ≤ φGi−zGφi≤φG,maxi(1−zGi);∀i∈ΩG (41) φD,minizDi ≤ zDφi≤φD,maxizDi;∀i∈ΩD (42) φD,mini(1−zDi) ≤ φDi−zDφi≤φD,maxi(1−zDi);∀i∈ΩD (43) αminizDi ≤ zDαi≤αmaxizDi;∀i∈ΩD (44) αmini(1−zDi) ≤ αDi−zDαi≤αmaxi(1−zDi);∀i∈ΩD, (45)

where parameters , , , , and are the lower and upper limits of the dual variables, which can be replaced by and , respectively, with being a positive and large enough constant. Note that in comparison with the subproblem proposed by Jabr in (Jabr:13, ), we take advantage of the fact that we know in advance what uncertain variables tend, respectively, to the upper and lower limit of the uncertainty set. This strategy allows the number of binary variables to be halved.

Finally, the iterative scheme set out is described step by step on the following algorithm:

Iterative method

Input:

Selection of uncertainty budgets and for each region and the tolerance of the process. This data is selected by the decision maker.

Step 1: Initialization.

Initialize the iteration counter to , and upper and lower limits of the objective function and .

Step 2: Solving the master problem at iteration .

Solve the master problem (19) subject to constraints (20)-(27). The result provides values of the decision variables and . Update the lower limit of the optimal objective function . Note that at the first iteration the optimal solution corresponds to the no investment scenario, alternatively, we could start with any other vector for decision variables.

Step 3: Solving subproblem at iteration .

For given values of the decision variables , we calculate the least optimum operating costs within the uncertainty set , also obtaining the corresponding uncertain parameters . This is achieved by solving problem (10) subject to constraints (7), (8), (9) and the uncertainty set definition (31)-(36).

Update the upper limit of the optimal objective function upper bound .

Step 4: Convergence checking.

If go to Step 5, else update the iteration counter and continue in Step 2.

Step 5: Output.

The solution for a given tolerance corresponds to .

### 3.3 Difference in computational complexity

To show the difference in computational complexity associated with the subproblems set out by Jabr (Jabr:13, ), by Cheng et al. (ChenWWHW:14, ) and Ruiz and Conejo (RuizC:14, ), and that set out in this paper (10)-(11), we calculate the number of variables and constraints required for each formulation in terms of the number of generators (), demand (), buses () and branches ():

Jabr subproblem:

In this case the number of variables and constraints is equal to:

1. binary variables associated with uncertain parameters. Note that each pair of binary variables corresponds to values above and below the nominal value, respectively.

2. continuous variables associated with the dual of the third-level problem, which optimizes operational cost.

3. constraints associated with the third-level dual problem, objective function, robust budgets related to demand and maximum power generation, and the linearization formulas.

Cheng et al. and Ruiz and Conejo subproblem:

In this case the number of variables and constraints is equal to:

1. binary variables related to the Fortuny-Amat transformation.

2. continuous variables associated with the primal and dual of the third-level problem, which optimizes operational costs.

3. constraints related to the third-level primal and dual problem, objective function, and robust budgets related to demand and maximum power generation.

4. constraints associated with the Fortuny-Amat transformation.

Proposed subproblem:

In this case the number of variables and constraints is equal to:

1. binary variables associated with uncertain parameters.

2. continuous variables associated with the dual of the third-level problem, which optimizes operational cost.

3. constraints associated with the third-level dual problem, objective function, robust budgets related to demand and maximum power generation, and the linearization formulas.

According to these figures, and considering that the higher number of elements in real networks are lines and buses, it is clear that the approach using the KKT conditions increases the number of binary variables with respect to other methods by almost doubling the number of lines and buses. Moreover, the number of constraints is also considerably higher. This is why the computational cost for subproblems increases substantially. If we compare the formulation proposed by Jabar and the proposal in this paper, the number of variables reduces by half for the latter, thereby reducing the computational time required to solve the proposed formulation.

Note that we do not compare computational complexity for master problems because it has been demonstrated in the current literature that the column and constraint generation method is more efficient than Benders decomposition.

## 4 Numerical case studies

In this section, we set out numerical experiments of our model and compare them with the method proposed in (ChenWWHW:14, ) and RuizC:14 (). We use an illustrative example i.e. the Garver system (Garver:70, ), and two realistic case studies: the IEEE 24-bus system (IEEE_RTS:99, ) and the IEEE 118-bus test system PSTCA:13 (). It must be stressed that in this paper we just focus on numerical performance as interpretation of robust solutions is beyond its scope but in references Jabr:13 (), (ChenWWHW:14, ) and RuizC:14 ()further details are given.

All examples have been implemented and solved using GAMS (GAMS, ; BrookeKMR:98, ) and CPLEX 12, 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 total running time used by the solvers to work out the masters and subproblems until the final solution is attained. We use the same tolerance for all problems that are equal to .

### 4.1 Illustrative example. Garver system

The model set out is illustrated with the Garver 6-bus system, depicted in Figure 1. This system is made up of 6 buses, 3 generators, 5 types of inelastic demand and 6 lines. Nominal values of generation capacities and demand and their prices and nominal costs can be found in Table 1. The load-shedding cost is equal to the nominal cost of each demand type.

Let us consider that a maximum of three lines can be installed between each pair of buses. Line data is obtained from Table I of GarcesCGR:09 () including construction costs, and the maximum funds available for investment is 40 million euros.

The return on investment period for each line is considered to be 25 years, and the discount rate is 10%, which results in an annual amortization rate of 11%. Since the investment cost is annualized, the weighted factor is equal to the number of hours in a year, i.e. 8760, by which an annualized load-shedding and power generation cost can be obtained.

Regarding the uncertainty set, power generation capacity can oscillate by a maximum of 50% of their nominal values, while demand can change up to a maximum of 20%.

We have solved the robust TNEP problem for three different combinations of uncertainty budgets associated with generation capacities and demand. We have used both the method proposed by (ChenWWHW:14, ) and RuizC:14 () and that set out in this paper. 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. The solution times from these numerical experiments are given in Table 2.

From Table 2 the following observations are pertinent:

1. The method set out is computationally faster on average than that put forward by (ChenWWHW:14, ) and RuizC:14 () for all situations, from the least desirable case which corresponds to Soyster’s solution (Soyster:73, ) through to the deterministic case when uncertainty budgets are null. Regarding the standard deviation of computing times the values are comparable. Computational time in the method set out in this paper is approximately 30% less than the other method with which it has been compared.

2. Both methods provide the same solution and converge in the same number of iterations. Note that if the uncertainty budgets are integer, both approaches are equivalent because the master problem is the same and the subproblems are equivalent.

3. The maximum number of iterations required is four.

The difference in computing time is associated with the solution to the subproblem. That set out by (ChenWWHW:14, ) and RuizC:14 () has 472 equations, 245 continuous and 119 binary variables, while the subproblem set out in this paper contains 119 equations, 187 continuous and 8 binary variables. Even though we are dealing with a small example, the difference in complexity between both subproblems is considerable, especially in terms of binary variables.

### 4.2 IEEE 24-bus Reliability Test System

The following case study is based on the IEEE 24-bus Reliability Test System (RTS) (IEEE_RTS:99, ), depicted in Fig. 2. The system comprises 24 buses, 34 existing corridors which can accept a maximum of three equal lines and 7 new corridors, generating units and loads. Data for lines in existing corridors has been taken from (IEEE_RTS:99, ), while line data for new corridors has been obtained from Table I in (FangH:03, ). Investment costs are 20 million euros, and analogously to the Garver illustrative example, the return on investment period for each line is considered to be 25 years, and the discount rate is 10%, which results in an annual amortization rate of 10%.

Table 3 shows where generators and demand can be located in the network, the maximum power and the corresponding costs for generating units, and the load and the nominal load-shedding cost for demand levels. The load-shedding cost is equal to 100 times the nominal value of each level of demand.

Regarding the uncertainty set, power generation capacity can oscillate by a maximum of 50% of its nominal values, while demand can change by up to a maximum of 20%.

We have solved the same computational experiment as in the previous illustrative example, also comparing it with the method set out by (ChenWWHW:14, ) and RuizC:14 (). Solution times from these numerical experiments are given in Table 4.

From Table 4 we can extract the same conclusions as in the previous illustrative example:

1. The method set out is considerably faster in terms of average computing time than that described by (ChenWWHW:14, ) and RuizC:14 () for every situation, from the least desirable scenario which corresponds to Soyster’s solution (Soyster:73, ) through to the deterministic case, with the latter being that having the lowest computing times. This is because the values of the binary variables of the subproblems are known in advance and only 2 iterations are required. This has made the average computational time 235 times faster.

2. The maximum number of iterations required is three.

The difference in computing time is associated with the solution to the subproblem. That set out by (ChenWWHW:14, ) and RuizC:14 () has 1397 equations, 757 continuous and 346 binary variables, while the subproblem set out in this paper contains 370 equations, 556 continuous and 27 binary variables. The difference in complexity between both subproblems is considerable.

### 4.3 IEEE 118-bus test system

FFinally, we ran additional computational tests using the IEEE 118-bus test system (PSTCA:13, ). The system was made up of 118 buses, 186 existing lines, generating units and loads. In addition, it was 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 was taken from (PSTCA:13, ). The investment cost came to 100 million euros, and as in both previous examples, the return on investment period for each line was considered to be 25 years, and the discount rate was 10%, which resulted in a 10% annual amortization rate.

Table 6 in B shows where the generators and demand is located in the network, the maximum power and the corresponding cost for generating units, and the load and nominal load-shedding cost for demand levels. The load-shedding cost is equal to 1.2 times the nominal load-shedding cost for each level of demand.

Regarding the uncertainty set, power generation capacity and demand can oscillate up to a maximum of 50% of their nominal values.

We have solved the same computational experiment as in the previous illustrative example, and have compared it with the method set out by (ChenWWHW:14, ) and RuizC:14 (). The solution times from these numerical experiments are given in Table 5.

From Table 5 the following observations are pertinent:

1. The method set out in this paper is far faster in average computing time than that proposed in (ChenWWHW:14, ) and RuizC:14 () for every situation. In fact, with the latter method the optimal solution was not reached within the set time limit of 36000 s (10 hours). Note that if any of the solvers (for master or subproblems) does not find an optimal solution within that time the process is stopped. Conversely, with the method set out in this paper the optimal solution was reached in less than one minute for all instances of budget uncertainties.

2. The maximum number of iterations required was 5.

The difference in computing time is associated with the solution to the subproblem. That set out by (ChenWWHW:14, ) and RuizC:14 () has 4115 equations, 2367 continuous and 1018 binary variables, while the subproblem set out in this paper contains 1548 equations, 1712 continuous and 145 binary variables. The difference in complexity between both subproblems is apparent and explains the large difference in computing time between both methods.

### 4.4 Discussion

Numerical experiments demonstrate that the method set out in this paper is computationally more efficient than that by (ChenWWHW:14, ) and RuizC:14 (). Nevertheless, this result is evident considering the difference in complexity of the subproblem. Note that Ruiz and Conejo RuizC:14 () already claimed that the most time-consuming stage of their process was the subproblem resolution, and it is the subproblems solved by our method that really makes a difference in terms of saving computational time, while master problems remain unchanged.

We have not reported results with respect to the method proposed by Jabr (Jabr:13, ). Nevertheless, it is also evident that our method is more efficient because we use the same subproblem whilst halving the number of binary variables. This difference might imply large computing time differences for large problems. In fact, the larger the problem, the larger the difference which is due to how standard branch-and-cut solvers function. The only uncertainty as to whether our method is faster might stem from the master problem. However, the number of iterations reported in this paper and in references (ChenWWHW:14, ) and RuizC:14 () are always lower than or equal to 5, while the number of iterations reported by (Jabr:13, ) varied considerably depending if additional cuts were included or not, but they were usually above 5 or more. Furthermore, this strategy of including additional primal cuts could be used in our master problem by just considering additional outcomes for uncertain parameters as Jabr (Jabr:13, ) does in his work. Nevertheless, because of the convergence behavior this is not necessary.

## 5 Conclusions

In this paper a new decomposition algorithm has been set out to attain the exact solution for the TNEP problem stemming from the use of a two-stage adaptive robust strategy. Although there is nothing novel with respect to the formulation of the problem, we have managed to combine formulations and findings from different authors in the proper manner. Therefore, in terms of computational efficiency, our method leads, occasionally by a large margin, than other methodologies in use.

We have run several computational experiments to show how the algorithm set out in this paper functions with respect to the approach set out by (ChenWWHW:14, ) and RuizC:14 (). Nevertheless, it has been demonstrated that the computational time required by this algorithm is always lower than those set out by (Jabr:13, ), (ChenWWHW:14, ) and (RuizC:14, ). The subproblem set out in this paper is less complex and has fewer binary variables than as is the case in other research, while the primal master problem is equal to that set out by (ChenWWHW:14, ) and (RuizC:14, ), which is at an advantage in terms of computing with respect to the Benders master problem used in (Jabr:13, ), and, moreover, the inclusion of additional cuts from given realizations of the random parameters is not required.

Bearing in mind that with the method set out all features as regards robust optimization design with respect to their predecessors were maintained, the algorithm set out in this paper turned out to be the most efficient method for solving the robust transmission network expansion planning to date, and takes advantage of state-of-the-art mixed-integer mathematical programming solvers. The only limitation is the integrality constraint of budget uncertainty, which simplifies the subproblem considerably although without detracting from the benefits of using robust optimization.

## Acknowledgements

This work set out in this paper has been partially funded by projects by the Junta de Comunidades of 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.

## Appendix A Detailed formulation of the TNEP problem

Before describing the detailed formulation, the main notation used and not defined previously is stated below for quick reference.

Constants:

Susceptance of line (S).

Generation cost for generator (€/MWh).

Load-shedding cost for consumer (€/MWh).

Investment cost of building line (€).

Percentage of load shed by the -th demand.

Capacity of line (MW).

Sending-end bus of line .

Receiving-end bus of line .

Primal variables:

Power consumed by the -th demand (MW).

Power flow through line (MW).

Power produced by the -th generating unit (MW).

Load shed by the -th demand (MW).

Binary variable that is equal to if line is built and otherwise.

Voltage angle at bus (radians).

Dual variables: Note that dual variables are provided after the corresponding equalities or inequalities separated by a colon.

Indices and Sets:

Bus index where the -th generating unit is located.

Bus index where the -th demand is located.

Set of indices of the demand located at bus .

Set of indices of the generating units located at bus .

Set of indices of the demand.

Set of indices of the generating units.

Set of all prospective and existing transmission lines.

Set of all prospective transmission lines.

Set of all networks buses.

TNEP detailed problem:

The decision-making problem pertaining to a transmission planner simultaneously minimizes network investment, generation and load-shedding costs. Problem (1)-(3) for given values of the random parameters is as follows:

 Minimize\boldmathx\unboldmath,\boldmathy\unboldmathR⎛⎜⎝∑k∈ΩL+ckxk⎞⎟⎠+σ⎛⎜⎝∑i∈ΩGcGigi+∑j∈ΩDcUjrj⎞⎟⎠ (46)

subject to

 xk = 1;∀k∈ΩL∖ΩL+ (47) xk ∈ {0,1};∀k∈ΩL (48) ∑k∈ΩL+ckxk ≤ Π (49) ∑i∈ΨGsgi−∑k∣o(k)=sfk +∑k∣r(k)=sfk+∑j∈ΨDsrj = ∑j∈ΨDsdj:λs;∀s∈ΩN (50) fk = (51) θs = 0:χs;s:slack (52)
 fk ≤ fmaxk:ϕmaxk;∀k∈ΩL (54) fk ≥ −fmaxk:ϕmink;∀k∈ΩL (55) θs ≤ π:ξmaxs;∀s∈ΩN∖s:slack (56) θs ≥ −π:ξmins;∀s∈ΩN∖s:slack (57) gi ≥ 0;∀i∈ΩG (58) rj ≥ 0;∀j∈ΩD (59) dj = dDj:αDj;∀j∈ΩD (60) gi ≤ dGi:φGi;∀i∈ΩG (61) rj ≤ ejdDj:φDj;∀j∈ΩD. (62)

The generation and load-shedding costs are multiplied by the weighting factor to make the annual investment cost and the generation and load-shedding cost comparable quantities. Investment costs in (46) are multiplied by the capital recovery factor , calculated as , where is the interest rate and is the number of years considered. Note that dual variables are provided after the corresponding equalities or inequalities separated by a colon. Constraint (49) sets an upper limit on the investment cost and corresponds to (2). Constraints (50) set the power balance at every bus. Constraints (51) represent the power flow through each line. Each of these constraints is multiplied by a binary variable, thus, if the corresponding line is not physically connected to the network, the power flow through it is zero. Equation (52) fixes the voltage angle of the reference bus. Note that (50)-(52) corresponds to the linear constraint set associated with in the first row of (4).

Constraints (54)-(55) set the line flow limits. Constraints (56)-(57) set limits on the voltage angles at every bus, and (58)-(59) ensures power generation and load-shedding are both positive. All these inequality constraints corresponds to the inequality set related to in the second row of (4).

Constraints (60) force demand levels to be equal to the uncertain demand variable, and corresponds to the equality constraint set related to in the third row of (4). Finally, (61) and (62) set power generation and load-shedding to be lower than the generation capacity and a percentage of demand, respectively. These last two sets of constraints correspond to the inequality constraint set related to in the fourth row of (4).

Detailed subproblem:

For the sake of completeness, we also provide the detailed formulation of the proposed subproblem, made up of the objective function (10) subject to constraints (7), (8), (9), the uncertainty set definition (31)-(36) and the linearization (40)-(45). It is as follows:

 (63)

subject to:

 λs(i)+φGi ≤ cGiσ;∀i∈ΩG (64) −λs(j)+αDj ≤ 0;∀j∈ΩD (65) λs(j)+φDi ≤ cUjσ;∀j∈ΩD (66) −λo(k)+λr(k)+ϕk+ϕmaxk+ϕmink = 0;∀k∈ΩL (67) −∑k∣o(k)=sbkxkϕk+∑k∣r(k)=sbkxkϕk +ξmaxs+ξmins = 0∀s∈ΩN∖s:slack (68) −∑k∣o(k)=sbkxkϕk+∑k∣r(k)=sbkxkϕk
 +χs = 0s:slack (69) −∞≤ λs ≤∞;∀s∈ΩN (70) −∞≤ ϕk ≤∞;∀k∈ΩL (71) −∞≤ χs ≤∞;s:slack (72)
 ϕmax