Applications of optimization under uncertainty methods on power system

planning problems

by

Bokan Chen

A dissertation submitted to the graduate faculty

in partial fulfillment of the requirements for the degree of

DOCTOR OF PHILOSOPHY

Major: Industrial Engineering

Program of Study Committee:

Lizhi Wang, Major Professor

Mingyi Hong

James D. McCalley

Sarah M. Ryan

Lily Wang

Iowa State University

Ames, Iowa

2016

Copyright © Bokan Chen, 2016. All rights reserved.

 ii

TABLE OF CONTENTS

LIST OF TABLES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . v

LIST OF FIGURES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vi

ACKNOWLEDGEMENTS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vii

ABSTRACT . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viii

CHAPTER 1. INTRODUCTION . . . . . . . . . . . . . . . . . . . . . . . . . . 1

1.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1

1.2 Optimization Under Uncertainty Methods . . . . . . . . . . . . . . . . . . . . . 1

1.3 Transmission Expansion Planning . . . . . . . . . . . . . . . . . . . . . . . . . . 3

1.4 Power System Resiliency Planning . . . . . . . . . . . . . . . . . . . . . . . . . 5

CHAPTER 2. ROBUST OPTIMIZATION FOR TRANSMISSION EXPAN-

SION PLANNING: MINIMAX COST VS. MINIMAX REGRET . . . . . 7

2.1 Nomenclature . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8

2.2 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9

2.3 Model Formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14

2.3.1 Deterministic model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14

2.3.2 The MMC model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16

2.3.3 The MMR model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16

2.4 Algorithm Development . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17

2.4.1 The master problem for the MMC model . . . . . . . . . . . . . . . . . 18

2.4.2 The subproblem for the MMC model . . . . . . . . . . . . . . . . . . . . 18

2.4.3 Algorithm for the MMC model . . . . . . . . . . . . . . . . . . . . . . . 19

2.4.4 Algorithm for the MMR model . . . . . . . . . . . . . . . . . . . . . . . 20

 iii

2.5 Case Study . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22

2.6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27

CHAPTER 3. ROBUST TRANSMISSION PLANNING UNDER UNCER-

TAIN GENERATION INVESTMENT AND RETIREMENT . . . . . . . . 29

3.1 Nomenclature . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30

3.2 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31

3.3 Model Formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35

3.4 Solution Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38

3.4.1 Master problem and subproblem . . . . . . . . . . . . . . . . . . . . . . 39

3.4.2 Solving the subproblem . . . . . . . . . . . . . . . . . . . . . . . . . . . 40

3.4.3 Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42

3.5 Case Study . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43

3.6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51

CHAPTER 4. A STOCHASTIC PROGRAMMING FRAMEWORK FOR

OPTIMAL DISTRIBUTION NETWORK HARDENING . . . . . . . . . . 54

4.1 Nomenclature . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55

4.2 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56

4.3 Model Formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60

4.4 Solution Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62

4.4.1 Sample average approximation . . . . . . . . . . . . . . . . . . . . . . . 62

4.4.2 Single replication procedure . . . . . . . . . . . . . . . . . . . . . . . . . 63

4.4.3 L-shaped method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64

4.4.4 Greedy search heuristic . . . . . . . . . . . . . . . . . . . . . . . . . . . 66

4.5 Case Study . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67

4.5.1 Scenario generation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67

4.5.2 Experiment results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69

4.5.3 Solution validation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71

4.6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72
 iv

CHAPTER 5. CONCLUSION AND DISCUSSION . . . . . . . . . . . . . . . 75

BIBLIOGRAPHY . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78
 v

LIST OF TABLES

Table 2.1 Motivating Example for the Minimax Regret Model . . . . . . . . . . . 12

Table 2.2 Motivating Example for the Minimax Cost Model . . . . . . . . . . . . 12

Table 2.3 Generation Parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . 22

Table 2.4 Candidate Line Parameters . . . . . . . . . . . . . . . . . . . . . . . . 23

Table 2.5 Number of Iterations for Each Instance . . . . . . . . . . . . . . . . . . 24

Table 2.6 Expansion plan of the MMC approach . . . . . . . . . . . . . . . . . . 25

Table 2.7 Expansion plan of the MMR approach . . . . . . . . . . . . . . . . . . 26

Table 2.8 Comparison of the MMC, MMR and Deterministic Solutions for Uncer-

tainty Set U1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26

Table 2.9 Comparison of the MMC, MMR and Deterministic Solutions for Uncer-

tainty Set U2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27

Table 2.10 Investment and Operational costs for scenarios ($m) . . . . . . . . . . 27

Table 3.1 Comparison with other robust optimization transmission planning mod-

els in the literature . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35

Table 4.1 Hardening costs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70

Table 4.2 Avergae Load Shedding under different budgets . . . . . . . . . . . . . 72

Table 4.3 Confidence interval of the optimality gap under different budgets . . . 73
 vi

LIST OF FIGURES

Figure 3.1 Summary of the robust optimization modeling framework . . . . . . . 36

Figure 3.2 Flow chart of the decomposition algorithm . . . . . . . . . . . . . . . . 43

Figure 3.3 WECC 240-bus test system topology from [52] . . . . . . . . . . . . . . 44

Figure 3.4 Visualization of the transmission planning problem . . . . . . . . . . . 47

Figure 3.5 Six transmission plans. . . . . . . . . . . . . . . . . . . . . . . . . . . . 49

Figure 3.6 Twelve scenarios. The label ‘W1’ is for the worst-case scenario for Plan

1, ‘B2’ for the best-case scenario for Plan 2, and so on. . . . . . . . . . 50

Figure 3.7 Performance of six plans. The horizontal axis of each sub-graph is the

20-year planning horizon, containing four discrete 5-year periods. . . . 52

Figure 4.1 Typical Scenario . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69

Figure 4.2 IEEE 33-Bus Test Case . . . . . . . . . . . . . . . . . . . . . . . . . . . 70

Figure 4.3 24-Hour Load Multiplier . . . . . . . . . . . . . . . . . . . . . . . . . . 71

Figure 4.4 Hardening Plan (Budget = 7m$) . . . . . . . . . . . . . . . . . . . . . 72

Figure 4.5 Histogram of Load Shedding Under Different Plans (Budget = 7m$) . 73
 vii

ACKNOWLEDGEMENTS

I would like to express my profound gratitude towards those people who helped me with

various aspects of conducting research and the writing of this dissertation. First and foremost, I

would like to thank my advisor Dr. Lizhi Wang for his guidance, patience and support through-

out my PhD study. He is a great teacher who not only taught me the fundamental technical

skills in operations research, but also guided me towards conducting research. He provided me

with his insights, inspiration and encouragement to help me complete this dissertation. I could

not have done this without him.

I would also like to thank my committee members: Dr. Mingyi Hong, Dr. James D.

McCalley, Dr. Sarah M. Ryan and Dr. Li Wang for their tremendous help. I appreciate their

time and effort in reading this dissertation and their helpful comments in improving this work.

In addition, I would like to express my gratitude towards Dr. Jianhui Wang from Argonne

National Laboratory for his effort and contribution to this work. His advice and input to

my research is invaluable. Needless to say, I am solely responsible for any errors that this

dissertation might contain.

Last but not least, I would like to thank my friends and colleagues for their companionships.

They make the office a fun place to be. I enjoyed our discussions on life, the universe, and

everything.
 viii

ABSTRACT

This dissertation consists of two published journal paper, both on transmission expansion

planning, and a report on distribution network hardening.

We first discuss our studies of two optimization criteria for the transmission planning prob-

lem with a simplified representation of load and the forecast generation investment additions

within the robust optimization paradigm. The objective is to determine either the minimum of

the maximum investment requirement or the maximum regret with all sources of uncertainty

explicitly represented. In this way, transmission planners can determine optimal planning de-

cisions that are robust against all sources of uncertainty. We use a two layer algorithm to

solve the resulting trilevel optimization problems. We also construct a new robust transmission

planning model that considers generation investment more realistically to improve the quan-

tification and visualization of uncertainty and the impacts of environmental policies. With

this model, we can explore the effect of uncertainty in both the size and the location of can-

didate generation additions. The corresponding algorithm we develop takes advantage of the

structural characteristics of the model so as to obtain a computationally efficient methodology.

The two robust optimization tools provide new capabilities to transmission planners for the

development of strategies that explicitly account for various sources of uncertainty.

We illustrate the application of the two optimization models and solution schemes on a

set of representative case studies. These studies give a good idea of the usefulness of these

tools and show their practical worth in the assessment of “what if” cases. We compare the

performance of the minimax cost approach and the minimax regret approach under different

characterizations of uncertain parameters. In addition, we also present extensive numerical

studies on an IEEE 118-bus test system and the WECC 240-bus system to illustrate the

effectiveness of the proposed decision support methods. The case study results are particularly

useful to understand the impacts of each individual investment plan on the power system’s
 ix

overall transmission adequacy in meeting the demand of the trade with the power output units

without violation of the physical limits of the grid.

In the report on distribution network hardening, a two-stage stochastic optimization model

is proposed. Transmission and distribution networks are essential infrastructures to modern

society. In the United States alone, there are there are more than 200,000 miles of high voltage

transmission lines and numerous distribution lines. The power network spans the whole country.

Such vast networks are vulnerable to disruptions caused by natural disasters. Hardening of

distribution lines could significantly reduce the impact of natural disasters on the operation

of power systems. However, due to the limited budget, it is impossible to upgrade the whole

power network. Thus, intelligent allocation of resources is crucial. Optimal allocation of limited

budget between different hardening methods on different distribution lines is explored.

 1

CHAPTER 1. INTRODUCTION

1.1 Overview

In this proposal, I present research during my PhD career. My research has been focusing

on the application of optimization under uncertainty techniques on problems arise in power

systems. Two main topics are covered in this proposal:

(1) The application of robust optimization on transmission expansion planning problems.

(2) The application of stochastic programming in the planning of the hardening strategies of

distribution networks.

This chapter presents a brief introduction to the discussed topics. Chapter 2 and Chapter 3

contains my two papers on robust transmission expansion planning. My report for the third

part of my research— power system resiliency planning, is presented in Chapter 4.

1.2 Optimization Under Uncertainty Methods

For many optimization problems, the input parameters of a developed model are assumed

to be deterministic. However, in more and more real world applications, this is no longer true.

In such cases, the optimal solution obtained by a deterministic model may not be optimal

or even feasible. As such, methods of optimization under uncertainty have been garnering a

lot of attention in fields including power system. Currently, two of the most popular ways

to deal with uncertainty in mathematical programming is robust optimization and stochastic

programming. Both methods have been widely applied in power system [61].

The problems solved in this dissertation all fall roughly into the category of two-stage

optimization models. The two stages refer to the two decision making process before and after
 2

uncertainties are observed. If no uncertainty is considered, and suppose b denote the sources of

uncertainty, x denote the first-stage decisions that should be made before uncertainty realization

and y denote the second-stage decisions that can only be made after uncertainty is observed, a

typical deterministic linear programming model can be formulated as follows:

minx,y c> x + d> y (1.1)

s. t. x∈X (1.2)

Ax + By ≤ b̄ (1.3)

where b̄ represent the mean of the uncertainty parameters.

In robust optimization, it is usually assumed that uncertain parameters in an optimization

model can be described by polyhedral sets [11, 25]. Such uncertainty sets can be constructed

by using information like the lower and upper bounds of the uncertain parameters. Then

the objective value under the worst-case scenario is optimized. Although robust optimization

is relatively conservative and can lead to worse performance in average, it can guarantee the

feasibility of the optimal solution and provide a bound on the objective value. Typically, robust

optimization can be formulated as follows:

min c> x + maxb∈U miny∈Y(x,b) d> y (1.4)

x∈X

where U is the uncertainty set and Y(x, b) = {y|Ax+By ≤ b} is the feasible set for second-stage

decision variables. From this formulation we can see that the uncertainty parameter is treated

as a variable to identify the worst-case scenario.

On the other hand, stochastic programming assumes the accurate probability distribution

of the uncertain parameters can be obtained. The expected value of the objective function

over the aforementioned distribution is optimized. Scenarios are usually generated based on

the probability distributions to represent the uncertain data. In order to guarantee feasibility

and good performance, a reasonable number of scenarios should be generated, which leads to

large problem size. Many algorithms can be applied to solve such large-scale optimization prob-

lems, including Benders’ decomposition and progressive hedging [43]. Stochastic programming

models can be formulated as follows:

 3

min +E[Q(x, b̃)] (1.5)

x∈§

where

Q(x, b) = min c> y (1.6)

s. t. Ax + By ≤ b̃ (1.7)

Both robust optimization and stochastic programming have their strengths and weaknesses.

In this dissertation they are applied to two different types of problems based on the requirements

of the applications.

1.3 Transmission Expansion Planning

Transmission planning (TP) is one of the key decision processes in the power system produc-

tion pipeline. Electricity transmission networks are responsible for reliably and economically

delivering power from generators to consumers. Thus, a robust and resilient transmission net-

work is essential to the operation of power systems for decades to come. Good transmission

investments have many benefits, including satisfying increasing demand, promoting social wel-

fare, improving system reliability and resource adequacy, etc. Transmission planning is also

very challenging due to various sources of uncertainty that planners need to consider. Besides

uncertainty sources, such as demand variations and renewable energy intermittency faced by

system operators in short-term scheduling, planners also need to take into consideration uncer-

tainty of policy changes, technological advancements and natural disasters [61]. Such sources

of uncertainty cannot be characterized in terms of analytical probability distributions. For

example, to cope with challenges of climate change, the power generation industry is facing

increasing pressure to reduce greenhouse gas emissions. In addition, a large amount of coal

plants are anticipated to retire in response to the implementation of the EPA clean power plan

[19], and their replacement in part by gas-fired plants. Moreover, after the restructuring of the

power system, certain behavior by generation companies introduces additional uncertainty in

power system planning in general and TP in particular. However, research on the impacts of
 4

uncertain future generation developments on TP is scarce in the decentralized decision making

environment in competitive market. As such, the study of this topic is warranted.

Conventional transmission planning methodologies are typically deterministic and some

make use of ad-hoc approaches to deal with uncertainty. These methodologies have served the

industry relatively well in the vertically-structured industry. Present realities brought about,

by the open access regime and the advent of competitive electricity markets, increased wind

and solar resource outputs. Their highly time varying, uncertain and intermittent patterns,

combined with the more dynamically varying and uncertain loads and a world-ranging envi-

ronmental policy initiative have resulted in a more volatile utilization of transmission facilities.

Critically needed is the improved planning methodologies that can effectively accommodate

the realities of the new environment. Transmission planning, by its very nature, is subject to a

wide range of sources of uncertainty that must be taken into account. The new regime indicates

many additional sources of uncertainty and numerous complications that must be considered.

Over a 10-20 year planning horizon, major sources of uncertainty faced by transmission

planning decision makers include load growth, fuel costs, climate change events, atmospheric

conditions, variable generation outputs, generation expansion/retirement events, regulatory

and legislative developments, technology breakthroughs and market outcomes. When the im-

pacts of transmission plans are assessed, such studies must be carefully performed with the

explicit representation of the sources of uncertainty, in light of the overarching consequences on

the power system. These sources of uncertainty may be classified into two categories: aleatoric

and epistemic sources of uncertainty. Sources of aleatoric uncertainty cannot be replaced by

more accurate measurements but they can be statistically quantified. For instance, the coinci-

dental uncertainty with respect to the electrical system state and network topological state can

be both probabilistically and stochastically modeled and quantified if the relevant sets of input

are made available. On the other hand, sources of epistemic uncertainty are due to information

we may in principle know, but do not in actual practice. Often, such sources are also called

sources of systematic uncertainty. For instance, it is hard to assign meaningful probabilities

to sources of uncertainty associated with government policies, technology breakthroughs, and

investment in renewable energy generation.

 5

Electricity demand is expected to grow continually and steadily for decades to come. In the

EIA Annual Energy Outlook 2013, forecasts indicate that electricity consumption will increase

by 28% from 2011 to 2040 at an average of about 1% per annum rate[24]. To accommodate such

growth in electricity demand, new capacity of electricity generation must be built at a pace that

meets or exceeds the growth of demand and retirement of old generators to maintain power

system reliability. The additional demand and capacity require commensurate transmission

capacity.

With the additional sources of uncertainty and complications facing transmission planners

and the need to accommodate increasing demand as well as more variable generation capacity,

there is a critical need for new TP methodologies that can address the aforementioned concerns.

The first two papers in this dissertation address this problem.

1.4 Power System Resiliency Planning

In addition to transmission expansion planning, resilience planning is also an important

step to improve the security of power systems. In the United States, there are more than

200,000 miles of high voltage transmission lines and numerous distribution lines that span the

whole country. Such vast networks are vulnerable to disruptions caused by natural disasters

including seismic activities and extreme climate events. Due to climate change, frequency of

natural disasters such as floods, droughts, hurricanes and other weather activities is expected

to rise. Improving power system security and resilience not only would reduce the probability

of disrupted operations, but also could be instrumental to the recovery of communities from

catastrophes.

Resilience is defined as a combination of two capabilities: the inherent ability of a system

to cope with disruption with the networks topology and redundancies, and its ability to recover

from disruptions. Both abilities need to be considered in order to improve system resilience.

One illustrative example is underground power lines, even though they are invulnerable to

many natural disasters including hurricanes, wild fires and snow storms, they are more difficult

to repair. In a region prone to earthquakes, even with low probability, undergrounding power

lines may not be a good choice.

 6

According to [54], hardening is one of the most effective approaches to improve the resilience

of power grids. Methods of hardening include upgrading power poles and structures with

stronger materials, vegetation management and undergrounding utility lines. In addition to

hardening, investment can also be devoted to reduce response time and expedite infrastructure

recovery after the occurrence of disruptions. Options include response planning, training of

personnel and improved outage notification ability. Such investments are usually constrained

by a limited budget. Thus, hardening of the entire grid simultaneously could be prohibitively

expensive and therefore is unrealistic. It is crucial to allocate the limited resources to the most

vulnerable components or to the power lines with the largest impact if affected. The last paper

uses stochastic programming to solve the budget allocation problem in power system hardening.
 7

CHAPTER 2. ROBUST OPTIMIZATION FOR TRANSMISSION

EXPANSION PLANNING: MINIMAX COST VS. MINIMAX REGRET

Published in IEEE Transactions on Power Systems

Bokan Chen, Jianhui Wang, Lizhi Wang, Yanyi He and Zhaoyu Wang

Abstract

Due to the long planning horizon, transmission expansion planning is typically subjected to a

lot of uncertainties including load growth, renewable energy penetration, policy changes, etc.

In addition, deregulation of the power industry and pressure from climate change introduced

new sources of uncertainties on the generation side of the system. Generation expansion and

retirement become highly uncertain as well. Some of the uncertainties do not have probability

distributions, making it difficult to use stochastic programming. Techniques like robust op-

timization that do not require a probability distribution became desirable. To address these

challenges, we study two optimization criteria for the transmission expansion planning prob-

lem under the robust optimization paradigm, where the maximum cost and maximum regret

of the expansion plan over all uncertainties are minimized respectively. With these models,

our objective is to make planning decisions that are robust against all scenarios. We use a two

layer algorithm to solve the resulting tri-level optimization problems. Then, in our case studies,

we compare the performance of the minimax cost approach and the minimax regret approach

under different characterizations of uncertainties.

 8

2.1 Nomenclature

Sets and indices

U The polyhedron uncertainty set of demand and new generation capacity profile

I Set of nodes

L Set of existing transmission lines

N Set of candidate transmission lines

T Set of years in the planning horizon

M Set of load blocks

K Set of technology types

Parameters

Pi,k,t Capacity of existing generator of technology k at node i at time t

cT
ij,t Cost of building the new transmission line ij at time t

cLi,t Cost of load curtailment at node i at time t

cP
i,k,t Cost of power production of technology k at node i
C
fk,t,m Average capacity factor of generation technology k at year t load block m

Fijmax The maximum power flow on transmission line ij

Bij Susceptance of transmission line ij

M A big constant used to linearize the power flow constraint

d¯i,t,m The average amount of demand at year t load block m at node i

N
P̄i,k,t The average amount of generation expansion of technology k at node i at time t

θmin The lower bound of voltage angles

θmax The upper bound of voltage angles

λ Market interest rate (Inflation included)

N,min
Pi,k,t Minimum amount of new generation at a node
N,max
Pi,k,t Maximum amount of new generation at a node

P N,min Lower bound on the total amount of generation at all the nodes

P N,max Upper bound on the total amount of generation at all the nodes
 9

Decision Variables

xij Binary variables indicating whether a transmission line is built

N
Pi,k,t N is
The amount of new generation capacity of technology k at node i at time t. Pi,k,t

negative in the case of power plant retirement

fij,t,m Power flow from node i to node j at year t load block m

pi,k,t,m Power production of technology k at node i at year t load block m

ri,t,m The amount of load shedding at node i at year t load block m

θi,t,m Voltage angle at node i at year t load block m

di,t,m Demand at year t load block m at node i

2.2 Introduction

Transmission expansion planning is very challenging due to various uncertainties planners

need to consider. Besides typical high-frequency uncertainties including demand variations

and renewable energy intermittency faced by system operators in short-term scheduling, plan-

ners also need to take into consideration low-frequency uncertain events like policy changes,

technological advancements, natural disasters, etc. [61]. Such uncertainties cannot be charac-

terized by probability distributions. For example, to cope with challenges of climate change,

the power generation industry is facing increasing pressure to reduce greenhouse gas emissions.

In addition, a large amount of coal plants are anticipated to retire in response to government

regulations and fuel price changes [19], replaced partially by gas-fired plants. Many of such

retirements could be announced on relative short notice and unexpected by the system oper-

ator. Moreover, after the deregulation of the power system, strategic behavior of generation

companies in generation expansion becomes an uncertain factor as well.

Currently, the most common practices in dealing with uncertainties in optimization include

stochastic programming and robust optimization. In stochastic programming, scenarios are

generated based on a certain probability distribution of the uncertain data. The weighted sum of

the total cost under different scenarios is usually optimized. Stochastic programming has been

successfully applied to power system capacity expansion planning problems. In [59, 36, 44, 26],

uncertainties including load prediction inaccuracies, transmission line and generator outages,
 10

generation and transmission line capacity factors are considered. All of those uncertainties

can be described by probability distributions and can be effectively modeled with stochastic

programming techniques. However, most of those works only focus on uncertainties in the

operation phase and do not address uncertainties in the planning phase. The reason is that it

is difficult to obtain probability distributions of the non-random [17] or epistemic uncertainties

caused by factors including policy changes, investment behavior of market players, etc. In this

paper, besides load uncertainty, we also take into consideration uncertain generation expansion

behavior of generating companies and coal power plants retirement.

As an alternative tool to address uncertainties, robust optimization [11, 25, 14] can avoid

some of the difficulties arising from stochastic programming approaches. With robust optimiza-

tion, uncertainty is described by parametric sets, which contain an infinite number of scenarios.

Such uncertainty sets can be constructed with information like the lower and upper bounds

of a random variable, which are much easier to derive than probability distributions. This

approach can identify a set of decisions that is robust under the worst-case scenario contained

in the uncertainty sets, which is desirable in planning problems where reliability is important.

The conservativeness of the solution can be adjusted by changing the uncertainty sets [13], de-

pending on how much uncertainty is desired to capture. Robust optimization has been applied

to many problems in power systems. In [62], robust unit commitment with the n − k security

criterion is studied. The problem is then reformulated into a single-level problem with the help

of dual variables. In [12, 32], the two-stage robust unit commitment problems under uncer-

tainty are studied. Both papers propose to use Bender’s Decomposition to solve the problem.

A similar model is used in [74] where demand response is considered. A robust minimax regret

model is proposed in [31] to solve the unit commitment problem under uncertainty. However,

all those works study operation problems. In long term power system planning problems where

more uncertainties need to be considered, application of robust optimization is limited. In [29],

a robust transmission expansion planning model is proposed considering load and generation

uncertainty. Bender’s decomposition is also used.

In this paper, we propose two robust optimization models to address two main sources of

uncertainty: load and generation expansion behavior of generating companies. Two criteria,
 11

minimax cost (MMC) and minimax regret (MMR), are used as the objective of our models.

The MMC criterion has been used widely in robust optimization applications [12]. The MMR

criterion is considered in [31] for the unit commitment problem. In comparison with the MMC

criterion, it is concluded that MMR outperforms MMC for certain unit commitment problems.

However, the same conclusion may not apply to transmission expansion planning problems due

to the different structures of such problems. In [46], regret is considered as one of the objectives

in a multi-objective optimization framework. It is applied to handle non-random uncertainties

in [23, 7]. Both criteria use the performance of a decision under the worst possible scenario as

the objective for optimization, but their main difference is how the “worst scenario” is defined.

The MMC criterion focuses on the cost associated with a decision under a scenario, so the

scenario that results in the highest cost is identified as the worst scenario. On the other hand,

the MMR criterion defines the worst scenario as the one that leads to the highest regret for

the decision maker. For a given decision d0 and a given scenario s0 , the regret is the highest

potential cost savings had the decision maker known that scenario s0 would occur and made a

decision accordingly. More rigorously,

R(d0 , s0 ) = C(d0 , s0 ) − min C(d, s0 ),

d∈D

where C(d0 , s0 ) is the cost associated with d0 and s0 , D is the set of all feasible decisions, and

R(d0 , s0 ) is the regret associated with d0 and s0 . Using these notations, the MMC and MMR

criteria can be respectively formulated as

 
min max C(d, s) and
d∈D s∈U
    
0
min max R(d, s) = min max C(d, s) − min
0
C(d , s) .
d∈D s∈U d∈D s∈U d ∈D

We use two simple examples in Tables 2.1 and 2.2 to demonstrate the differences between

MMC and MMR. In the first example, under MMC, D2 is the optimal decision because its

worst scenario cost, $8, is lower than that of D1 , $9. Under MMR, D1 is the optimal decision

because its worst scenario regret, $1, is lower than that of D2 , $5. The argument for MMR
 12

is that since scenario S1 is a “bad” scenario anyway because D1 and D2 both lead to higher

costs in S1 than S2 , the difference between the costs associated with the two decisions, which

is the regret, may provide more information for decision making than the absolute value of the

cost itself. In the second example, decision D3 is obviously a bad choice because of its high

Table 2.1 Motivating Example for the Minimax Regret Model

Cost / Regret Decision D1 Decision D2
Scenario S1 $9 / $1 $8 / $0
Scenario S2 $2 / $0 $7 / $5

cost in scenario S4 . Decision D5 will be selected under MMC because its worst cost, $18, is

lower than that of D3 , $40, and D4 , $19. Under MMR, decision D4 will be selected since its

worst-case regret is $13 while the regret of D5 is $14. We argue the MMC solution D5 is better

in this example because it is only slightly worse than the MMR decision D4 in terms of regret

in scenario S3 only because of the existence of decision D3 , which cannot be selected anyways,

but has a much lower cost in scenario S4 .

Table 2.2 Motivating Example for the Minimax Cost Model

Cost / Regret Decision D3 Decision D4 Decision D5
Scenario S3 $4 / $0 $16 / $12 $18 / $14
Scenario S4 $40 / $34 $19 / $13 $6 / $0

From the previous two examples, we can see that there are no clear cut answers as to which

criterion is superior. Each of them has advantages and disadvantages. The examples above can

shed some light on which criterion may perform better in what situations. In the first example,

both decisions perform better in one scenario and worse in the other. In this case, it makes

sense to use MMR as the criterion because the MMC criterion is too conservative and does not

consider non-extreme scenarios. On the other hand, in the second example, there exists a very

risky decision that performs well under one scenario and extremely poorly under the other. In

a planning problem where risks should be controlled, such decisions are usually not desirable,

but they may affect the maximum regret of other decisions and distort the final decision.

The two-stage structure of our robust optimization models can capture both the planning
 13

and operation stages of the transmission expansion planning problem very well. They can

be formulated as special cases of trilevel optimization problems. However, due to their non-

linear, non-convex structure, they are very difficult to solve. In previous researches [12, 32],

the authors use Bender’s decomposition to reformulate the problem into a master problem

and a bilinear subproblem, which is then solved with outer approximation. However, the

outer approximation approach cannot handle the binary variables in the subproblem when the

MMR criterion is used. In [31], statistical upper bounds are used to complement the outer

approximation approach. We propose a two-layer algorithm where we decompose our problem

into a master problem and a bilevel subproblem. The master problem is updated with a branch

and cut type procedure, where new constraints and variables are iteratively generated and then

solved as a mixed-integer program. This algorithm is a special case of the bilevel optimization

algorithm [67]. Similar algorithms are proposed in [74, 73]. It works faster than the traditional

Bender’s decomposition approach with the use of primal information instead of dual variables.

The subproblem is a mixed-integer bilevel optimization problem, which is more difficult to

solve. In [49], the difficulty of solving a bilevel linear optimization program is discussed and

several heuristics are proposed. We use the Karush-Kuhn-Tucker (KKT) conditions [15] to

reformulate the bilevel problem into a single level problem with complementarity constraints,

which is then reformulated into a mixed-integer programming problem [28].

The contribution of this paper can be summarized as follows. Firstly, we propose two

robust optimization models that use two criteria to assess decisions. Both load uncertainties

and generation expansion uncertainties are considered. Then, effective algorithms are proposed

to solve the resulting trilevel optimization problems. Finally, in our case studies, the two models

and their corresponding algorithms are tested with a modified IEEE 118-bus test system. We

then analyze the results and compare the performances of the expansion plans under different

scenarios.

The rest of the paper is structured as follows. Section 2.3 presents the mathematical

formulation of the transmission expansion planning problems. In section 2.4, we present the

trilevel optimization algorithm. Numerical results are presented in section 2.5. Finally, section

2.6 concludes this paper.

 14

2.3 Model Formulation

Transmission expansion planning problems are usually modeled as two-stage problems to

account for the long planning horizon, where the expansion planning decisions are made in

the first stage when there is limited information on uncertain parameters and the operational

decisions are made in the second stage after uncertainty realizations are observed. In this

section, we first present the deterministic model, and then introduce two robust optimization

models under the MMC and MMR criteria.

2.3.1 Deterministic model

In the deterministic model, consideration of uncertainty is avoided by assuming perfect

information for all parameters. For example, in the following deterministic model, the load is

fixed as d¯ and the new generation capacity is fixed as P̄ N .

X X
min cT
ij,t xij + (1 + λ)t (cP L
i,k,t pi,k,t,m + ci,t ri,t,m ) (2.1)
ij,t i,k,t,m
X X X
s. t. pi,k,t,m + fji,t,m − fij,t,m = d¯i,t,m − ri,t,m , ∀i ∈ I, t ∈ T , m ∈ M (2.2)
k j j
fij,t,m − Bij (θi,t,m − θj,t,m ) − (1 − xij )M ≤ 0, ∀ij ∈ N (2.3)

Bij (θi,t,m − θj,t,m ) − fij,t,m − (1 − xij )M ≤ 0, ∀ij ∈ N (2.4)

fij,t,m = Bij (θi,t,m − θj,t,m ), ∀ij ∈ L (2.5)

fij,t,m ≤ Fijmax xij , ∀ij ∈ N (2.6)

−fij,t,m ≤ Fijmax xij , ∀ij ∈ N (2.7)

fij,t,m ≤ Fijmax , ∀ij ∈ L (2.8)

−fij,t,m ≤ Fijmax , ∀ij ∈ L (2.9)

C N
pi,k,t,m ≤ fk,t,m (Pi,k,t + P̄i,k,t )∀i ∈ I, t ∈ T , m ∈ M, k ∈ K (2.10)

θmin ≤ θi,t,m ≤ θmax , ∀i ∈ I, t ∈ T , m ∈ M (2.11)

x binary (2.12)
 15

The objective function (2.1) is the transmission capital investment cost and total opera-

tional cost (including cost of power production and load shedding) over the planning horizon.

This model is a static model, in which the total operational cost over the planning horizon

is estimated by extrapolating from |T | years. A similar approach has been used by several

other related studies [36, 26, 35]. Constraint (2.2) requires that the net influx at a node

should be equal to the net outflow. Constraints (2.3) and (2.4) are equivalent to the equation

fij,t,m = xij Bij (θi,t,m − θj,t,m ), which is nonlinear and complicates the model. We introduce the

constant M to linearize this equation [35]. When xij = 1, then the two constraints are reduced

to fij,t,m = Bij (θi,t,m − θj,t,m ), where the value of M does not matter. When xij = 0, then we

need to make sure that M is large enough so that no additional constraints are imposed. On

the other hand, if M is too large, it may cause computational difficulties. In our experiments,

we set it to be ten times the largest value of Fijmax . Equation (2.5) calculates the power flow on

existing transmission lines. Constraints (2.6)- (2.9) dictate that the power flow on transmission

lines does not exceed their limits. Constraint (2.10) specifies the generation capacity on each

node. Constraint (2.11) limits the range of phase angles at a node.

To facilitate algorithmic development and simplify the notations, we abstract the determin-

istic model as follows:

minx,z c> x + b> z (2.13)

s. t. Ax + C1 z ≤ g1 (2.14)

B ȳ + C2 z ≤ g2 (2.15)

Jz = d¯ (2.16)

In this more concise abstract formulation, we use x to represent the binary variable indicating

whether or not a transmission line should be built, y to represent the amount of new generation

and z to represent operational variables including power production, phase angles, power flow

and load curtailment. Vectors c and b represent coefficients of variables in the objective function.

Matrices A, B, C1 , C2 , J are the coefficients of variables in the constraints. Vectors g1 , g2 are

the right-hand-side parameters in the constraints. Constraint (2.14) corresponds to equations

(2.3)-(2.11). Constraint (2.15) corresponds to (2.10). Constraint (2.16) corresponds to (2.2).

 16

2.3.2 The MMC model

In the two-stage MMC model, given the first-stage decisions, the second-stage problem

is commonly known as the recourse problem [38], where the optimal operation decisions are

identified. The feasible set of the recourse problem is defined as follows:

Z(x, d, y) = {z : C1 z ≤ g1 − Ax, C2 z ≤ g2 − By, Jz = d}

The uncertainty set is defined as

U = {(d, y) : Q1 d ≤ q1 , Q2 y ≤ q2 }

The matrices Q1 , Q2 are the coefficients of d and y in the uncertainty set. Vectors q1 , q2 are

the right-hand-side parameters. They can contain information including the lower and upper

bounds of the uncertain parameters, the lower and upper bounds of the linear combination

of the uncertain parameters, etc. Such information can be obtained from historical data or

statistical tests on historical data. In this paper we consider both the uncertainty caused by

load forecast and the uncertainty caused by future generation expansion. In addition, other

types of uncertainties can be easily plugged into the model without affecting the algorithm.

The MMC model can be formulated as:

 
> >
min c x + max min b z . (2.17)
x binary (d,y)∈U z∈Z(x,d,y)

2.3.3 The MMR model

Unlike the MMC model, the MMR model aims to minimize the worst-case regret under all

possible scenarios. Before presenting the MMR model, we first define the feasible set of the

perfect information solution G(d, y) and the perfect information cost G(d, y) as follows:

G(d, y) = {(x̂, ẑ) : Ax̂ + C1 ẑ ≤ g1 , C2 ẑ ≤ g2 − By, J ẑ = d}.

n o
G(d, y) = min c> x̂ + b> ẑ .
x̂,ẑ∈G(d,y)
 17

We can see that G(d, y) is only dependent on the uncertain parameters (d, y) and can only

be known after the uncertainty realizations are observed. We call it the perfect information

solution because G(d, y) can only be achieved if perfect information about the uncertainties is

available.

Then we can define the MMR model as follows:

  
> >
min c x + max min b z − G(d, y) . (2.18)
x binary (d,y)∈U z∈Z(x,d,y)

Comparing the MMC model and MMR model side by side, their similarities are very no-

ticeable. The difference between them lies in their definition of the “worst-case scenario”. With

the MMC criterion, the worst-case scenario is defined as the scenario with the highest cost,

while the MMR criterion defines the worst-case scenario where regret is the highest.

To shed more light on which criterion is more appropriate under different situations, we

can classify scenarios into two categories: regretful vs. regretless. We use xC to denote the

MMC solution and xR to denote the MMR solution. R(x, s) is the regret of decision x under

scenario s. If R(xC , s) ≥ R(xR , s), then we call scenario s a regretful scenario for decision

xC . Otherwise, we say it is regretless. When MMR is used, regret is redistributed among the

scenarios. The regretful ones become less regretful and the costs in the regretless scenarios

increase. As an uncertainty set consists of both regretful scenarios and regretless scenarios, it

is unpractical to predict accurately which type of scenario will occur in the future. However,

the classification of scenarios compares and illustrates the advantages of the MMR and MMC

approaches for decision-makers to choose the criterion more appropriately for their specific

problem. For example, if they are more confident that regretful scenarios will occur, they can

choose the MMR criterion. Otherwise, they may choose the MMC criterion.

2.4 Algorithm Development

In this section, we develop a new trilevel optimization algorithm to solve the two robust

optimization problems, which we decompose into two levels: the master problem and the

subproblem. We first present the algorithm for the MMC model. This algorithm is then

modified for the MMR model. Since we use a cutting plane procedure that does not require
 18

duality information, we can reformulate the sub-problem as a mixed-integer linear programming

problem. In [31], the worst-case scenarios are identified via statistical upper bounds with

Monte Carlo simulation. In contrast, our algorithm provides a theoretical global optimality

guarantee to find the worst-case scenarios as the entire problem is solved as a mixed-integer

linear programming problem after reformulating the sub-problem.

2.4.1 The master problem for the MMC model

The master problem is designed to provide a relaxation of the MMC model (2.17), in which

the search for the worst-case scenario is restricted to be within a given finite set of scenarios,

ΩC = {(di , y i ), ∀i = 1, ..., |ΩC |}, rather than the complete set of scenarios, U. As such, the

master problem yields a lower bound of the MMC model (2.17). We denote the master problem

as MC (ΩC ), and it is formulated as the following single level mixed integer linear program.

min c> x + ξ (2.19)

x,ξ,z i

s. t. ξ ≥ b> z i ∀i = 1, . . . , |ΩC | (2.20)

Ax + C1 z i ≤ g1 ∀i = 1, . . . , |ΩC | (2.21)

By i + C2 z i ≤ g2 ∀i = 1, . . . , |ΩC | (2.22)

Jz i = di ∀i = 1, . . . , |ΩC | (2.23)

x binary. (2.24)

2.4.2 The subproblem for the MMC model

The subproblem is defined as the MMC model (2.17) with a given first-stage decision, x.

As such, the subproblem yields an upper bound of the MMC model (2.17). We denote the

subproblem as SC (x), and it is formulated as the following bilevel linear program.

max min b> z. (2.25)

(d,y)∈U z∈Z(x,d,y)
 19

This model can be further reformulated as the following linear program with complementarity

constraints (LPCC).

max b> z (2.26)

d,y,z,α,β,γ

s. t. Q1 d ≤ q1 (2.27)

Q2 y ≤ q2 (2.28)

0 ≤ g1 − Ax − C1 z ⊥ α ≥ 0 (2.29)

0 ≤ g2 − By − C2 z ⊥ β ≥ 0 (2.30)

Jz = d (2.31)

C1> α + C2> β + J > γ + b = 0 (2.32)

LPCC problems can be solved by several algorithms ([28] Branch-and-Bound, Bender’s, Big-

M). The big-M approach [28] was found to be one of the most computationally efficient in our

computational experiments. This approach reformulates (2.26)-(2.32) as the following mixed-

integer linear program (MILP).

max b> z (2.33)

d,y,z,α,β,γ,ω

s. t. Constraints (2.27), (2.28), (2.31), (2.32) (2.34)

0 ≤ g1 − Ax − C1 z ≤ M ω1 (2.35)

0 ≤ α ≤ M (1 − ω1 ) (2.36)

0 ≤ g2 − By − C2 z ≤ M ω2 (2.37)

0 ≤ β ≤ M (1 − ω2 ) (2.38)

Here, M is a sufficiently large constant (big-M) and ω1 and ω2 are auxiliary binary variables

that are introduced to enforce the complementarity conditions in (2.29) and (2.30).

2.4.3 Algorithm for the MMC model

The proposed algorithm for the MMC model, which we call AlgMMC , is an iterative one,

in which the master problem is solved to provide an increasing series of lower bound solutions,

and then the subproblem is solved to provide a series of decreasing upper bound solutions using
 20

the solution from the master problem, x, as an input. The input for the master problem, ΩC , is

iteratively enriched by the solutions from the subproblem until the gap between the lower and

upper bounds falls below a tolerance, . Detailed steps of this algorithm are described as follows:

AlgMMC (c, b, A, B, C1 , C2 , g1 , g2 , J, Q1 , q1 , Q2 , q2 )

Step 0 : Initialization. Create ΩC that contains at least one selected scenario. Set LB = −∞,

U B = ∞, and k = 1. Go to Step 1.

Step 1 : Update k ← k + 1. Solve the master problem MC (ΩC ) and let (xk , ξ k ) denote its

optimal solution. Update the lower bound as LB ← c> xk + ξ k and go to Step 2.

Step 2 : Solve the sub-problem SC (xk ) and let (dk , y k , z k ) denote its optimal solution. Update

ΩC ← ΩC ∪ {(dk , y k )}, and U B ← c> xk + b> z k .

Step 3 : If U B − LB > , go to Step 1; otherwise return (xk , dk , y k , z k ) as the optimal solution

to (2.17) and LB as the optimal value.

2.4.4 Algorithm for the MMR model

The MMR model can be solved using a similar algorithmic framework to AlgMMC after the

following simplifying yet equivalent reformulation.

  
> >
min c x + max min b z − G(d, y)
x binary (d,y)∈U z∈Z(x,d,y)
  n o
> > > >
= min c x+ max min b z− min c x̂ + b ẑ
x binary (d,y)∈U z∈Z(x,d,y) (x̂,ẑ)∈G(d,y)
 

  

= min c> x + max min b> z − (c> x̂ + b> ẑ)
x binary
 (x̂,ẑ)∈G(d,y) z∈Z(x,d,y) 

(d,y)∈U

To solve this reformulation of the MMR model, which is structurally similar to the MMC

model (2.17), only slight modifications to the master and sub-problems are required. The

master problem is defined for a different set of input scenarios, ΩR = {(di , y i , x̂i , ẑ i ), ∀i =

1, ..., |ΩR |}, in which the two additional variables, x̂i and ẑ i , represent the optimal investment

and recourse decisions with hindsight of the uncertainty realization (di , y i ). We denote the
 21

master problem as MR (ΩR ), and it is formulated as the following single level mixed integer

linear program.

min c> x + ξ (2.39)

x,ξ,z i

s. t.ξ ≥ b> z i − (c> x̂i + b> ẑ i )∀i = 1, . . . , |ΩR | (2.40)

Ax + C1 z i ≤ g1 ∀i = 1, . . . , |ΩR | (2.41)

By i + C2 z i ≤ g2 ∀i = 1, . . . , |ΩR | (2.42)

Jz i = di ∀i = 1, . . . , |ΩR | (2.43)

x binary. (2.44)

We denote the subproblem as SR (x), and it is formulated as the following bilevel linear

program.
 
> > >
max min b z − (c x̂ + b ẑ) ,
(d,y)∈U ;(x̂,ẑ)∈G(d,y) z∈Z(x,d,y)

which can be solved using the same big-M approach with the following MILP.

max b> z − (c> x̂ + b> ẑ) (2.45)

d,y,z,x̂,ẑ,α,β,γ,ω

s. t. Constraints (2.34)-(2.38) (2.46)

Ax̂ + C1 ẑ ≤ g1 (2.47)

By + C2 ẑ ≤ g2 (2.48)

J ẑ = d (2.49)

x̂ binary. (2.50)

With the new definitions of master and sub-problems, the same algorithm AlgMMC can be

used to solve the MMR model with the following minor modification to Step 2, besides the

apparent need to change the superscript “C” to “R”:

Step 2 : Solve the sub-problem SR (xk ) and let (dk , y k , z k , x̂k , ẑ k ) denote its optimal solution.

Update ΩR ← ΩR ∪ {(dk , y k , x̂k , ẑ k )}, and U B ← c> xk + b> z k − (c> x̂k + b> ẑ k ).
 22

2.5 Case Study

In this section, we present numerical experiments of our model and algorithm on an IEEE

118-bus test system, which consists of 186 transmission lines, 5 wind farms, 5 coal plants, 5 gas

plants and 33 loads. The network data is available in [2].

We consider 10 candidate lines. The operation costs are calculated based on the data of 4

load blocks. We consider a planning horizon of 20 years, with the operation cost extrapolated

from the cost of year 1. Then the operation cost is assumed to increase at the same rate each

year. The characteristics of generation and candidate lines are summarized in Table 2.3 and

Table 2.4 respectively. In our case study, generation capacity data in the system is set to be

able to satisfy all demand levels if there is no network congestion.

Table 2.3 Generation Parameters

Current Fuel Min New Max New
Bus
Type Capacity Cost Capacity Capacity
NO.
(MW) ($/MWh) (MW) (MW)
25 Wind 1,000 0.0 2,000 3,000
31 Wind 500 0.0 150 400
32 Coal 300 43.0 –160 –80
36 Wind 1,100 0.0 800 1,500
49 Wind 1,000 0.0 1,200 2,000
61 Coal 400 28.2 –200 –70
65 Coal 1,500 52.7 –800 –600
66 Coal 300 28.3 –150 –100
76 Gas 200 31.0 80 100
85 Gas 300 63.9 150 200
87 Wind 900 0.0 1,300 2,500
89 Gas 150 49.1 100 150
92 Gas 200 32.0 200 250
99 Coal 300 27.5 –150 –80
113 Gas 200 65.0 200 300
Total 8,350 6,000 8,500

Uncertainty in future generation capacity consists of two parts, expansion and retirement.

For wind and natural gas fired plants, we set the lower and upper bounds for new capacity. For

coal plants, the range of reduced capacity is also provided. We use negative capacity to depict

coal retirement. In addition to bounds on individual plants, we also set the lower bound and

upper bound on total new generation capacity to control the randomness of the uncertainty set.
¯ and the capacity factor are modified based on the real data from
The mean of our demand (d)
 23

Table 2.4 Candidate Line Parameters

Transmission Construction
Susceptance
From To Capacity Cost
(Ω−1 )
(MW) ($m)
25 4 30 390 40.60
25 18 30 390 32.48
25 115 30 390 40.60
32 6 30 390 50.75
36 34 30 390 28.42
36 77 30 390 44.66
70 25 30 390 97.44
86 82 30 390 36.54
87 106 30 390 62.93
87 108 30 390 52.78

WECC [58]. We generate four instances by changing the uncertainty sets. Their definitions

are listed as follows:

U1 = {0.95d¯i,t,m ≤ di,t,m ≤ 1.05d¯i,t,m (2.51)

N,min N N,max
Pi,k,t ≤ Pi,k,t ≤ Pi,k,t (2.52)
X
P N,min ≤ N
Pi,k,t ≤ P N,max } (2.53)
i,k,t
U2 = {0.85d¯i,t,m ≤ di,t,m ≤ 1.15d¯i,t,m (2.54)

Equations(2.52) − (2.53)} (2.55)

U3 = {0.95d¯i,t,m ≤ di,t,m ≤ 1.05d¯i,t,m (2.56)

N,min N,min N
[1.25 − 0.5 sgn(Pi,k,t )]Pi,k,t ≤ Pi,k,t
N,max N,max
≤ [0.75 + 0.5 sgn(Pi,k,t )]Pi,k,t (2.57)
X
0.75P N,min ≤ N
Pi,k,t ≤ 1.15P N,max } (2.58)
i,k,t
U4 = {0.85d¯i,t,m ≤ di,t,m ≤ 1.15d¯i,t,m (2.59)

Equations(2.57) − (2.58)} (2.60)

N,min N,max
where (Pi,k,t , Pi,k,t ) and (P N,min , P N,max ) are listed in the last two columns in Table 2.3,

with (P N,min , P N,max ) in the last row.

The load curtailment cost is set as 2000$/MWh. The interest rate is set as 0.1. The

experiment is implemented on a computer with Intel Core i5 3.30GHz with 4GB memory and
 24

CPLEX 12.5. The computation time of each instance is around 9 hours. The numbers of

iterations for solving each instance are summarized in Table 2.5.

Table 2.5 Number of Iterations for Each Instance

U1 U2 U3 U4
MMC 3 3 2 2
MMR 3 6 2 3

The transmission expansion plans, investment costs and objective values of each criterion

under the four uncertainty sets are summarized in Table 2.6 and 2.7. We then compare the

performances of the MMC solution and the MMR solution under various scenarios in Tables

2.8 and 2.9, where we use Dc , Dr and Dd to denote the optimal MMC solution, the optimal

MMR solution and the optimal deterministic solution. The lower cost and regret between Dc

and Dr are highlighted. The deterministic solution is derived by setting the mean demand as

the load levels and the median of new capacity as the future expansion plans. The scenarios

are generated by our algorithms when solving the MMR problems and MMC problems. Each

scenario corresponds to an optimal solution to a sub-problem at an iteration of our algorithm

and is the worst-case scenario for the first-stage solution obtained at the same iteration. Sce-

narios S1 − S3 and S6 − S10 are generated by solving the MMR problems. Scenarios S4 , S5 , S11

and S12 are generated when solving the MMC problems. Those scenarios typically have very

high costs or regrets, thus are representative of bad scenarios that robust optimization tries to

hedge against. The investment and operational costs of both the MMC and MMR decisions

under the above scenarios are summarized in Table 2.10.

From Table 2.6 and Table 2.7, we can see that as the uncertainty in demand increases,

although the numerical value of the maximum regret and worst-case cost increases, the change

in the transmission expansion plan is not very substantial. It means many of the candidate

lines are necessary regardless of the demand levels with our unchanged depiction of generation

capacity uncertainties. The reason is that those candidate lines connect regions with very high

locational marginal price differences due to the presence of large amount of wind energy. On

the other hand, when uncertainty in generation expansion is increased, although the total cost

also increases, fewer lines are actually built with the MMC criterion. That is because in the
 25

Table 2.6 Expansion plan of the MMC approach

Uncertainty Set
Lines (from bus, to bus)
U1 U2 U3 U4
Candidate Line (25, 4) 1 0 1 1
Candidate Line (25, 18) 1 0 1 1
Candidate Line (25, 115) 0 1 0 0
Candidate Line (32, 6) 0 1 0 0
Candidate Line (36, 34) 1 1 1 1
Candidate Line (36, 77) 1 1 1 1
Candidate Line (70, 25) 0 1 0 0
Candidate Line (86, 82) 1 1 0 0
Candidate Line (87, 106) 1 1 1 1
Candidate Line (87, 108) 1 1 1 1
Investment Cost ($m) 298 414 262 262
Maximum Cost ($m) 1,233 1,960 1,927 2,802

worst-case scenarios, the system contains less renewable energy capacity and the differences

in the locational marginal prices between the otherwise connected regions are not substantial

enough to justify new transmission lines. When the MMR criterion is used, however, the

final expansion plan does not seem to be sensitive to the change of uncertainty in generation

expansion. One possible explanation is that since the MMR criterion does not make decisions

only based on the boundary scenarios, it is less sensitive to the changes in uncertainty sets.

From the more detailed comparisons of results from the MMC and MMR approaches in

Tables 2.8 and 2.9, we can gain more insights about which criterion is more appropriate under

different situations. Both robust optimization solutions outperform the deterministic solution

under most of the scenarios. When the uncertainty set is U1 , the MMC solution outperforms

the MMR solution under most of the listed scenarios. When the uncertainty set is U2 , on the

other hand, the MMR solution has a lower total cost under more listed scenarios. The solutions

of the cases when the uncertainty sets are U3 and U4 yield similar results. According to our

definition at the end of section 2.3, scenarios S2 − S5 in Table 2.8 and scenario S12 in Table 2.9

are regretless scenarios for the MMC decision, while scenarios S1 and S6 −S11 are regretful ones.

Which criterion should be used depends on the decision-maker’s perception on the uncertainty

sets. In typical regretless scenarios for MMC decisions, there usually exists high demand and
 26

Table 2.7 Expansion plan of the MMR approach

Uncertainty Set
Lines (from bus, to bus)
U1 U2 U3 U4
Candidate Line (25, 4) 1 1 1 1
Candidate Line (25, 18) 1 1 1 1
Candidate Line (25, 115) 1 0 1 0
Candidate Line (32, 6) 0 0 0 0
Candidate Line (36, 34) 1 1 1 1
Candidate Line (36, 77) 1 1 1 1
Candidate Line (70, 25) 0 1 0 1
Candidate Line (86, 82) 1 1 1 1
Candidate Line (87, 106) 1 1 1 1
Candidate Line (87, 108) 1 1 1 1
Investment Cost ($m) 339 395 339 395
Maximum Regret ($m) 89 234 110 235

Table 2.8 Comparison of the MMC, MMR and Deterministic Solutions for Uncertainty Set
U1
Cost/Regret ($m) Dc Dr Dd
Scenario S1 1,137 / 117 1,109 / 89 1,085 / 65
Scenario S2 1,039 / 0 1,080 / 41 1,263 / 224
Scenario S3 803 / 0 844 / 41 1,125 / 422
Scenario S4 1,126 / 0 1,167 / 41 1,325 / 199
Scenario S5 1,233 / 27 1,272 / 66 1,367 / 161

low renewable energy penetration. If decision-makers care more about such scenarios or believe

they are more likely, then MMC should be used. Otherwise, choosing MMR might be better.

Both criteria provide good upper bounds for the total costs under scenarios contained in an

uncertainty set. The MMC criterion provides a smaller upper bound with higher average costs

while the costs of MMR decisions are lower on average but have higher variability.

From the above results, it is obvious that both the future generation expansion behavior of

generation companies and demand uncertainty play important roles in transmission expansion

planning. In addition, we can also conclude that both criteria have their merits and can yield

relatively reliable expansion plans that guarantee zero curtailment for uncertainty realizations

contained in the uncertainty sets. However, depending on the characteristics of uncertainty

sets and the preference of decision-makers, they may outperform each other under different
 27

Table 2.9 Comparison of the MMC, MMR and Deterministic Solutions for Uncertainty Set
U2
Cost/Regret ($m) Dc Dr Dd
Scenario S6 1,730 / 170 1,597 / 37 1,861 / 301
Scenario S7 1,375 / 133 1,354 / 112 1,519 / 277
Scenario S8 1,309 / 497 1,046 / 234 836 / 24
Scenario S9 1,398 / 288 1,266 / 156 1,473 / 363
Scenario S10 776 / 246 701 / 171 666 / 136
Scenario S11 1,959 / 201 1,862 / 104 2,082 / 324
Scenario S12 1,960 / 31 1,962 / 33 2,173 / 244

Table 2.10 Investment and Operational costs for scenarios ($m)

Dc Dr Dc Dr
Invest S1 − S5 298 339 Invest S6 − S12 414 395
Operation S1 839 770 Operation S6 1,316 1,202
Operation S2 741 741 Operation S7 961 959
Operation S3 505 505 Operation S8 895 651
Operation S4 828 828 Operation S9 984 871
Operation S5 935 933 Operation S10 362 306
1 Operation S11 1,545 1,467
1 Operation S12 1,546 1,567

situations. Thus, a comparative analysis of the MMC and MMR criteria can shed more light

on better utilization of both approaches.

2.6 Conclusion

In this paper, we propose two robust optimization models for the transmission expansion

planning problem under uncertainty, where we take into consideration both the high-frequency

uncertainty caused by load forecast errors and the low-frequency uncertainty caused by future

generation expansion and retirement. We use two criteria: minimax cost and minimax regret,

and compare their performances. The uncertain parameters are described by a polyhedral

uncertainty set. With this approach, we can derive an expansion plan that is robust under

all scenarios. The resulting models can be formulated as trilevel mixed-integer problems. We

use a branch and cut type mechanism to decompose the problem into a master problem and

a subproblem. The subproblem generates scenarios and returns them to the master problem
 28

to cut off sub-optimal solutions. The bilevel mixed-integer subproblem is reformulated into a

single level mixed-integer-programming problem with the KKT conditions to obtain the global

optimal solution. Our model and algorithm are then tested on an IEEE 118-bus system,

where we compare the results of our MMR and MMC models and analyze their differences.

We conclude that the MMR and MMC criteria may outperform each other depending on the

uncertainty set and decision-maker’s preference. Interesting topics for future research include

developing effective heuristics and parallel computing mechanisms to speed up our algorithms

and implementing our algorithms on high performance machines to test larger systems with

more candidate lines.

 29

CHAPTER 3. ROBUST TRANSMISSION PLANNING UNDER

UNCERTAIN GENERATION INVESTMENT AND RETIREMENT

Published in IEEE Transactions on Power Systems

Bokan Chen and Lizhi Wang

Abstract

Transmission planning is typically faced with a wide range of uncertainty including growth in

demand, renewable energy generation and fuel price. The restructuring of the electric power

industry, the drive for energy independence and the push for a cleaner environment have led to

additional sources of uncertainty in all aspects of power system operations and planning. The

more decentralized decision-making implicates uncertainty in future generation investment and

retirement of existing units. Such uncertainties are among the biggest challenges in transmis-

sion planning in power systems. However, existing approaches in the literature do not fully

address this problem. In this paper, we construct a new robust transmission planning model

representing generation investment and retirement uncertainty more realistically to improve

the quantification and visualization of uncertainty and the impacts of environmental policies.

With this model, we can manage the impacts of uncertainty in the time, the size and the lo-

cation of candidate generation additions as well as the retirement of units. The corresponding

algorithm we develop takes advantage of the structural characteristics of the model so as to

obtain a computationally efficient solution methodology. We present an extensive numerical

study on the WECC 240-bus system [58] and illustrate the effectiveness of our approach.
 30

3.1 Nomenclature

Sets and indices

I Set of nodes, i ∈ I

Q Set of all generators (including existing and candidate generators), q ∈ Q

W Set of candidate generators, q ∈ W

L Set of transmission lines existing before any new investment, (i, j) ∈ L

N Set of candidate transmission lines, (i, j) ∈ N

T Set of periods in the planning horizon, t ∈ T

M Set of load blocks, m ∈ M

K Set of generation technology types, k ∈ K

R Set of renewable generation technology types, k ∈ R

Parameters

Pq,i,k Capacity of generator q of technology k at bus i

cT
i,j,t Cost of building the new transmission line (i, j) at period t

cLi,t Cost of load curtailment at node i at period t

cP
q,i,k,t Cost of power production by generator q of technology k at node i at period t
C
Fk,t,m Average capacity factor of generation technology k at period t for load block

m
max
Fi,j Thermal capacity of transmission line (i, j)

Bi,j Susceptance of transmission line (i, j)

M, M0 Big constants used to linearize the power flow constraints and the products

between two variables

di,t,m Demand at node i at period t for load block m

θmin The lower bound of voltage angles

θmax The upper bound of voltage angles

λ Cash flow interest rate (with inflation)

gq,i,k,0 Availability of generator q of technology type k at bus i at period 0

 31

Variables

xi,j,t Binary variable indicating whether transmission line (i, j) is built at period t

(xi,j,t = 1) or not (xi,j,t = 0)

gq,i,k,t Binary variable indicating whether a generator q of technology type k exists

at bus i in period t (gq,i,k,t = 1) or not (gq,i,k,t = 0)

fi,j,t,m Power flow from node i to node j at period t for load block m

pq,i,k,t,m Power production by generator q of technology k at node i at period t for load

block m

ri,t,m The amount of load shedding at node i at period t for load block m

θi,t,m Voltage angle at node i at period t for load block m

z An aggregate variable representing all dispatch variables, including fij,t,m ,

pq,i,k,t,m , ri,t,m , and θi,t,m

C I (x) Net present value of the investment cost

C O (z) Net present value of the operation cost

3.2 Introduction

Transmission planning is one of the key decision processes in the power system production

pipeline. Electricity transmission networks are responsible for reliably and economically deliv-

ering power from generators to consumers. Thus, a robust and resilient transmission network

is essential to the operation of power systems for decades to come. Good transmission invest-

ments have many benefits, including satisfying increasing demand, promoting social welfare,

improving system reliability and resource adequacy. In addition, in the face of challenges of

climate change, the power generation industry is under increasing pressure to reduce its carbon

footprint, and to introduce more renewable energy into the power system. Transmission enables

a large quantity of renewable energy to be integrated into the power system by bridging the

distances between renewable energy resources and load centers.

Transmission planning is very challenging due to the long planning horizon and various

sources of uncertainty planners need to consider. Traditionally, planners take into consideration

sources of uncertainty such as demand variations and renewable energy intermittency [70,
 32

3], forced outages of transmission lines and generators [40, 18], and transmission capacity

factor [45]. Since the restructuring of the power system, generation expansion planning and

transmission planning are conducted by separate entities, which makes the future generation

investment uncertain to transmission planners. In addition, a large amount of coal plants are

anticipated to retire in response to the EPA clean power plan, to be replaced in part by gas-fired

plants. Such retirement could be unexpected by system operators. All of the above factors

introduce additional sources of uncertainty to the transmission planning process. However,

research on the impacts of uncertain future generation developments on transmission planning

in the decentralized decision making environment is scarce in the literature.

In existing literature, the most popular methods to manage uncertainty in optimization in-

clude stochastic programming and robust optimization. Stochastic programming assumes that

the uncertain parameters follow an estimated probability distribution and generates scenarios

accordingly. It has been applied to many transmission planning problems [3, 40, 18, 45]. Ryan

et al. [61] classify sources of uncertainties in the power system into two categories—high fre-

quency and low frequency. Most operational level uncertainties are of high-frequency, whose

probability distributions are readily available by utilizing historical data. Low frequencies un-

certainties, such as climate change, natural disaster, technology breakthroughs and future gen-

eration investment and retirement, cannot be easily characterized by probability distributions.

Stochastic programming, with its dependence on probability distributions, is more equipped to

deal with high-frequency uncertainties. On the other hand, robust optimization uses paramet-

ric sets to describe uncertainty [74, 13], which can contain infinitely many scenarios without

specific knowledge of the probability distributions. Therefore, robust optimization can also

deal with low-frequency uncertainty. However, the focus in literature has been using robust

optimization to ensure feasibility against any operational level high-frequency uncertainty re-

alizations. In [29], a robust transmission planning model considering demand and renewable

energy uncertainty is proposed. The model is solved by a two-layer Benders’ decomposition

algorithm. Ruiz and Conejo [60] propose a similar model with the uncertainty in equipment

outages also represented. Robust optimization is also used in [50] with the explicit represen-

tation of the system with n − k contingencies for transmission planning. In [5], uncertainty
 33

in construction costs and demand is considered. The cited references are limited to the man-

agement of operational level sources of uncertainty. The impacts of generation investment on

transmission planning are usually explored with other methods.

Generation investment in transmission planning is usually studied with three different

perspectives—deterministic, stochastic and game theoretic. The references cited in the pre-

vious paragraph fall into the deterministic category, which assume generation capacity avail-

able in the system to be known information throughout the transmission planning horizon.

Stochastic methods include scenario analysis and robust optimization. Similar to stochastic

programming, scenario analysis optimizes the average objective among several scenarios with

equal weights that represent various sources of uncertainty. However, such scenarios do not

incorporate any information on probability distribution. Munoz et al. [52] construct 3 scenarios

based on EIA forecasts and some educated assumptions to describe different policy mandates

and fuel costs. Generation investments are determined according to the three scenarios and

are assumed to have the same objective as transmission investments. Buygi et al. [17] study

several cases with different scenarios to explore the effect of low frequency uncertainties. In

[42], multi-period scenario trees are constructed to describe the uncertainty in the size, location

and types of renewable generation. The disadvantages of this approach is that only a limited

number of scenarios are explored. In addition, the constructed scenarios are more of a “big

picture” rather than detailed depictions of possible futures. Unlike scenario analysis, robust

optimization can describe uncertainty with more details. Two robust transmission planning

models with different optimization criteria—minimax cost and minimax regret are proposed

in [20]. Demand and generation capacity uncertainty is considered, where generation capacity

at specific locations is assumed to be in a polyhedral set. Limited by the way uncertainty is

modeled, this work does not explore the impacts of location and time of generation investment

and retirement.

In game theoretic models, generation investment is not a source of uncertainty, but the op-

timal behavior of generation companies (GENCOs), which can be anticipated based on market

condition and the behavior of other players. A trilevel model is proposed in [57]. In the first

level transmission investment decisions are made. In the second level, multiple GENCOs make
 34

generation investments in response to the transmission investment. Operational decisions are

made in the third level. Perfect information competition between GENCOs is assumed. A

similar trilevel model is proposed in [33] and [34], where a Cournot game between GENCOs is

solved. Similar approaches are also used in [59, 51, 27]. Game theoretical models are useful in

that they can provide us some insights into the behavior of different participants in the power

market. However, some strong assumptions about the players are needed. One example is per-

fect information between different players, which usually is not possible. In addition, the high

complexity of the models means that they are difficult to implement for large realistic-sized

systems.

Given the limitations of the methods in the literature, there is a need for practical ap-

proaches to address the challenges faced by transmission planners in the current environment

of the power electric industry. In this paper, we propose a new multi-period robust optimiza-

tion framework for transmission planning with explicit representation of uncertainty in future

generation investment and the retirement of existing units. According to the EIA forecast

[24], the average annual demand growth in electricity is less than 1% in the next 20 years.

Since the uncertainty in demand growth is relatively low, we do not consider demand growth

uncertainty in this paper and focus on the uncertainty in future generation instead. However,

demand uncertainty can be easily incorporated into the uncertainty set for future research. A

customized solution method is developed for the model. The model and solution method is

then tested on the WECC 240-bus system. The impacts of policy and natural gas price change

is also explored in the case study.

Table 3.1 summaries the comparison between our new model and some previous works using

robust optimization. The structure of our new model is more similar to the models in [20], but

a few key aspects are different, including uncertainty modeling, number of decision periods and

solution algorithm.

The contributions of this paper are summarized as follows:

1. A novel method to model generation investment and retirement uncertainty is proposed.

With the more flexible uncertainty modeling, the uncertainty set can describe the uncer-

tainty in not only the capacity of new generation investments but also their locations.
 35

Table 3.1 Comparison with other robust optimization transmission planning models in the
literature
Reference Generation investment Uncertainty Case Study # periods
[29] none load, renewable in- IEEE 118-bus 1
termittency
[60] none load, renewable in- RTS 24-bus 1
termittency, equip-
ment outage
[50] none n − k contingency IEEE 300-bus 1
[5] none investment costs, IEEE 118-bus 1
demand
[20] uncertainty load, generation ca- IEEE 118-bus 1
pacity
Our model uncertainty size, location and WECC 240-bus 4
time of generation
investment and re-
tirement

2. We presented a tri-level robust optimization model that contains multiple decision pe-

riods in the middle and bottom levels, representing when generation and transmission

investments could be made. This feature allows a more detailed schedule of transmission

line construction to be obtained by the model. It also adds the dimension of time into

the description of generation investment and retirement uncertainty.

3. The proposed solution method is more efficient than the existing Karush-Kuhn-Tucker

(KKT) reformulation method [35, 49], which made case studies on larger-sized systems

with multiple decision periods more computationally tractable.

The rest of the paper is organized as follows. Section 3.3 describes a detailed formulation

of the robust transmission planning model. Section 3.4 explains the solution method used in

this paper. Section 3.5 presents the results of our numeric experiments. Finally, Section 5

concludes the paper.

3.3 Model Formulation

Transmission planning is a two-stage decision making problem. The two stages refer to

the decision making processes before and after uncertainty realization. The first-stage deci-

sions usually cannot be changed after uncertainty is observed. The second-stage decisions are
 36

Transmission plan proposition

6
plan max cost
?
Worst scenario identification
6
scenario cost
?
Optimal power flow

Figure 3.1 Summary of the robust optimization modeling framework

heavily influenced by the particular realization of uncertainty as well as the first-stage deci-

sions. In our robust optimization model, new transmission lines are selected for construction

with limited information on uncertainty during the first decision making stage. In the second

stage, operation decisions are made after uncertainty realization. A summary of the two-stage

modeling framework is shown in Figure 3.1.

In addition, multiple decision periods are captured in our model. Transmission planning

usually spans a long planning horizon. Due to budget and other physical constraints, construc-

tion of all the needed transmission lines at the same time might be difficult or unrealistic. It

stands to reason that the planned transmission lines may be available at different points along

the planning horizon. In this model, we assume the planned transmission lines are available

at the beginning of each decision period. The first-stage decision consists of transmission line

construction plans for the multiple decision periods in the planning horizon.

The goal of robust optimization is to identify a schedule to construct new transmission

lines that can achieve the lowest total cost, including construction cost and estimated operation

cost. In this section, we present a two-stage multi-period robust transmission planning problem

formulated as a trilevel optimization model.

 
I O
min C (x) + max min C (z) (3.1)
x∈X g∈G z∈Z(x,g)

where
 37

 X
X = x| xi,j,t ≤ 1, ∀(i, j) ∈ N ,
t

xi,j,t Binary, ∀(i, j) ∈ N , t ∈ T (3.2)
X 1
C I (x) = cT xi,j,t (3.3)
(1 + λ)t i,j,t
t,i,j
X 
O
X 1 P L
C (z) = c pq,i,k,t,m + ci,t ri,t,m (3.4)
(1 + λ)t q q,i,k,t
i,t,m


Z(x, g) = z : fi,j,t,m = Bi,j (θi,t,m − θj,t,m ),

∀t, m, (i, j) ∈ L [µ(1) ] (3.5)

t
X
|fi,j,t,m −Bi,j (θi,t,m − θj,t,m )| ≤ (1 − xi,j,k )M,
k=1

∀t, m, (i, j) ∈ N [µ(1) ] (3.6)

t
X t
X
max max
− Fi,j ( xi,j,k ) ≤ fi,j,t,m ≤ Fi,j ( xi,j,k ),
k=1 k=1

∀t, m, (i, j) ∈ N [µ(2) , µ(3) ] (3.7)

max max
− Fi,j ≤ fi,j,t,m ≤ Fi,j ,

∀t, m, (i, j) ∈ L [µ(2) , µ(3) ] (3.8)

C
0 ≤ pq,i,k,t,m ≤ Fk,t,m Pq,i,k gq,i,k,t ,

∀q, i, k, t, m [µ(4) ] (3.9)

θmin ≤ θi,t,m ≤ θmax , ∀i, t, m [µ(5) , µ(6) ] (3.10)

X X X
pq,i,k,t,m + fj,i,t,m − fi,j,t,m = di,t,m
q j j

− ri,t,m , ∀i, t, m [µ(7) ] (3.11)

and
 
G= g|Hg ≤ h, gq,i,k,t Binary, ∀q, i, k, t . (3.12)

From the model formulation (3.1), we can see the two-stage structure of the robust transmis-

sion planning model. Transmission construction schedule x is determined in the first stage. In
 38

the second stage, after generation investment and retirement g is observed, operational decision

z is made. The objective is to minimize the investment cost and the projected operational cost

under the most costly scenario. Investment and operation costs are represented by Equations

(3.3) and (3.4) respectively. The operation cost includes generation cost and load curtailment

cost (the two terms in the parentheses in Equation (3.4)). Constraint (3.2) means that a

transmission line only needs to be built once in the planning horizon. The uncertainty set G

is defined in Equation (3.12), which can be modified based on planners’ belief on uncertain

parameters. Thus the definition of the uncertainty sets is case-specific. The symbols H and h

represent, respectively, the coefficient matrix of the constraints and the right-hand-side vector.

The detailed formulation of the uncertainty set used in our case study can be found in Section

3.5.

Constraints (3.5)-(3.11) define Z(x, g), the set of feasible power flow solutions after new

transmission and generation is fixed. Equation (3.5) defines the power flow on existing trans-

mission lines, while Constraint (3.6) defines the power flow on candidate lines. It is equivalent

to the constraint fi,j,t,m = ( tk=1 xi,j,k )Bi,j (θi,t,m − θj,t,m ), which is nonlinear and can com-
P

plicate the computation. By using the auxiliary parameter M , it is linearized. Constraints

(3.7)-(3.8) impose power flow limits on candidate and existing transmission lines. Constraint

(3.9) requires power generation to be within the capacity limit of generators after capacity

factor is considered. Constraint (3.10) limits the range of difference between phase angles at

each node and the phase angle of the reference node. Constraint (3.11) requires the net influx

at a node to be equal to the net outflow. The symbols in the brackets represent the dual

variables corresponding to the constraints. The loads here are modeled by load blocks. The

dual variables for Constraints (3.5) and (3.6) only differ in the indices, so do the dual variables

for (3.7) and (3.8). Constraints (3.7), (3.8) and (3.10) each have two dual variables, one for

each side of the inequalities.

3.4 Solution Methods

From the decision making perspective, the proposed model is a two-stage multi-period

model, where transmission investment decisions are made in the first stage and the operational
 39

decisions are made in the second stage. From the solution technique perspective, the model

is also a trilevel robust optimization model, in which the two-stage decisions are made at the

top and bottom levels and the middle level identifies the worst-case scenario of uncertainty

realization. Trilevel optimization problems can be difficult to solve. In this section, we present

a column and constraint generation (CCG) algorithm [74] to decompose the trilevel model

into a master problem that makes the first level decisions and a subproblem that solves the

middle and bottom levels. As such the subproblem is a bilevel optimization problem. The

CCG method is similar to Benders’ decomposition [12] in that both generate cutting planes

iteratively. The difference is that Benders’ decomposition relies on the dual formulation of the

linear program subproblems to generate cuts, whereas the CCG method introduces both new

cuts and new variables by solving the bilevel subproblem.

One of the most popular methods for solving bilevel optimization problems is reformulating

the lower level as complementarity constraints using its KKT conditions and then linearizing it

using the big-M method [20]. However, such method can be inefficient due to the large number

of auxiliary binary variables introduced in the linearization step. In our solution method, we

first dualize the subproblem by introducing bilinear terms, and then linearize them by exploiting

the discrete structure of the uncertainty set. As such, the bilevel subproblem is reformulated

as a single level mixed integer linear program. Since less binary variables are introduced, the

solution method is more efficient, which will be demonstrated in Section 3.5.

3.4.1 Master problem and subproblem

The master problem is a relaxation of the robust transmission planning problem and pro-

vides a first-stage solution. The original trilevel problem requires the most costly scenario in

the entire uncertainty set G. In the master problem, the search for the worst-case scenario is

restricted in a smaller set of scenarios S ⊆ G. We denote the master problem as M(S), and it

is formulated as the following mixed-integer linear program:

 40

minx,z s C I (x) + ζ (3.13)

s. t. x ∈ X (3.14)

ζ ≥ C O (z s ), ∀s ∈ S (3.15)

z s ∈ Z(x, g s ), ∀s ∈ S. (3.16)

When S = G, the master problem M(S) is equivalent to the trilevel formulation (3.1). How-

ever, due to the enormous scenario space in G, it is unrealistic to search the entire uncertainty

set. Instead, we search the subset S. As such, the master problem provides a lower bound to

the actual optimal objective value. New scenarios are iteratively added to the set S and the

lower bound increases at each iteration.

The subproblem identifies the worst-case scenario given a first-stage decision x. Since

the performance of such solution can only be worse than the optimal first-stage solution, the

subproblem provides an upper bound to the optimal objective value for the original trilevel

problem. We denote the subproblem as R(x) and it can be formulated as the following bilevel

linear program:

max min C O (z) (3.17)

g∈G z∈Z(x,g)

3.4.2 Solving the subproblem

The subproblem defined in (3.17) is a bilevel optimization problem. In order to solve

it, we write out the dual of the operation cost minimization problem min C O (z), and the
z∈Z(x,g)
subproblem R(x) can be replaced by the following formulation.
 41

X X (2) (3)
max
max − Fi,j (µi,j,t,m + µi,j,t,m )
µ,g
t,m (i,j)∈L∪N

(4)
X
C
− Fk,t,m Pq,i,k gq,i,k,t µq,i,k,t,m
q,i,k

(5) (6)
X
− (θmax µi,t,m − θmin µi,t,m )
i

(7)
X
− di,t,m µi,t,m (3.18)
i
(4) (7)
s. t. µq,i,k,t,m + µi,t,m ≥ −(1 + λ)−t cP
q,i,k,t ,

∀q, i, k, t, m [p] (3.19)

(7)
µi,t,m ≥ −(1 + λ)−t cLi,t , ∀i, t, m [r] (3.20)
(1) (2) (3) (7) (7)
µi,j,t,m + µi,j,t,m − µi,j,t,m − µi,t,m + µj,t,m = 0,

∀(i, j) ∈ L ∪ N , t, m [f ] (3.21)
(5) (6) (1)
X
µi,t,m − µi,t,m − Bi,j µi,j,t,m
(i,j)∈L
t
(1) (1)
X X X
+ Bj,i µj,i,t,m − Bi,j ( xi,j,k )µi,j,t,m
(j,i)∈L (i,j)∈N k=1
t
(1)
X X
+ Bj,i ( xj,i,k )µj,i,t,m = 0,
(j,i)∈N k=1

∀i, t, m [θ] (3.22)

µ(k) ≥ 0, k = 2, 3, 4, 5, 6. (3.23)

In this formulation, µ(2) to µ(6) are the dual variables of Constraints (3.5)-(3.11). Note that

because the first-stage decision variable xi,j,t is a parameter in the subproblem, we can replace

(3.6) with its equivalent form fi,j,t,m = ( tk=1 xi,j,k )Bi,j (θi,t,m − θj,t,m ) to simply the formula-
P

tion. In addition, with the previous replacement, we can combine Constraints (3.7) and (3.8)

by removing the term tk=1 xi,j,k . This explains why we can use the same dual variables µ(2)
P

and µ(3) for Constraints (3.7) and (3.8).

C (4)
In the objective function (3.18), the terms Fk,t,m Pq,i,k gq,i,k,t µq,i,k,t,m are bilinear. Since

they are multiplications between a continuous variable and a binary variable, we linearize
 42

(4)
them by replacing gq,i,k,t µq,i,k,t,m with variable ηq,i,k,t,m and add two sets of constraints to the

subproblem. The constraints are as follows.

(4)
ηq,i,k,t,m ≥ µq,i,k,t,m − M 0 (1 − gq,i,k,t ), ∀q, i, k, t, m (3.24)

ηq,i,k,t,m ≥ 0, ∀q, i, k, t, m (3.25)

Then the subproblem can be easily solved as a mixed-integer linear program.

3.4.3 Algorithm

The algorithm uses the column and constraint generation procedure to iteratively generate

new scenarios to be included in the subset S. The master problem is solved in each iter-

ation to provide an increasing series of lower bounds and a first-stage solution. Given the

aforementioned first-stage solution, the subproblem is solved to provide an upper bound and

a worst-case scenario to be included in set S. The algorithm terminates after the optimality

gap (difference between the upper and lower bounds) falls below a user defined threshold .

Since the uncertainty set is a finite discrete one, the algorithm converges in a finite number of

iterations. More detailed theoretic proof of the correctness and convergence of the algorithm

can be found in [28]. A flow chart of the algorithm is summarized in Figure 3.2. Detailed steps

of this algorithm are described as follows:

Step 0 : Initialize by creating S that contains at least one selected scenario. Set lower bound

LB = −∞, U B = ∞, and s = 1. Go to Step 1.

Step 1 : Update s ← s+1. Solve the master problem M(S) and let (xs , ζ s ) denote its optimal

solution. Update the lower bound as LB ← ζ s . Go to Step 2.

Step 2 : Solve the subproblem R(x) and let (g s , z s ) denote the optimal solution. Update

S ← S ∪ {g s }, and U B ← C O (z s ).

Step 3 : If U B − LB > , go to Step 1; otherwise return xs as the optimal investment plan

and C I (xs ) + ζ s as the optimal value.

 43

Initialization

Master Problem

Subproblem yes

U B − LB > ?

no
Stop

Figure 3.2 Flow chart of the decomposition algorithm

3.5 Case Study

In this section, we present a case study of our model and algorithm on the modified WECC

240-bus test system [58]. The topology of the test system is shown in Figure 3.3, which is from

[52]. The system has 240 buses, 448 transmission lines, 139 load centers and 124 generation

units, including coal, gas-fired, nuclear, hydro, geothermal, biomass, wind, solar and other

renewable generators. The definition for generator type “renewables” can be found in [58].

We consider 18 candidate lines. The planning horizon of 20 years is divided into four five-

year periods. As such, the solution space consists of 518 possible solutions. We select the load

of two representative hours from the typical week market data to be the baseline of our demand

data. An annual growth rate of 1% is assumed.

We assume all of the 17 coal power plants will retire over the study period. Uncertainty

lies in the time of their retirement. We also consider 53 candidate gas-fired units, wind and

solar farms, and geothermal generators at various locations in the system, a subset of which

will receive investment and become available during the planning horizon. When and which

generators will be built is uncertain. The total number of generators with uncertainty adds up

to 70. As such, the uncertainty set consists of 417 × 553 possible scenarios. In the uncertainty

set, we specify several requirements for the uncertain generators. Firstly, 50% of the new

capacity needs to be included at the end of the planning horizon. Secondly, the total additional

capacity at a decision point cannot be lower than 50% of the total additional capacity at the

next decision point. Finally, we require that renewable energy generation account for at least
 44

Figure 3.3 WECC 240-bus test system topology from [52]

 45

a certain percentage of the total additional capacity. These restrictions on the uncertainty set

can be adjusted based on available information as well as planners’ belief on future generation

profiles. The load curtailment cost is set as $2000/MWh. The discount rate is set to be

λ = 10%.

We visualize the test system as well as the transmission planning problem in Figure 3.4.

The 240 buses are arranged in a circle, with bus #1 starting at the three o’clock position and

going counter-clockwise; bus #240 also comes back to the three o’clock position. Although the

spatial information of the buses is distorted by this arrangement, it allows for more informative

visualization of the big picture for the transmission planning problem. The 448 transmission

lines are represented by black line segments connecting the corresponding bus pairs. The

relative lengths of the line segments in the figure do not necessarily represent the actual lengths

of the transmission lines. The 18 candidate lines are also plotted in green. These lines do not

exist in the test system yet, and it is up to the planner to decide which (if any) lines should

be built and when. The squares surrounding the inner circle represent the 139 load centers,

with the areas of the squares proportional to the magnitudes of the loads at the corresponding

nodes. The colorful dots surrounding the squares represent the existing 124 generators, with

the areas of the dots proportional to the capacities of the generators at the corresponding

nodes. The colorful stars represent the potential new generators that may be added to the

corresponding nodes. The color map on the left sub-figure is for both existing and potential

new generators. Multiple generators at the same nodes are plotted at different orbits to avoid

overlap. The uncertainty set used in this case study is defined in Equation (3.26)-(3.31).

Equation (3.26) means that once a coal power plant is retired, it remains retired throughout

the rest of the planning horizon. Equation (3.27) means all coal power plants should be retired

by the end of the planning horizon. Equation (3.28) means once a new generator is brought

into the system, it remains there. Equation (3.29) means a certain percentage, denoted as w1 ,

of the total new generation capacity should be available by the end of the planning horizon.

Equation (3.30) requires that the capacity change in the system be gradual, where w2 is the

highest allowed ratio between the variable capacity in the system in one period and that in

the previous period. Equation (3.31) means certain percentage of the investment should be
 46

renewable resources, where w3 is the percentage. In this case study, we set w1 = 0.7, w2 = 0.6,

w3 = 0.2 or 0.4.

G = g|gq,i,k,t ≥ gq,i,k,t+1 , k ∈ {Coal}∀q, i, t (3.26)

gq,i,k,tmax = 0, k ∈ {Coal}, ∀q, i (3.27)

gq,i,k,t−1 ≤ gq,i,k,t , k ∈ K \ {Coal}, ∀q, i, t (3.28)

X X
gq,i,k,tmax Pq,i,k ≥ w1 Pq,i,k (3.29)
q,i,k q∈W,i,k∈K
X X
gq,i,k,t Pq,i,k ≥ w2 gq,i,k,t−1 Pq,i,k , ∀t (3.30)
q,i,k q,i,k
X X 
gq,i,k,t Pq,i,k ≥ w3 gq,i,k,t Pq,i,k , ∀t . (3.31)
q,i,k∈R q,i,k∈K\{Coal}

We define four futures with respect to different scenarios of natural gas prices and policy

requirement on renewables in order to analyze the sensitivity of the transmission planning

solutions from our model.

ˆ Future 1: High natural gas prices (double the current price) with a mandate of 20%

additional renewables.

ˆ Future 2: Low natural gas prices (at the current level throughout the next twenty years)

with a mandate of 20% additional renewables.

ˆ Future 3: High natural gas prices with a mandate of 40% additional renewables.

ˆ Future 4: Low natural gas prices with a mandate of 40% additional renewables.

The methodology was implemented on a desktop computer with Intel Core i7 3.4GHz CPU,

8 GB memory and CPLEX 12.5. We set the optimality tolerance threshold  as 10−3 , which

means that the solution we found was at most $0.001 more costly than the global optimal

solution. The algorithm usually terminated within 3 to 7 iterations, with each iteration taking

approximately two hours. In comparison, we also tried solving the case study using the KKT

based algorithm from [49], but it was not able to terminate within days.
 47

biomass
natural gas
geothermal
hydro
nuclear
renewables
solar

wind
coal

Figure 3.4 Visualization of the transmission planning problem

 48

We obtained six transmission plans. Plans 1 to 4 are the optimal transmission plans under

Futures 1 to 4, respectively using the proposed robust optimization model. Plans 5 and 6

are two other transmission investment plans that were constructed to represent non-optimal

planning decisions. Plan 5 represents the under-investment case when the planner fails to add

sufficient new transmission capacity, whereas Plan 6 represents the over-investment case, which

is to invest in 15 out of the 18 candidate lines at once in the first period. These six transmission

plans are visualized in Figure 3.5. Each column represents a plan, and each row represents a

five-year period in the planning horizon. The green lines are to be added to the system in the

period according to the transmission plan. Although the added new lines are cumulative, only

new additions are shown in green. It can be seen that the first four transmission plans are

similar and all build new lines gradually in the first three periods. Plan 5 has only a small

number of lines, and plan 6 adds 15 out of 18 candidate transmission lines at once in the first

period.

A by-product of the algorithm from Section 3.4 is a set of scenarios. We have collected

twelve of these scenarios under the four futures and use them to illustrate the robustness

performances of the six plans under these scenarios. The scenarios are shown in Figure 3.6.

The four orbits of the stars indicate the periods in which the new generators are added to the

generation portfolio. The higher the orbit, the later in time. The timing of the retirement of

the coal power plants is not visualized in those scenarios. The label ‘W1’ is for the worst-case

scenario (with highest cost) for Plan 1, ‘B2’ for the best-case scenario (with lowest cost) for

Plan 2, and so on. It is noticeable in the figure that in the worst-case scenarios new generators

tend to emerge later and in smaller amounts, most of which are expensive renewables; whereas

in the best-case scenarios more of them start to serve the demand earlier in time, and they

include more natural gas generators, which are less costly.

Finally, the robustness of the six transmission plans under the twelve scenarios was evaluated

with respect to the investment cost, the operations cost, and load curtailment (as a percentage

of total load) over the planning horizon in Figure 3.7. The horizontal axis of each sub-graph

is the 20-year planning horizon, containing four discrete 5-year periods. The blue horizontal

lines in the first row represent the net present values of the total investments of the plans. The
 49

T1

T2

T3

T4

plan 1 plan 2 plan 3 plan 4 plan 5 plan 6

Figure 3.5 Six transmission plans.

 50

Figure 3.6 Twelve scenarios. The label ‘W1’ is for the worst-case scenario for Plan 1, ‘B2’ for
the best-case scenario for Plan 2, and so on.
 51

colored regions in the bottom two rows in Figure 3.7 are outlined by the upper and lower bounds

performances under the twelve scenarios. The white curves inside the colored region are the

intermediate scenarios. The figure suggests that the first four plans have similar investment

costs, with Plan 4 being the least costly. Their performances in operational costs and load

curtailments are also comparable. Plan 5 has a negligible amount of investment, which is a

$1.4B savings from Plan 4. However, as a consequence of the lack of enough investment, its

operational cost is almost $20B more than Plan 4 over the planning horizon, and its load

curtailment in the worst case is also much severer than Plan 4. On the other extreme, Plan

6 makes about 20% more investment than Plan 4, but its performance in operations cost and

load curtailment are similar, if not even worse than the cheaper counterpart. These results

demonstrate the need for making enough and smart investment in transmission planning in

order to reduce the long-term operations cost as well as enhancing the reliability of the grid.

Another observation we would like to make is the subtle differences among the first four

plans. According to their performances in Figure 7, Plans 1 and 2 appear to be more similar

and Plans 3 and 4 look more alike. This observation suggests that the transmission plans

might be more sensitive to the renewable mandates (that set Plans 1 and 2 apart from 3 and 4)

than they are to the natural gas prices. This observation is consistent with the intuition that

transmission investment decisions should be driven by the availability of generation resources

(or the lack of them) more than their fuel prices.

3.6 Conclusion

In this paper, we propose a multi-period robust transmission planning model focusing on

the management of the impact of generation investment and retirement uncertainty. By us-

ing a multi-period model, a more detailed construction schedule for transmission lines can be

derived compared to the static single period models. We propose a flexible modeling frame-

work to represent generation investment and retirement uncertainty in more detail. Under

the new modeling framework, uncertainty concerning the time, location and size of generation

investment can be described in the uncertainty set. To solve the model, we use a column and

constraint generation procedure to decompose the two-stage problem into a master problem
 52

plan 1 plan 2 plan 3 plan 4 plan 5 plan 6

Figure 3.7 Performance of six plans. The horizontal axis of each sub-graph is the 20-year
planning horizon, containing four discrete 5-year periods.

and a subproblem. By using the structural characteristic of our model, the bilinear subproblem

can be reformulated as a mixed-integer linear problem without utilizing the KKT conditions,

which also makes a larger test case more tractable. We apply our model and solution method

to the WECC 240-bus test system. Then we compare the performance of the robust transmis-

sion plans derived by our model to the performance of plans derived from solving deterministic

models or other heuristic rules under four different scenarios, which are obtained by adjusting

natural gas price and renewable energy policy mandate. The results show that our transmis-

sion plans not only out-perform plans obtained from other means but also remain robust under

different scenarios.

Several directions of possible extensions to this work are worth investigating for future
 53

research. One way to extend this work is to combine robust optimization and stochastic pro-

gramming by modeling various sources of uncertainty with different methods. For example, we

can generate scenarios for operational level uncertainty such as load and renewable energy gen-

eration and use uncertainty sets for other types of uncertainty. To solve such a problem, and to

improve the computational efficiency of current models, new algorithms need to be developed.

Another way to extend the work done to date is to incorporate additional details into the model,

including explicit representations of various existing and possible policies, changes in demand

caused by demographic development, climate change and recovery from natural disasters. As

a caveat, our model was based on DC power flow, which is known to be less accurate than AC

power flow model, although the latter would introduce non-linearity and non-convexity and

make the model much harder to solve. It would be a meaningful future research to examine the

extent to which the simpler DC power flow could compromise the quality of the transmission

planning decisions. We also note that the proposed model is able to represent different types

of transmission lines that differ in cost and thermal capacity; but other types of technological

variations are not reflected.

 54

CHAPTER 4. A STOCHASTIC PROGRAMMING FRAMEWORK FOR

OPTIMAL DISTRIBUTION NETWORK HARDENING

Bokan Chen, Jianhui Wang, Lizhi Wang

Abstract

The power industry is paying more attention to improving system resilience by hardening

the infrastructure in face of increasingly frequent natural disasters. In this paper, we propose

a stochastic programming framework to optimally allocate limited budget to the hardening of

the power distribution networks. The stochastic programming model is formulated as a two-

stage problem, where first-stage decisions are investment strategy selections, and second-stage

decisions include optimal power flow and load shedding. Different hardening methods and their

impact on the failure probability of the distribution lines are considered. To consider all possible

scenarios would make the stochastic program intractable. Sample average approximation is

used to reduce the size of the problem. A reasonably large number of scenarios are simulated

to improve the accuracy of the approximation. In order to solve the resulting large scale

mixed-integer program, the L-shaped method (Benders Decomposition) is used to decompose

the model into a master problem and a group of subproblems. A heuristic algorithm is also

developed to get a reasonably good solution in a short time. The model and algorithms are

tested on an IEEE 33-bus distribution system to demonstrate their effectiveness. The results

are then validated through the single replication procedure.

 55

4.1 Nomenclature

Sets and indices

I Set of nodes, i ∈ I

L Set of distribution lines, (i, j) ∈ L

T Set of time periods in consideration, t ∈ T

K Set of hardening techniques, k ∈ K

S Set of scenarios, s ∈ S

Parameters

cH
i,j,k Cost of hardening line (i, j) with technique k

B Total hardening budget

xi,j Reactance of power line (i, j)

ri,j Resistance of power line (i, j)

dP
i,t Real power demand at node i at time t

dQ
i,t Reactive power demand at node i at time t

Pimax Real capacity of distributed generator i

Qmax
i Reactive capacity of distributed generator i

z̃i,j,t,k Binary random variables indicating whether power line (i, j) hardened with
s
strategy k is damaged at time t. We use zi,j,t,k to represent a sample of z̃i,j,t,k .

z̃i,j,t,k = 1 means the power line is damaged; z̃i,j,t,k = 0 otherwise

M Large constant used to limit the power flow based on the condition of the

distribution line

V min The lower bound on voltage levels

V max The upper bound on voltage levels

V0 Reference voltage magnitude

 56

Decision Variables
P
fi,j,t Real power flow on line (i, j) at time t
Q
fi,j,t Reactive power flow on line (i, j) at time t

pi,t Real power generated by generator at node i at time t

qi,t Reactive power generated at node i at time t

vi,t Voltage magnitude at node i at time t

yi,j,t Binary variable indicating whether power line (i, j) is operational at time t.

yi,j,t = 1 means the line is operational; yi,j,t = 0 otherwise

hi,t The percentage of load not served at node i at time t

ui,j,k Binary variable indicating whether line (i, j) is hardened with strategy k.

ui,j,k = 1 means the line is hardened; ui,j,k = 0 otherwise

4.2 Introduction

Electricity has become an integral part of people’s day-to-day life. Transmission and distri-

bution (T&D) networks are essential infrastructures for the power grid. In the United States,

there are more than 200,000 miles of high voltage transmission lines and numerous distribution

lines that span the entire country. The vast network of T&D lines enables the transportation

and distribution of cheap electricity to countless consumers. In addition to the existing infras-

tructure, the networks are expanding to accommodate increasing demands and integration of

renewable energy resources. However, such vast networks are vulnerable to disruptions caused

by natural disasters including seismic activities and extreme climate events. The President’s

Council of Economic Advisers and the U.S. Department of Energy reported an estimated 679

widespread power outages related to severe weather conditions between 2003 and 2012, leading

to an annual average of 18 to 33 billion dollars of cost to the economy [54]. Due to climate

change, the frequency and intensity of natural disasters such as floods, droughts, hurricanes

and other extreme weather activities are expected to rise [63]. Improving power system se-

curity and resilience not only would reduce the severity of operation disruptions and save the

economy billions of dollars, but also could be instrumental to the recovery of communities from
 57

catastrophes.

The resilience concept includes two levels of capabilities: the inherent ability of a system to

cope with disruptions and the ability to recover from such disruptions [48]. Both abilities of the

power system could be improved to reduce the impact of natural disasters through investment

in power grid modernization and resilience. One of the most effective ways to improve grid’s

ability to resist disruptions is by the hardening of existing components [54]. Investment in

hardening can also help to reduce response time and expedite recovery after the occurrence of

disruptions. Methods of hardening include upgrading power poles and structures with stronger

material, vegetation management and undergrounding utility lines. Hardening of the entire

grid simultaneously could be prohibitively expensive and therefore is unrealistic. Thus, such

investments are usually constrained by a limited budget. It is crucial to allocate the limited

resources to the most vulnerable components or to the power lines with the largest impact to

system performance if affected.

Research on power system resilience planning is relatively limited in scope, with a lot of work

focusing on transmission networks. Several previous studies use the Stackelberg game theoretic

framework that models terrorists and natural disasters as adversaries aim to maximize damage

on the power system. A bilevel programming approach was proposed in [6] to identify the

most critical components in the power system. Other studies take the problem one step further

by adding another level to the model that optimizes the preventive measures taken by the

network planner before disruption occurs. Yao et al. [69] first proposed the defender-attacker-

defender trilevel model for the resource allocation problem in power system defense. Alguacil

et al. [4] proposed a two-stage trilevel model for power grid defense planning, where planners

make moves before and after an attack occurs to minimize the maximum damage attackers can

inflict. A similar model was proposed in [72] with improved solution methods using the column

and constraint generation algorithm. These models can provide valuable insights on how to

identify vulnerable components in a system and what to do with such information, but they

do have some shortcomings. One drawback of the game theoretic models is that some strong

assumptions (perfect information between players, for example) need to be made, which might

not be true in some cases. In addition, due to the complexity of trilevel models, details giving
 58

nuances to the problem might be lost.

Designing adaptive network topology that can isolate the faulty components or reconfigure

the larger network into smaller functioning networks is another popular direction in the research

of improving power system resilience. Baran et al. [8] developed two power flow formulations

for network reconfiguration to reduce loss and balance load. Chen et al. [21] proposed to break

a large network into several small microgrids in the aftermath of disasters to pick up as much

critical loads as possible. In [53, 56], intentional islanding was proposed as a method to minimize

blackouts by creating a collection of operational or nonoperational islands to prevent blackouts

in a larger scale. These works focus on corrective actions after damages to the network have

been caused by natural disasters. The work in this paper, on the other hand, represents efforts

aim to prevent such damages from occurring.

Due to the stochastic nature of natural disasters, uncertainty is an important facet of

the resilience planning problem. As an optimization framework for decision making under

uncertainty, robust optimization has been used to improve the resilience of power systems.

The aforementioned research in the trilevel defender-attacker-defender models [4, 72], etc. can

be seen as special cases of robust optimization. In addition, Yuan et al. [71] used robust

optimization to solve the resilient distribution network planning problem. However, robust

optimization uses polyhedral sets to describe uncertainty and optimize the objective function

under the worst-case scenario. Due to the structure limitations of the model, it is very difficult

to capture the dependency between the first-stage decisions and the uncertain parameters in a

two-stage model. For example, in [71], it is assumed that hardened distribution lines would be

invulnerable to the effect of natural disasters, which might not be true in many cases.

Stochastic programming, on the other hand, uses scenarios generated according to the

probability distributions of uncertain parameters to characterize uncertainty, which is very

suitable for the description of the impact of natural disasters on the power grid. Although

its application in power system resilience is relatively limited, stochastic programming has

been successfully applied to freight transportation network resilience planning [48, 22], which

shares several similarities with power transmission and distribution networks. For example,

both types of networks have high levels of interconnectedness that make them vulnerable to
 59

disruptions caused by natural disasters. Yamangil et al. [68] proposed a two-stage stochastic

programming model for designing resilient electrical distribution grids by altering the topology

of the networks, in which power flow is approximated by a multi-commodity network flow

model. Damage scenarios from natural disasters are modeled as stochastic events to estimate

the performance of the network in the aftermath of a disaster.

In this paper, we propose to use a stochastic programming approach to study the resilience

planning problem. The effect of disasters on distribution networks can be easily described

by probability distributions (Bernoulli distribution is used in this paper). The treatment of

uncertainty in stochastic programming is more flexible than robust optimization. As a result,

we are able to capture the fact that hardened lines are still susceptible to damages, albeit with

lower probabilities. The resilience planning problem is modeled as a two-stage stochastic mixed-

integer programming problem. The objective of this model is to identify the best strategy

to allocate limited budget to critical links in the power grid and improve the resilience of

distribution networks. In this problem, the first-stage variables are planning decisions on

the selection of hardening strategies. The second-stage decisions include the optimal power

flow, generation and load shedding. Sample average approximation and the L-shaped method

(Benders decomposition) are used to reduce the size of the optimization model and speed up

the computation. A heuristic algorithm based on the Knapsack problem is also developed to

produce a reasonably good solution quickly. The model and algorithms are then applied to

an IEEE 33-bus distribution system to test the effectiveness of the proposed methods. The

contributions of this paper can be summarized as follows:

1) The proposed model can take the recovery time of the damaged components into considera-

tion. In addition, the effects of hardening different components in the gird are considered.

2) The proposed approach models the dependencies between uncertainty parameters and de-

cision variables.

3) A new heuristic method is developed to solve optimization problems with Knapsack con-

straints, which can identify a good solution quickly.

 60

The rest of the paper is organized as follows. Section 4.3 presents the two-stage stochas-

tic resilience planning model formulation in detail. Section 4.4 describes the sample average

approximation and Benders decomposition methods used in this paper for computation. The

single replication procedure and the Knapsack based heuristic is also presented. Section 4.5

outlines the data preparation (scenario generation) procedure for the case study and presents

the numerical results. Finally, Section 4.6 provides a conclusion and some discussions on future

research.

4.3 Model Formulation

The optimal resilience planning problem is formulated as a two-stage stochastic program-

ming model. In the two-stage model, network hardening decisions need to be made before ex-

treme weather events happens, without the information on the network damage status caused

by such events. Then after the damage in the network is observed, the optimal operational

decisions can be made.

One difficulty in the modeling process is taking the impact of investment decisions on the

uncertain parameters into consideration. In most stochastic programs, uncertain parameters do

not depend on decisions, which is not the case in this paper. Stochastic programs with decision

dependent uncertainty is explored by Jonsbråten et al. [37] for a specific class of problems. In

this paper, the decision-dependency of random variables is achieved by using binary variables

to select the appropriate set of uncertainty realizations, which simplifies the solution process.

The two-stage model can be formulated as follows:

min Ez̃ [Q(u, z̃)] (4.1)

u∈U

where
 X
U= u| cH
i,j,k ui,j,k ≤ B (4.2)
i,j,k
X
ui,j,k = 1, ∀(i, j) ∈ L (4.3)
k

ui,j,k Binary, ∀(i, j) ∈ L, k ∈ K (4.4)
 61

X
Q(u, z̃) = min dP
i,t hi,t (4.5)
y,f,p,q,v,h
i,t
X
s. t. yi,j,t = 1 − z̃i,j,k,t ui,j,k , ∀(i, j) ∈ L, t (4.6)
k
X X
P P
fi,j,t = fj,i,t + pi,t − hi,t dP
i,t , ∀i, t (4.7)
j j

Q Q
+ qi,t − hi,t dQ
X X
fi,j,t = fj,i,t i,t , ∀i, t (4.8)
j j
P + x fQ
ri,j fi,j,t
max i,j i,j,t
−V (1 − yi,j,t ) ≤ vj,t − vi,t +
V0
≤ V max (1 − yi,j,t ), ∀(i, j) ∈ L, t (4.9)

0 ≤ pi,t ≤ Pimax , ∀i, t (4.10)

0 ≤ qi,t ≤ Qmax
i , ∀i, t (4.11)
P
− M yi,j,t ≤ fi,j,t ≤ M yi,j,t , ∀(i, j), t (4.12)
Q
− M yi,j,t ≤ fi,j,t ≤ M yi,j,t , ∀(i, j), t (4.13)

V min ≤ vi,t ≤ V max , ∀i, t (4.14)

0 ≤ hi,t ≤ 1 ∀i, t (4.15)

yi,j,t Binary, ∀(i, j) ∈ L, t (4.16)

In this two-stage model, the objective is to minimize the expected value of the second-

stage objective function with respect to the probability distribution of the random vector z̃,

as expressed in Equation (4.1). In the first-stage, hardening strategies are selected from the

feasibility set U defined by Equations (4.2)-(4.4). Equation (4.2) is the budget constraint. It

means the investment on hardening distribution lines should be within a predetermined budget.

Equation (4.3) means that exactly one hardening strategy should be selected. Doing nothing

is also considered a strategy in this context.

The second-stage problem, which is defined by Equations (4.5)-(4.16), is a function of the

first-stage decisions u and the uncertainty parameter z. Equation (4.5) represents the second-

stage objective function—total load shedding. Equation (4.6) connects the operational status

of power lines with the hardening strategy and the damaging status of power lines. It selects the

appropriate status of a power line based on the first-stage solution and uncertainty realization.
 62

Equation (4.7) and (4.8) means the real and reactive power flowing out of a node must be equal

to the power flowing into a node. Equation (4.9) enforces the voltage differences caused by

power flow between two connected nodes. Equations (4.7)-(4.9) are called DistFlow equations

and have been widely used in distribution network planning and operation problems [66, 65].

Equations (4.10) and (4.11) are the generator capacity constraints. Equations (4.12) and (4.13)

limit the power flow on distribution lines based on the operation status of lines (If a line is

damaged, there would be no power flow on it). Equation (4.14) limits the range of voltage

levels. Equation (4.15) limits the percentage levels of unserved load.

4.4 Solution Method

In the two-stage stochastic programming model represented by Equation (4.1), the random

data vector z̃ has finite support. However, the possible number of scenarios for the uncertainty

parameter is 2|L| and it is computationally intractable to solve the exact two-stage stochastic

model. In this paper, we use sample average approximation (SAA) [41] to reduce the problem

size through random sampling of the uncertainty parameters. In order for the approximation to

achieve reasonable accuracy, a decent number of scenarios are needed. The resulting problem

would still be large, which makes solving the problem directly difficult. We use the L-shaped

method (based on Benders Decomposition) to further reduce the computation time.

4.4.1 Sample average approximation

The idea of SAA is rather simple. When SAA is applied to solve the two-stage stochastic

program in Equation (4.1), instead of solving it directly, we estimate the optimal solution

by sampling. More specifically, a sample z 1 , . . . , z |S| of |S| realizations of the random vector

z̃ is generated through Monte Carlo simulation. Assuming the sample is independent and

identically distributed (i.i.d.), the expected value function E[Q(u, z̃)] is approximated by the
1 P s
sample average function |S| s∈S Q(u, z ). As a result, the original two-stage problem can be

approximated by the following problem.

1 X
min Q(u, z s ) (4.17)
u∈U |S|
s∈S
 63

Given the scenario set S, the model in Equation (4.17) is deterministic and can be solved

with a deterministic optimization algorithm.

If we use (û, Q̂) to represent the optimal solution and objective value of the SAA model

and use (u∗ , Q∗ ) to represent the optimal solution and objective value of the original stochastic

model in Equation (4.1), we have Q̂ → Q∗ as |S| → ∞. Detailed theoretical discussion on SAA

is available in [41].

4.4.2 Single replication procedure

SAA is a simulation based method. As the name suggests, it can only provide an ap-

proximation of the optimal solution. Thus, it is important to test the solution quality to see

how good the approximation is. One good way to test the quality of stochastic programming

solutions is the multiple replication procedure [10]. The drawback of the multiple replication

procedure is that many sets of samples need to be generated (usually more than 30) and the

corresponding stochastic program solved, which is undesirable when even solving one program

takes a long time. To circumvent this problem, the single replication procedure outlined by

Bayraksan et al. [9] is used instead.

Use û to denote the current optimal solution, which is the object for test; and g(û) to

denote the optimality gap between E[Q(û, z̃)] and minE[Q(u, z̃)] The procedure is described as
u∈U
follows:

Step 1 : Generate a new sample of scenarios S̄ of the random vector z̃, with z̄ = {z̄ 1 , z̄ 2 , . . . , z̄ |S̄| }.

Step 2 : Solve the following stochastic program, and denote the optimal solution as u∗ .

1 X
min Q(u, z̄ s ) (4.18)
u∈U |S̄|
s∈S̄

Step 3 : Calculate the sample mean of the optimality gap ĝ and sample variance s2 (u∗ ) as

follows:

1 X
ĝ(û) = Q(u∗ , z̄ s ) − Q(û, z̄ s ) (4.19)
|S̄|
s∈S̄
 2
2 1 X s ∗ s
σ̂ (û) = (Q(û, z̄ ) − Q(u , z̄ )) − ĝ(û) (4.20)
|S̄| − 1
s∈S̄
 64


Step 4 : Calculate the one-sided confidence interval of the optimality gap g(û) as 0, ĝ(û) +
p 
t|S̄|−1,α σ̂(û)/ |S̄| , where t|S̄|−1,α is the 1 − α quantile of the Students’ t distribution

with |S̄| − 1 degrees of freedom.

This single replication procedure can provide an interval estimation of the optimality gap

between the SAA approximate and the optimal solution, based on which we can judge the

solution quality.

4.4.3 L-shaped method

Even with the reduced scenario number brought by SAA, the optimization problem can

be very large and difficult to solve. Currently, the most frequently used algorithms to solve

large scale stochastic programming include progressive hedging [43] and traditional L-shaped

method [64] based on Benders decomposition. Since Benders decomposition relies on the dual

information of the second-stage variables, the second-stage problem should be continuous. In

this paper, the integer constraints on the second-stage variable y can be relaxed, which makes

Benders decomposition a viable option. Another advantage of the L-shaped method is that it

is relatively easy to implement and performs pretty well.

In the L-shaped method, the large scale optimization problem is decomposed into a master

problem and a series of subproblems. The master problem and subproblems are formulated as

follows:

Master Problem:

min θ (4.21)

s. t. u ∈ U (4.22)

θ ≥ f ω (u), ∀ω = 1, . . . , Ω (4.23)

Subproblems:

Q(û, z s ), ∀s ∈ S (4.24)
 65

Optimality cuts represented by Equation (4.23) are iteratively generated until optimal so-

lution is identified. The optimality cuts can be formulated as follows:


1 X X ω
ui,j,k ) − µ−,ω
X
s max
θ≥ [αi,j,t,s (1 − zi,j,k i,j,t,s V
|S| s
i,j,t k

−,ω
− µ+,ω
X
max ω
i,j,t,s V ]− (βi,t,s Pimax + γi,t,s
ω
Qmax
i − ηi,t,s
i,t

min +,ω max
V + ηi,t,s V + λωi,t,s ) (4.25)

ω
where αi,j,t,s is the dual variable for Equation (4.6) corresponding to scenario s at iteration ω.

µ−,+,ω ω ω
i,j,t,s represents the dual variables of Equation (4.9). βi,t,s and γi,t,s are the dual variables
−,+,ω
for Equations (4.10) and (4.11). ηi,t,s and corresponds to Constraint (4.14). λωi,t,s is the dual

variable for Equation (4.15).

The master problem is equivalent to the full problem if the optimality cuts in Equation

(4.23) include cuts generated from all possible combinations of feasible dual variables from

the subproblems. Since the optimality cuts are iteratively generated, the master problem

represents a relaxation of the SAA approximation model and can provide an “optimistic”

estimation (lower bound) of the optimal solution as well as a candidate first-stage solution û.

The candidate solution is then passed to the subproblems. After solving the subproblems, we

can check if the estimation from the master problem is realistic or not. If it is not realistic, then

the subproblems can send a feedback message in the form of optimality cuts in Equation (4.23)

to the master problem. In addition, the average of the subproblem objective values provides

an upper bound of the optimal solution. No feasibility cut is generated since the problem has

complete recourse, which means the second-stage problems are always feasible regardless of

which first-stage solution is selected.

Detailed steps of the algorithm are described as follows.

Step 0 : Initialize by setting lower bound LB = −∞, U B = ∞, and iteration number ω = 0.

Go to Step 1.

Step 1 : Solve the master problem (4.21)-(4.23) and let (û, θ̂) denote its optimal solution.

Update the lower bound as LB = θ̂. Go to Step 2.

 66

Step 2 : Solve the subproblems Q(û, z s ), ∀s ∈ S. Update the iteration number ω ← ω + 1,

1 P s
and the upper bound U B = |S| s∈S Q(û, z ). Go to Step 3.

Step 3 : If U B − LB ≥ , generate an optimality cut based on solutions of the subproblems

and go to Step 1; otherwise return û as the optimal hardening strategy and U B as the

optimal value.

4.4.4 Greedy search heuristic

With the Knapsack constraints, large number of scenarios and 24 hour operation consider-

ation, this problem still takes a large amount of time to solve. In this section, we develop a

heuristic method as an alternative to exact methods like the Benders Decomposition. In this

algorithm, the two-stage model is simplified into a Knapsack problem, which then can be solved

by dynamic programming [47]. The algorithm is briefly described as follows:

Step 0 : Normalize the budget and hardening costs: divide the budget and hardening costs

by min cH
i,j,k . Denote the normalized costs as wi,j,k and normalized budget as W .

1 P s )−
Step 1 : Calculate the benefit of hardening each individual lines defined as vi,j,k = |S| s [Q(0, z

Q(u0i,j,k , z s )], where u0i,j,k means only line (i, j) is hardened with method k.

Step 2 : Solve the following Knapsack problem:

X
max vi,j,k ui,j,k (4.26)
i,j,k
X
s. t. wi,j,k ui,j,k ≤ W (4.27)
i,j,k
Equations(4.3) − (4.4) (4.28)

This heuristic does not consider the additional benefits the system can obtain through

coordination between different lines. Therefore the solution it yields is usually not optimal.

However, it can provide a reasonably good solution almost instantly, which can be desirable

when solving industry scale problems. In addition, this algorithm can be widely applied to

problems with Knapsack constraints.

 67

4.5 Case Study

In this section, we present our case study and the numerical results. The model and

algorithm is tested on an IEEE-33 bus system. The performance of the distribution network

in the 24-hour period during a hurricane is used to evaluate the hardening plans. As such,

the effect of hurricane decay, different time points a component could fail, and component

repair/recovery can be accounted for.

4.5.1 Scenario generation

For overhead power lines, severe weather poses the most significant threat [16]. In this

paper, we focus on hurricanes for illustration. Note scenarios are treated as input data for the

model. Other natural disasters can be added without affecting the model.

For power lines in distribution networks hit by a hurricane, failure typically occurs for

several reasons:

ˆ Distribution poles are toppled by strong winds.

ˆ The conductors between distribution poles are damaged by direct wind force impact.

ˆ The conductors between distribution poles are damaged by trees uprooted or broken by

strong winds.

Due to the lack of data, we only consider the first two reasons. In addition, since in this

paper we mainly consider the hardening of distribution lines, damages to other components of

a distribution networks including substations and distribution transformers are not considered.

We use ppole and pconductor to represent the probability of failure caused by the afore-

mentioned two reasons. Using the component fragility model proposed in [55] and [30], the

probabilities can be formulated as follows:

ppole =min{aebw , 1} (4.29)


0,


 w ≤ wmin

pconductor = w−wmin (4.30)
 wmax −wmin ,wmin ≤ w ≤ wmax


w ≥ wmax

1,
 68

where w represents wind speed. (a, b) are parameters related to pole property. wmax represents

the maximum wind speed the conductor material can withstand and wmin represents the min-

imum wind speed that can affect the conductor material. The restoration time for failed poles

and conductors are assumed to be 6 and 4 hours respectively.

The model can handle multiple hardening strategies. However, to simplify computation and

analysis, we only consider the following hardening strategy—updating the pole and conductor

of a distribution line. Hardening strategies affect one or more parameters in the fragility model

in (4.29)-(4.30), with the effect of reducing failure probability of the affected components. For

example, with new conductor material, both wmin and wmax could increase.

Wind speed information can be obtained from the ASCE-7 windspeed website [1]. A distri-

bution network typically covers a relatively small geographical area. Therefore in this paper,

we assume the wind speed experienced by different segments of the entire distribution network

is the same at a single moment. Based on the model proposed in [39], the wind speed of a

hurricane decays after landing. The decay pattern is described by Equation (4.31).

V (t) = Vb + (RV0 − Vb )e−αt (4.31)

where R = 0.9 accounts for the sea-land wind speed reduction, Vb = 13.75m/s, α = 0.095h−1 ,

and V0 is the wind speed at the time of landfall.

With the wind speed information and component fragility probability, we can simulate the

network failure occurrence via Monte Carlo simulation. For each pole or segment of conductor,

if they are not damaged in the previous hours, then their failure statuses follow independent

Bernoulli distributions with probabilities calculated from (4.29) and (4.30). If they are already

damaged, then they will remain damaged until having been repaired. Since each distribution

line in the test case consists of multiple poles and segments of conductors, their failure status

can be obtained by aggregating the failure status of individual poles and conductors. The

reparation of multiple components in the system is assumed to be conducted simultaneously.

In addition, the repair time for each pole and conductor is assumed to be constant. However,

due to the different failure time of the components in a distribution line, the repair time

experienced by each line is variable. The typical damage condition scenario of a distribution
 69

Figure 4.1 Typical Scenario

line in 24 hours is demonstrated in Figure 4.1. In this scenario, the line remains damaged for

13 hours if it is not hardened. If it is hardened, then it will only be damaged for 6 hours.

4.5.2 Experiment results

The IEEE 33-bus distribution system used in this paper is shown in Figure 4.2. In this

system, there are 37 distribution lines, five of which are switchable tie lines. Detailed network

and load data can be found in [8]. The length of each power line is assumed to be proportional

to the resistance of the line. Then we can estimate the number of poles and conductors

represented by each line in the network. Assuming the span between 2 poles is 150 feet, the

cost of hardening each component is summarized in Table 4.1. With this information we can

calculate the cost of hardening each line accordingly. The 24-hour load profile is displayed in

Figure 4.3. In addition, we assume distributed generators (DG) exist at nodes 7, 10, 13, 19,

23, 26 and 30. The complete scenario space has 237 scenarios. Based on the procedure outlined

in Section 4.5.1, 100 scenarios are randomly generated. These scenarios represent the damage

conditions of the power lines under different states of hardening.

The solution methods were implemented with AMPL and MATLAB on a desktop computer

with Intel Core i7 3.4 GHz CPU, 8 GB memory. The Benders Decomposition algorithm takes

several hours to converge, while the heuristic can yield a solution instantaneously.
 70

Figure 4.2 IEEE 33-Bus Test Case

Table 4.1 Hardening costs

Strategy Cost ($/unit) Unit
Update poles 6,000 N/A
Update conductors 700,000 mile
Update both 910,000 mile

When budget is 7 million dollars, 19 lines are hardened. The hardened lines are highlighted

in Figure 4.4. We can see many lines on the main branch are strengthened. In addition, one

tie-line connecting two branches is hardened, forming a large loop of strengthened lines. DGs

also have noticeable effect on the selection of hardening strategies—most of the lines adjacent

to DGs are hardened.

In Figure 4.5, the performances of three solutions under different scenarios are summarized.

We call the three solutions Wait and See (WS) solution, Optimal Recourse Problem (RP)

solution, and “None” solution. WS solution is obtained by solving the two-stage problem for

each scenario and then taking the expected value of the objectives. WS solutions represent

the best possible investment decisions a planner can make if the perfect information about

future uncertainty is available. In other words, planners know exactly which lines are going to

be damaged and when such damages would occur, and can change their investment decisions

accordingly. This solution is unattainable and can only provide a lower bound on the objective
 71

Figure 4.3 24-Hour Load Multiplier

value. The RP solution is obtained through the solution procedure outlined in Section 4.4.

The “None” solution means no investment. We can see that although the RP solution has

worse performances than the WS solutions, it is considerably better than no investment, which

highlights the importance of investment in network hardening.

By varying the amount of available budget, we can examine the effect of available budget on

the performance of resulting hardening plans. It is obvious that the amount of load shedding

would decrease and budget increases. However, by looking at the trend, we can gauge the

marginal benefit of increasing the investment budget. The performance of the hardening plans

under different budgets are summarized in Table 4.2. The difference between WS and RP

values is called the Expected Value of Perfect Information (EVPI). It represents the benefit of

having access to perfect information about uncertainty. We can see EPVI is roughly decreasing

with the increase of budgets. The heuristic solution is also included in the table. We can see

the heuristic method can provide pretty good solutions.

4.5.3 Solution validation

In this paper, we test the quality of solution by calculating the confidence interval of opti-

mality gaps under different budgets. The results are summarized in Table 4.3. The width of the

confidence intervals vary, but are typically small (less than 6% of the objective value). This sug-
 72

Figure 4.4 Hardening Plan (Budget = 7m$)

Table 4.2 Avergae Load Shedding under different budgets

Budget ($m) WS RP EVPI Heuristic
5 2.4445 3.0355 0.5910 3.1126
6 2.3105 2.9164 0.6059 2.9578
7 2.2029 2.8245 0.6216 2.8869
8 2.1187 2.6659 0.5472 2.6869
9 2.0548 2.5672 0.5124 2.5927
10 2.0050 2.4816 0.4766 2.5439
11 1.9753 2.4455 0.4702 2.5187
12 1.9577 2.4287 0.4710 2.5111

gests that with 100 scenarios, the SAA method can provide a reasonably good approximation

to the optimal solution.

4.6 Conclusion

In this paper, a new two-stage stochastic programming model is proposed for improving the

resilience of a distribution network with hardening. In this problem, the uncertain parameters

concerning the condition of a network affected by natural disasters are considered. The average

load shedding under different scenarios are optimized. To simplify the computation, sample

average approximation is used to reduce the problem size. In addition, Benders decomposition

is applied to speed up the computation. A heuristic method is developed to yield a good solu-
 73

Figure 4.5 Histogram of Load Shedding Under Different Plans (Budget = 7m$)

Table 4.3 Confidence interval of the optimality gap under different budgets
Budget ($m) 95% C. I. Budget ($m) 95% C. I.
5 [0, 0.1004] 9 [0, 0.0576]
6 [0, 0.0281] 10 [0, 0.0783]
7 [0, 0.0739] 11 [0, 0.0741]
8 [0, 0.1514] 12 [0, 0.1145]

tion quickly. A case study on the IEEE-33 bus distribution system is used to test the model

and algorithms. The experiment results are then tested for quality with the single replica-

tion procedure. The solution validation procedure demonstrates that the SAA can provide a

reasonably good approximation to the actual optimal solution.

For future research, we will focus on the following directions. Firstly, better scenario genera-

tion methods can be developed to take into consideration the dependencies between components

and the randomness in the recovery time of damaged components. In addition, other disaster

scenarios like earthquakes and floods can be explored. Secondly, more hardening strategies

and preparedness measures such as training of personnel and emergency drill to improve the

response time of repair crew can be included in the model for a more flexible solution. Thirdly,

we can take into consideration other emergency operation measures such as intentional island-

ing to further improve the solution. Last but not least, more sophisticated solution methods
 74

might be needed to solve larger scale problems.

 75

CHAPTER 5. CONCLUSION AND DISCUSSION

In my dissertation I present three applications of mathematical programming under uncer-

tainty methods on different types of power system planning problems. The first two papers

focus on the application of robust optimization to transmission expansion planning. The third

paper presents an application of stochastic programming to power distribution network hard-

ening planning. The detailed contributions of these works and conclusions are discussed as

follows.

Transmission expansion planning facilitates generation expansion and is critical to the reli-

able and economic operations of power systems. The long planning horizon means there are a

various sources of uncertainty. In the first two papers, robust optimization is used to address

the challenges posed by generation expansion uncertainty. In the first paper, in addition to

generation expansion and retirement, load forecast is also considered as a source of uncertainty.

Two optimization criteria: minimax cost and minimax regret are compared under the robust

optimization paradigm. With this approach, a transmission plan that is robust and has good

performances under all scenarios can be derived. The resulting models are formulated as trilevel

mixed-integer programming models. Column and constraint generation is used to decompose

the model into a master problem and a subproblem, where new constraints and variables are

iteratively generated by the subproblems to be fed back to the master problem. The subprob-

lems themselves are bilevel mixed-integer programs, which are reformulated into single level

mixed-integer programs with the KKT optimality conditions. The models and algorithm are

tested on an IEEE 118-bus system. The results of the MMR and MMC models are analyzed

and compared. We conclude that both criteria may outperform each other depending on the

uncertainty set. Thus, the optimization criteria should be selected based on available data and

decision-makers’ preferences.
 76

As for the second paper, a robust transmission expansion planning model that can explore

the impact of generation expansion and retirement uncertainty in more detail is proposed. We

use binary variables to describe whether a generator exists after generation expansion or retire-

ment in the future. With the new modeling framework, the uncertainty set can contain detailed

description of the uncertain parameters including the location, size and time of commission.

The same column and constraint generation procedure is used to solve the trilevel model. How-

ever, by utilizing the special structure of the subproblem, the bilinear subproblem can be solved

to optimality relatively quickly. A WECC 240-bus test system is used to demonstrate the ef-

fectiveness of the model and algorithm. The performance of the robust transmission plans is

compared to that of the plans derived from the deterministic model and heuristic rules under

different scenarios of fuel price and environmental policies. The results show that the robust

plans not only out-perform plans obtained by other methods, but also remain robust under

different scenarios.

In the third paper, a distribution network planning problem is solved with stochastic pro-

gramming. To improve system resilience under extreme weather events, the power industry

is investing to hardening the distribution networks by upgrading the infrastructure. With the

limited budget, it is impossible to harden the entire network indiscriminately. A model is de-

veloped to optimally allocate the resources to the most critical or vulnerable components in

the network. In this problem, the uncertainty lies in the condition of the network after it is

affected by a natural disaster. Sample average approximation is used to reduce the problem

size and Benders decomposition is used to accelerate the computation. A heuristic algorithm

is also developed to yield a good solution quickly. A case study on the IEEE 33-bus system is

used to test the model and algorithms. The solution is then validated via the single replication

procedure, which shows that SAA can provide a reasonably good approximation to the optimal

solution.

There exist several interesting directions for future work. In terms of transmission expansion

planning, other sources of uncertainty and methods of dealing with them can be incorporated

into the model. For example, since the probability distributions of load and renewable genera-

tion might be obtained from historical data, it is possible to combine stochastic programming
 77

and robust optimization by modeling different sources of uncertainty with different methods.

Such a problem would be even more challenging to solve with the increased problem size. New

algorithms and heuristics needs to be developed. In addition, more details including more

explicit representations of various policies, demographic changes and climate changes can be

incorporated into the existing model. In terms of distribution network hardening and resilience

planning, more hardening strategies can be added to the model, including preparedness mea-

sures such as training of personnel and emergency drill to improve the response time of repair

crew. Other disaster scenarios like earthquakes and floods can also be considered.
 78

BIBLIOGRAPHY

[1] ASCE 7 windspeed.

[2] IEEE 118-bus, 54-unit, 24-hour system.

[3] Tohid Akbari, Ashkan Rahimi-Kian, and Mohammad Tavakoli Bina. Security-constrained

transmission expansion planning: A stochastic multi-objective approach. International

Journal of Electrical Power & Energy Systems, 43(1):444–453, 2012.

[4] Natalia Alguacil, Andrés Delgadillo, and José M Arroyo. A trilevel programming approach

for electric grid defense planning. Computers & Operations Research, 41:282–290, 2014.

[5] B. Alizadeh, S. Dehghan, N. Amjady, S. Jadid, and A. Kazemi. Robust transmission sys-

tem expansion considering planning uncertainties. Generation, Transmission Distribution,

IET, 7(11):1318–1331, November 2013.

[6] J.M. Arroyo. Bilevel programming applied to power system vulnerability analysis under

multiple contingencies. Generation, Transmission Distribution, IET, 4(2):178–190, Febru-

ary 2010.

[7] José Manuel Arroyo, Natalia Alguacil, and Miguel Carrión. A risk-based approach for

transmission network expansion planning under deliberate outages. Power Systems, IEEE

Transactions on, 25(3):1759–1766, 2010.

[8] M.E. Baran and F.F. Wu. Network reconfiguration in distribution systems for loss re-

duction and load balancing. Power Delivery, IEEE Transactions on, 4(2):1401–1407, Apr

1989.
 79

[9] Güzin Bayraksan and David P Morton. Testing solution quality in stochastic program-

ming: a single replication procedure. In Proceedings of the 16th Symposium of IASC on

Computational Statistics. Physica-Verlag/Springer, Prague, Czech Republic, 2004.

[10] Güzin Bayraksan and David P Morton. Assessing solution quality in stochastic programs

via sampling. Tutorials in Operations Research, 5:102–122, 2009.

[11] Aharon Ben-Tal and Arkadi Nemirovski. Robust solutions of uncertain linear programs.

Operations research letters, 25(1):1–13, 1999.

[12] D. Bertsimas, E. Litvinov, X.A. Sun, Jinye Zhao, and Tongxin Zheng. Adaptive robust

optimization for the security constrained unit commitment problem. Power Systems, IEEE

Transactions on, 28(1):52–63, 2013.

[13] D. Bertsimas and M. Sim. The price of robustness. Operations research, 52(1):35–53, 2004.

[14] Dimitris Bertsimas and Melvyn Sim. Robust discrete optimization and network flows.

Mathematical Programming, 98(1-3):49–71, 2003.

[15] Stephen Poythress Boyd and Lieven Vandenberghe. Convex optimization. Cambridge

university press, 2004.

[16] RE Brown. Cost-benefit analysis of the deployment of utility infrastructure upgrades and

storm hardening programs. Quanta Technology, Raleigh, 2009.

[17] M Oloomi Buygi, Gerd Balzer, H Modir Shanechi, and Mohammad Shahidehpour.

Market-based transmission expansion planning. Power Systems, IEEE Transactions on,

19(4):2060–2067, 2004.

[18] Miguel Carrión, José Manuel Arroyo, and Natalia Alguacil. Vulnerability-constrained

transmission expansion planning: A stochastic programming approach. Power Systems,

IEEE Transactions on, 22(4):1436–1445, 2007.

[19] Metin Celebi, Frank Graves, Gunjan Bathla, and Lucas Bressan. Potential coal plant

retirements under emerging environmental regulations. The Brattle Group, December, 8,

2010.
 80

[20] Bokan Chen, Jianhui Wang, Lizhi Wang, Yanyi He, and Zhaoyu Wang. Robust opti-

mization for transmission expansion planning: Minimax cost vs. minimax regret. Power

Systems, IEEE Transactions on, 29(6):3069–3077, November 2014.

[21] C. Chen, J. Wang, F. Qiu, and D. Zhao. Resilient distribution system by microgrids

formation after natural disasters. Smart Grid, IEEE Transactions on, PP(99):1–1, 2015.

[22] Lichun Chen and Elise Miller-Hooks. Resilience: An indicator of recovery capability in

intermodal freight transport. Transportation Science, 46(1):109–123, 2012.

[23] Teofilo De la Torre, James W Feltes, Tomas Gomez San Roman, and Hyde M Merrill.

Deregulation, privatization, and competition: Transmission planning under uncertainty.

Power Systems, IEEE Transactions on, 14(2):460–465, 1999.

[24] US EIA. Annual energy outlook 2015. US Energy Information Administration, Washing-

ton, DC, 2015.

[25] Laurent El Ghaoui and Hervé Lebret. Robust solutions to least-squares problems with

uncertain data. SIAM Journal on Matrix Analysis and Applications, 18(4):1035–1064,

1997.

[26] Lina P Garcés, Antonio J Conejo, Raquel Garcı́a-Bertrand, and Rubén Romero. A bilevel

approach to transmission expansion planning within a market environment. Power Sys-

tems, IEEE Transactions on, 24(3):1513–1522, 2009.

[27] Mohammad R Hesamzadeh, Darryl R Biggar, Nasser Hosseinzadeh, and Peter J Wolfs.

Transmission augmentation with mathematical modeling of market power and strategic

generation expansion: Part I. Power Systems, IEEE Transactions on, 26(4):2040–2048,

2011.

[28] Jing Hu, John E Mitchell, Jong-Shi Pang, Kristin P Bennett, and Gautam Kunapuli. On

the global solution of linear programs with linear complementarity constraints. SIAM

Journal on Optimization, 19(1):445–471, 2008.

 81

[29] R.A. Jabr. Robust transmission network expansion planning with uncertain renewable

generation and loads. Power Systems, IEEE Transactions on, 28(4):4558–4567, 2013.

[30] Pirooz Javanbakht and Salman Mohagheghi. A risk-averse security-constrained optimal

power flow for a power grid subject to hurricanes. Electric Power Systems Research,

116:408–418, 2014.

[31] R. Jiang, J. Wang, M. Zhang, and Y. Guan. Two-stage minimax regret robust unit

commitment. Power Systems, IEEE Transactions on, PP(99):1–12, 2013.

[32] Ruiwei Jiang, Jianhui Wang, and Yongpei Guan. Robust unit commitment with wind

power and pumped storage hydro. Power Systems, IEEE Transactions on, 27(2):800–810,

2012.

[33] Shan Jin and Sarah M Ryan. A tri-level model of centralized transmission and decentralized

generation expansion planning for an electricity market: Part I. Power Systems, IEEE

Transactions on, 29(1):132–141, 2014.

[34] Shan Jin and Sarah M Ryan. A tri-level model of centralized transmission and decentralized

generation expansion planning for an electricity market: Part II. Power Systems, IEEE

Transactions on, 29(1):142–148, 2014.

[35] Shan Jin and S.M. Ryan. Capacity expansion in the integrated supply network for an

electricity market. Power Systems, IEEE Transactions on, 26(4):2275–2284, 2011.

[36] Panida Jirutitijaroen and Chanan Singh. Reliability constrained multi-area adequacy plan-

ning using stochastic programming with sample-average approximations. Power Systems,

IEEE Transactions on, 23(2):504–513, 2008.

[37] Tore W Jonsbråten, Roger JB Wets, and David L Woodruff. A class of stochastic programs

withdecision dependent random elements. Annals of Operations Research, 82:83–106, 1998.

[38] Peter Kall and Stein W Wallace. Stochastic programming. John Wiley and Sons Ltd, 1994.
 82

[39] John Kaplan and Mark DeMaria. A simple empirical model for predicting the decay of

tropical cyclone winds after landfall. Journal of applied meteorology, 34(11):2499–2512,

1995.

[40] Amin Khodaei and Mohammad Shahidehpour. Microgrid-based co-optimization of gen-

eration and transmission planning in power systems. Power Systems, IEEE Transactions

on, 28(2):1582–1590, 2013.

[41] Anton J Kleywegt, Alexander Shapiro, and Tito Homem-de Mello. The sample average ap-

proximation method for stochastic discrete optimization. SIAM Journal on Optimization,

12(2):479–502, 2002.

[42] I. Konstantelos and G. Strbac. Valuation of flexible transmission investment options under

uncertainty. Power Systems, IEEE Transactions on, 30(2):1047–1055, March 2015.

[43] Arne Løkketangen and David L Woodruff. Progressive hedging and tabu search applied to

mixed integer (0, 1) multistage stochastic programming. Journal of Heuristics, 2(2):111–

128, 1996.

[44] Juan Álvarez López, Kumaraswamy Ponnambalam, and Vı́ctor H Quintana. Generation

and transmission expansion under risk using stochastic programming. Power Systems,

IEEE Transactions on, 22(3):1369–1378, 2007.

[45] Juan Álvarez López, Kumaraswamy Ponnambalam, and Vı́ctor H Quintana. Generation

and transmission expansion under risk using stochastic programming. Power Systems,

IEEE Transactions on, 22(3):1369–1378, Aug 2007.

[46] Pouria Maghouli, Seyed Hamid Hosseini, Majid Oloomi Buygi, and Mohammad Shahideh-

pour. A scenario-based multi-objective model for multi-stage transmission expansion plan-

ning. Power Systems, IEEE Transactions on, 26(1):470–478, 2011.

[47] Silvano Martello, David Pisinger, and Paolo Toth. Dynamic programming and strong

bounds for the 0-1 knapsack problem. Management Science, 45(3):414–424, 1999.
 83

[48] Elise Miller-Hooks, Xiaodong Zhang, and Reza Faturechi. Measuring and maximizing re-

silience of freight transportation networks. Computers & Operations Research, 39(7):1633–

1643, 2012.

[49] James T Moore and Jonathan F Bard. The mixed integer linear bilevel programming

problem. Operations Research, pages 911–921, 1990.

[50] A. Moreira, A. Street, and J.M. Arroyo. An adjustable robust optimization approach for

contingency-constrained transmission expansion planning. Power Systems, IEEE Trans-

actions on, 30(4):2013–2022, July 2015.

[51] Amir Motamedi, Hamidreza Zareipour, Majid Oloomi Buygi, and William D Rosehart. A

transmission planning framework considering future generation expansions in electricity

markets. Power Systems, IEEE Transactions on, 25(4):1987–1995, 2010.

[52] Francisco D Munoz, Benjamin F Hobbs, Jonathan L Ho, and Saamrat Kasina. An

engineering-economic approach to transmission planning under market and regulatory un-

certainties: WECC case study. Power Systems, IEEE Transactions on, 29(1):307–317,

2014.

[53] K.A. Nigim and Y.G. Hegazy. Intention islanding of distributed generation for reliability

enhancement. In Power Engineering Society General Meeting, 2003, IEEE, volume 4, page

2451 Vol. 4, July 2003.

[54] Executive Office of the President. Economic Benefits of Increasing Electric Grid Resilience

to Weather Outages - August 2013. IEEE USA Books & eBooks, 2013.

[55] Min Ouyang and Leonardo Dueñas-Osorio. Multi-dimensional hurricane resilience assess-

ment of electric power systems. Structural Safety, 48:15–24, 2014.

[56] F. Pilo, G. Celli, and S. Mocci. Improvement of reliability in active networks with inten-

tional islanding. In Electric Utility Deregulation, Restructuring and Power Technologies,

2004. (DRPT 2004). Proceedings of the 2004 IEEE International Conference on, volume 2,

pages 474–479 Vol.2, April 2004.

 84

[57] David Pozo, Enzo E Sauma, and Javier Contreras. A three-level static milp model for

generation and transmission expansion planning. Power Systems, IEEE Transactions on,

28(1):202–210, 2013.

[58] J.E. Price and J. Goodin. Reduced network modeling of wecc as a market design prototype.

In Power and Energy Society General Meeting, 2011 IEEE, pages 1–6, 2011.

[59] Jae Hyung Roh, Mohammad Shahidehpour, and Lei Wu. Market-based generation

and transmission planning with uncertainties. Power Systems, IEEE Transactions on,

24(3):1587–1598, 2009.

[60] C Ruiz and AJ Conejo. Robust transmission expansion planning. European Journal of

Operational Research, 242(2):390–401, 2015.

[61] S Ryan, J McCalley, and D Woodruff. Long term resource planning for electric power

systems under uncertainty. Technical report, Iowa State University, 2011.

[62] A. Street, F. Oliveira, and J.M. Arroyo. Contingency-constrained unit commitment with

security criterion: A robust optimization approach. Power Systems, IEEE Transactions

on, 26(3):1581–1590, 2011.

[63] Maarten K Van Aalst. The impacts of climate change on the risk of natural disasters.

Disasters, 30(1):5–18, 2006.

[64] Richard M Van Slyke and Roger Wets. L-shaped linear programs with applications to

optimal control and stochastic programming. SIAM Journal on Applied Mathematics,

17(4):638–663, 1969.

[65] Zhaoyu Wang, Bokan Chen, Jianhui Wang, and M.M. Begovic. Stochastic dg placement for

conservation voltage reduction based on multiple replications procedure. Power Delivery,

IEEE Transactions on, 30(3):1039–1047, June 2015.

[66] Zhaoyu Wang, Bokan Chen, Jianhui Wang, M.M. Begovic, and Chen Chen. Coordinated

energy management of networked microgrids in distribution systems. Smart Grid, IEEE

Transactions on, 6(1):45–53, Jan 2015.

 85

[67] Pan Xu and Lizhi Wang. An exact algorithm for the bilevel mixed integer linear program-

ming problem under three simplifying assumptions. Computers & Operations Research,

2013.

[68] Emre Yamangil, Russell Bent, and Scott Backhaus. Designing resilient electrical distribu-

tion grids. arXiv preprint arXiv:1409.4477, 2014.

[69] Yiming Yao, Thomas Edmunds, Dimitri Papageorgiou, and Rogelio Alvarez. Trilevel opti-

mization in power network defense. Systems, Man, and Cybernetics, Part C: Applications

and Reviews, IEEE Transactions on, 37(4):712–718, 2007.

[70] H Yu, CY Chung, KP Wong, and JH Zhang. A chance constrained transmission network

expansion planning method with consideration of load and wind farm uncertainties. Power

Systems, IEEE Transactions on, 24(3):1568–1576, 2009.

[71] W. Yuan, J. Wang, F. Qiu, C. Chen, C. Kang, and B. Zeng. Robust optimization-based

resilient distribution network planning against natural disasters. IEEE Transactions on

Smart Grid, PP(99):1–10, 2016.

[72] Wei Yuan, Long Zhao, and Bo Zeng. Optimal power grid protection through a defender–

attacker–defender model. Reliability Engineering & System Safety, 121:83–89, 2014.

[73] Bo Zeng and Long Zhao. Solving two-stage robust optimization problems using a column-

and-constraint generation method. Operations Research Letters, 2013.

[74] Long Zhao and Bo Zeng. Robust unit commitment problem with demand response and

wind energy. In Power and Energy Society General Meeting, 2012 IEEE, pages 1–8. IEEE,

2012.