APPLICATIONS OF SENSITIVITY ANALYSIS IN PLANNING AND OPERATION OF
MODERN POWER SYSTEMS
By
Mohammed Ben-Idris

A DISSERTATION
Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of
Electrical Engineering - Doctor of Philosophy
2014

ABSTRACT
APPLICATIONS OF SENSITIVITY ANALYSIS IN PLANNING AND
OPERATION OF MODERN POWER SYSTEMS
By
Mohammed Ben-Idris
More large blackouts in power systems have been encountered in the last few years than
in any other period in history. Although the evolution of a large blackout is a complex
process, it usually starts with a small triggering event that by successive events can cascade
and evolve into a complete blackout. Therefore, hardening power systems against such events
and improving power system reliability are of necessity.
Several reliability indices for composite systems have been introduced such as loss of load
probability, loss of load frequency, loss of load expectation, expected energy not supplied,
expected demand not supplied, etc. Although these indices are of a great importance in both
planning and operational power system reliability evaluation, they lack the ability to identify
the influence of each area or equipment on system reliability. Moreover, these indices cannot
determine the most effective way to improve the system in order to keep it reliable. Due to
these reasons, significant efforts have been devoted to this area in recent years to evaluate
what the system’s reliability justifications are and to determine the best location in which
to invest.
One of the promising methods that can identify the most vulnerable component(s),
bus(es), or area(s) is performing sensitivity analysis of the possible load shedding with respect to component parameters and operating limits. The sensitivity analysis provides a
measure of the amount of load to be curtailed in response to the violation of the operating
limits or component characteristics.

The task of performing sensitivity analysis is amply used in power system reliability planning and operational studies. However, the application of sensitivity analysis in the changing
operating environment such as the newly imposed emission constraints, penetration of the
renewable energy sources and cascading failures in power systems has not been given much
attention. Sensitivity analyses of the reliability indices can be conducted either analytically
such as using state space enumeration or by means of simulations. Therefore, performing
sensitivity analyses over a large set of variables for large systems is computationally expensive
and often intractable.
The work presented in thesis develops and proposes new expressions for sensitivity analysis that can be used to identify the most vulnerable components in power systems. Also,
this work proposes a heuristic technique in conjunction with the population-based intelligent search methods (PIS) to reduce the computation burden. The work in this thesis is
divided into three parts. Part I presents the concept of the sensitivity analyses and provides the mathematical derivations of the proposed expressions. Also, this part introduces
new techniques to reduce the computational time and burden based on PIS methods and
their applications in power system reliability evaluation. Furthermore, it introduces a new
technique to classify the search space into success, failure and unclassified subspaces. Part
II presents the application of the derived expressions on power system reliability planning
studies such as the effect of each parameter on system reliability and the effect of emission
constraints on power system reliability. Also, this part examines the proposed algorithms
to speed up the computational time against the existing techniques. Part III conducts the
application of the sensitivity analysis in power system reliability operation such as mitigation
of cascading failures and prevention of the unfolding events.

Copyright by
MOHAMMED BEN-IDRIS
2014

My parents, wife and my two little daughters

v

ACKNOWLEDGMENTS

I would like to express my sincere and deepest appreciation to my advisor, Prof. Joydeep
Mitra for his guidance, valuable comments, encouragement, support and direction for this
work. It has been an honor to be one of his Ph.D. students. I would also like to thank Prof.
Shanker Balasubramaniam, Prof. Jongeun Choi and Prof. Bingsen Wang for serving on my
committee and their valuable suggestions and comments.
A word of gratitude goes to Prof. Tim Hogan and Prof. Shanker Balasubramaniam for
their support and guidance in completion of my certificate of college teaching as well as
providing funds to complete my work.
I sincerely acknowledge the funding sources that helped me to complete my Ph.D. work.
I was funded by the Libyan North American Program and the Ministry of Higher Education
and Scientific Research - Libya, my advisor and the ECE department at Michigan State
University.
I would like to thank my friends at Michigan State University and the research group
at the ERISE lab. The group has been a good advice and collaboration. In particular, I
would like to acknowledge Salem Elsaiah and Samer Sulaeman for their valuable comments,
suggestions, and collaboration.
Most of all, I appreciate the support and encouragement of my family; my father, my
mother, my wife and my brothers and sisters. I would like to thank my wife and my two
little girls for being patient and not to complain about using their time.

vi

TABLE OF CONTENTS

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

xi

LIST OF FIGURES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiv
Chapter 1 Introduction . . . .
1.1 Motivation . . . . . . . . .
1.2 Challenges . . . . . . . . .
1.3 Approaches . . . . . . . .
1.4 Contributions . . . . . . .
1.5 Organization of the Thesis

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

1
3
3
4
5
6

Chapter 2 Power Flow Modeling and Optimization Framework . . . . . . .
2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2 Power Flow Models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2.1 AC Power Flow Model . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2.2 DC Power Flow model . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2.3 Linearized Power Flow model . . . . . . . . . . . . . . . . . . . . . .
2.3 Optimization Problem for Minimum Load Curtailment . . . . . . . . . . . .
2.3.1 Network Modeling using Nonlinear Programming and the AC Power
Flow Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.3.2 Network Modeling using Linear Programming and the DC Power Flow
Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.3.3 Network Modeling using Linear Programming and the Linearized Power
Flow Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.3.3.1 Power Balance Constraints . . . . . . . . . . . . . . . . . .
2.3.3.2 Real and Reactive Power Constraints of the Generators . . .
2.3.3.3 Voltage Limits Constraints . . . . . . . . . . . . . . . . . .
2.3.3.4 Line Current Limit . . . . . . . . . . . . . . . . . . . . . . .
2.3.3.5 Linear Programming Formulation . . . . . . . . . . . . . . .
2.3.4 Other Constraints . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

8
8
9
9
10
11
13

Chapter 3 Power System Reliability Indices and Sensitivity Analysis
3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.2 Calculation of the Reliability Indices . . . . . . . . . . . . . . . . . . .
3.2.1 Calculation of System Probability Indices . . . . . . . . . . . . .
3.2.2 Calculation of System Energy Indices . . . . . . . . . . . . . . .
3.2.3 Calculation of System Frequency and Duration Indices . . . . .
3.2.4 Calculation of Bus Indices . . . . . . . . . . . . . . . . . . . . .
3.2.5 Modeling of System Components . . . . . . . . . . . . . . . . .

25
25
25
26
27
29
32
32

vii

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

14
15
17
17
18
18
19
21
23

3.3
3.4

3.2.5.1 Modeling of Available Generation . . . . . . . . . . . . . . .
3.2.5.2 Modeling of System Load . . . . . . . . . . . . . . . . . . .
3.2.5.3 Modeling of Transmission Lines . . . . . . . . . . . . . . . .
3.2.6 A Stopping Criterion . . . . . . . . . . . . . . . . . . . . . . . . . . .
Concept of the Sensitivity Analysis . . . . . . . . . . . . . . . . . . . . . . .
Sensitivity Analysis of the Reliability Indices . . . . . . . . . . . . . . . . . .
3.4.1 Sensitivity Analysis of the LOLP and LOLF Indices with Respect to
Component Parameters . . . . . . . . . . . . . . . . . . . . . . . . . .
3.4.2 Sensitivity Analysis of the EDNS index with Respect to Component
Capacities . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Chapter 4 State Space Reduction and Population-based Intelligent Search
Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.2 State Space Reduction Technique . . . . . . . . . . . . . . . . . . . . . . . .
4.3 Power System Reliability Evaluation using Binary Particle Swarm Optimization Search Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.4 Non-Sequential Monte Carlo Simulation vs. Binary Particle Swarm Optimization Search Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.5 The complementary Concept . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.5.1 Theoretical Background . . . . . . . . . . . . . . . . . . . . . . . . .
4.5.2 Mathematical Justification . . . . . . . . . . . . . . . . . . . . . . . .
4.6 Particle Swarm Optimization . . . . . . . . . . . . . . . . . . . . . . . . . .
4.6.1 Multiple Objective Functions . . . . . . . . . . . . . . . . . . . . . .
4.6.2 Single Objective Function . . . . . . . . . . . . . . . . . . . . . . . .
4.7 Dynamically-directed Binary Particle Swarm Optimization Search Method .
4.7.1 State Space Description . . . . . . . . . . . . . . . . . . . . . . . . .
4.7.2 Directed BPSO . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

33
33
33
33
34
35
35
37

39
39
41
45
46
48
49
51
52
55
55
56
57
58

Chapter 5 Sensitivity Analysis in Composite System Reliability Using Weighted
Shadow Prices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62
5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62
5.2 Shadow Price Theory and Description . . . . . . . . . . . . . . . . . . . . . . 63
5.3 Solution Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66
5.4 Case Studies . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67
5.4.1 Case Study I: IEEE RTS Base Case . . . . . . . . . . . . . . . . . . . 69
5.4.2 Case Study I: Modifed IEEE RTS . . . . . . . . . . . . . . . . . . . . 69
5.4.3 Case Study III: Increasing Generation Capacity and Loading Level . 70
5.4.4 Case Study IV: Reducing Lines Capacity . . . . . . . . . . . . . . . . 73
5.4.5 Case Study V: Calculation of the Reliability Indices of Four Test Systems 73
5.5 Results and Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74
5.6 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75

viii

Chapter 6 Power System Reliability Evaluation and Sensitivity Analysis
Using the State Space Reduction Technique and Populationbased Intelligent Search Methods . . . . . . . . . . . . . . . . . . .
6.1 Discussion on the Reduction in the Computational Time . . . . . . . . . . .
6.1.1 Reduction in Computational Time in Case of Ignoring the Unclassified
Subspace . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
6.1.2 Reduction in Computational Time in Case of Including the Unclassified Subspace . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
6.2 Solution Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
6.3 Case Studies . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
6.4 Case Study I: Calculation of the Reliability Indices . . . . . . . . . . . . . .
6.5 Case Study II: Sensitivity Analyses . . . . . . . . . . . . . . . . . . . . . . .
6.6 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Chapter 7 Reliability and Sensitivity Analysis of Composite Power Systems Under Emission Constraints . . . . . . . . . . . . . . . . . .
7.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7.2 Emissions from Electric Generation . . . . . . . . . . . . . . . . . . . . . . .
7.3 CO2 , SO2 and N Ox Emission Constraints . . . . . . . . . . . . . . . . . . .
7.4 Sensitivity Analysis of EDNS Index with respect to the Emission Limits . . .
7.5 Solution Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7.6 Case Studies . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7.7 Results and Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7.8 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

76
77
78
80
81
83
84
86
87

91
91
95
96
98
101
102
104
105

Chapter 8 Reliability and Sensitivity Analysis of Composite Power Systems Considering Voltage and Reactive Power Constraints . . . 111
8.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111
8.1.1 Calculation of Voltage and Reactive Power Limit Indices . . . . . . . 113
8.2 Sensitivity Analysis of the EDNS Index with respect to the Voltage and Reactive Power Limit Indices . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114
8.2.1 Sensitivity of the EDNS Index with Respect to Voltage Limit Constraints115
8.2.2 Sensitivity of the EDNS Index with Respect to Reactive Power Limit
Constraints . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115
8.3 Case Studies . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 116
8.3.1 Reliability Assessment . . . . . . . . . . . . . . . . . . . . . . . . . . 116
8.3.1.1 Case Study I: Gradually Increasing the Generation and Load
Levels . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118
8.3.1.2 Case Study II: Gradually Increasing the Load Level and Keeping the Generation Level Constant . . . . . . . . . . . . . . 119
8.3.2 Comparison with the DC Power Flow Model . . . . . . . . . . . . . . 120
8.3.3 Effects of Voltage and Reactive Power Limits on the Reliability Indices 121
8.3.3.1 Contributions of the Voltage and Reactive Power Constraints 122
8.3.3.2 Indices of the Voltage and Reactive Power Limits . . . . . . 123
8.3.3.3 Relaxing Voltage Constraints . . . . . . . . . . . . . . . . . 123
ix

8.4

8.3.4 Sensitivity Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . 124
Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127

Chapter 9 Hardening Power Systems Against Cascading Failures Based
on Risk Sensitivity Analysis . . . . . . . . . . . . . . . . . . . . . .
9.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
9.2 Sensitivity Analysis-based Cascading Failure Mitigation . . . . . . . . . . . .
9.3 Remedial Actions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
9.4 Case Studies . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
9.5 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

128
128
129
130
132
135

Chapter 10 Conclusions and Future Work . . . . . . . . . . . . . . . . . . . . . 137
10.1 Concluding Remarks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 137
10.2 Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 140
APPENDIX

. . . . . . . . . . . . . . .

BIBLIOGRAPHY

. . . . . . . . . . . . . . . . . . . 142

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 149

x

LIST OF TABLES

Table 5.1

Generators Shadow Prices . . . . . . . . . . . . . . . . . . . . . . .

70

Table 5.2

Transmission Lines Shadow Prices (Forward) . . . . . . . . . . . . .

71

Table 5.3

Transmission Lines Shadow Prices (Reverse) . . . . . . . . . . . . .

72

Table 5.4

Annualized Indices of the Tested Systems. . . . . . . . . . . . . . . .

73

Table 6.1

Annual Indices of the IEEE RTS . . . . . . . . . . . . . . . . . . . .

85

Table 6.2

Annualized Indices of the IEEE RTS . . . . . . . . . . . . . . . . . .

85

Table 6.3

Annual Indices of the Modified IEEE RTS . . . . . . . . . . . . . .

85

Table 6.4

Annualized Indices of the Modified IEEE RTS . . . . . . . . . . . .

85

Table 6.5

Reduction in the computational time of the IEEE RTS . . . . . . .

86

Table 6.6

Reduction in the computational time of the Modified IEEE RTS . .

86

Table 6.7

Bus Annual Indices for the IEEE RTS . . . . . . . . . . . . . . . . .

87

Table 6.8

Sensitivity Analysis of the LOLP index . . . . . . . . . . . . . . . .

88

Table 6.9

Sensitivity Analysis of the LOLF index . . . . . . . . . . . . . . . .

89

Table 6.10

Sensitivity Analysis of the EDNS index . . . . . . . . . . . . . . . .

89

Table 7.1

Annualized Indices of the IEEE Reliability Test System . . . . . . . 105

Table 7.2

Sensitivity Analysis of LOLP with Respect to Reliability Parameters
for the Base Case . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105

Table 7.3

Sensitivity Analysis of LOLP with Respect to Reliability Parameters
for the Reduced Emission Case . . . . . . . . . . . . . . . . . . . . . 106

Table 7.4

Sensitivity Analysis of LOLF with Respect to Reliability Parameters
for the Base Case . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106
xi

Table 7.5

Sensitivity Analysis of LOLF with Respect to Reliability Parameters
for the Reduced Emission Case . . . . . . . . . . . . . . . . . . . . . 107

Table 7.6

Sensitivity Analysis of EDNS with Respect to Reliability Parameters
for the Base Case . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107

Table 7.7

Sensitivity Analysis of EDNS with Respect to Reliability Parameters
for the Reduced Emission Case . . . . . . . . . . . . . . . . . . . . . 108

Table 7.8

Sensitivity Analysis of EDNS with Respect to Emission Allowance
Limits for the Base Case . . . . . . . . . . . . . . . . . . . . . . . . 108

Table 7.9

Sensitivity Analysis of EDNS with Respect to Emission Allowance
Limits for the Reduced Emission Case . . . . . . . . . . . . . . . . . 109

Table 7.10

Percentages of Increase in the Sensitivities of EDNS with Respect to
Emission Allowance Limits Due to Reduced the Emissions . . . . . . 109

Table 8.1

System Annualized Indices of the IEEE RTS . . . . . . . . . . . . . 117

Table 8.2

System Annualized Indices of the Modified IEEE RTS . . . . . . . . 117

Table 8.3

System Annualized Indices for Case I: Effect of Gradually Increasing
the Generation and Load Levels . . . . . . . . . . . . . . . . . . . . 118

Table 8.4

System Annualized Indices for Case II: Effect of Gradually Increasing
the Load Levels and Fixing the Generation . . . . . . . . . . . . . . 119

Table 8.5

Contributions of Voltage and Reactive Power Constraints . . . . . . 124

Table 8.6

Reliability Indices with Relaxing Voltage Limit Constraints . . . . . 125

Table 8.7

Sensitivity Analyses of the Severity Index with respect to the Voltage
and Reactive Power Constraints . . . . . . . . . . . . . . . . . . . . 126

Table 9.1

Sensitivity of the risk index with respect to minimum reactive power
limits for the modified IEEE RTS . . . . . . . . . . . . . . . . . . . 133

Table 9.2

Sensitivity of the risk index with respect to maximum reactive power
limits for the modified IEEE RTS . . . . . . . . . . . . . . . . . . . 133

Table 9.3

Sensitivity of the risk index with respect to under voltage limits for
the modified IEEE RTS system . . . . . . . . . . . . . . . . . . . . 134
xii

Table 9.4

Sensitivity of the risk index with respect to over voltage limits for the
modified IEEE RTS system . . . . . . . . . . . . . . . . . . . . . . . 134

xiii

LIST OF FIGURES

Figure 2.1

Line Current Linearization . . . . . . . . . . . . . . . . . . . . . . .

20

Figure 4.1

A bus connecting two transmission lines. If the total generation at
this bus and the sum of line capacities are less than the load, then
this state is a failure state . . . . . . . . . . . . . . . . . . . . . . .

43

A hypothetical system to demonstrate the application of the state
space reduction technique . . . . . . . . . . . . . . . . . . . . . . . .

45

Figure 4.3

State space representation with two conflicting objective functions .

54

Figure 4.4

State space representation of three dimension system (three components with multiple outcomes) . . . . . . . . . . . . . . . . . . . . .

58

Reliability indices profile in case the particles are trapped to the upper
corner of the state space . . . . . . . . . . . . . . . . . . . . . . . .

59

Reliability indices profile in case the particles are trapped to the lower
corner of the state space . . . . . . . . . . . . . . . . . . . . . . . .

60

Figure 4.2

Figure 4.5

Figure 4.6

Figure 4.7

Redirecting the swarm in case it enters the success or failure subspaces 61

Figure 5.1

Flow chart of the proposed method . . . . . . . . . . . . . . . . . .

68

Figure 6.1

Reliability indices profile of the IEEE RTS System . . . . . . . . . .

88

Figure 7.1

Flowchart of Estimating the Sensitivity Analysis . . . . . . . . . . . 102

Figure 8.1

System loss of load probability (LOLP) index profile against the
change in the system loading level . . . . . . . . . . . . . . . . . . . 120

Figure 8.2

Contributions of the voltage limit violations on the system loss of load
probability (LOLP) index against the change in the system loading
level . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121

xiv

Figure 8.3

Contributions of the reactive power limit violations on the system loss
of load probability (LOLP) index against the change in the system
loading level . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122

Figure 8.4

The profile of system loss of load probability (LOLP) index with
relaxing voltage limit constraints . . . . . . . . . . . . . . . . . . . . 125

xv

Chapter 1
Introduction
Due to the ongoing changes in generation portfolios and environmental concerns, reliability
evaluation of composite generation and transmission systems has received a great attention.
The task of composite power system reliability evaluation has become more complicated due
to the rapid increase in electrical demand and the liberalization of the electricity markets.
Several reliability indices for composite systems have been introduced such as loss of load
probability, loss of load frequency, loss of load expectation, expected energy not supplied,
expected demand not supplied, etc. Although these indices are of a great importance in
both planning and operational power system reliability evaluation, they lack the ability of
identifying the influence of each area or equipment on system reliability. Moreover, these
indices cannot determine the most-effective way to improve the system in order to keep it
reliable.
One of the promising methods that can identify the most vulnerable component, bus, or
area is performing sensitivity analysis of the possible load shedding with respect to component parameters and operating limits. The sensitivity analysis provides a measure for the
amount of load to be curtailed in response to the violation of the operating limits or component characteristics. Sensitivity analysis of the power system reliability indices is intended to
determine the effect of the capacity and reliability parameters such as failure and repair rates
of each component on system reliability. Also, it examines the effect of power quality such
as voltage and reactive power reserve on system behavior. Furthermore, sensitivity analysis
1

can be used and is being used in mitigating the unfolding events in case of cascading failures.
The task of performing sensitivity analysis is amply used in power system reliability planning and operational studies. However, the application of sensitivity analysis in the changing
operating environment such as the newly imposed emission constraints, penetration of the
renewable energy sources and cascading failures in power systems has not been given much
attention. Sensitivity analyses of the reliability indices can be conducted either analytically
such as using state space enumeration or by means of simulations. Therefore, performing
sensitivity over a large set of variables for large systems is computationally expensive.
In response to the growing need to harden power systems against unfolding events with
efficient computation tools, this thesis provides a framework for performing sensitivity analysis for the risk indices with respect system variables. To make the analyses feasible in terms
of computation time and burden, a state space reduction technique in conjunction with the
population-based intelligent search methods is proposed.
In addition to component parameters such as component availability/unavailability, failure rate and repair rate, the operating limits are also included in this work. The operating
limits studied in this work are emission caps and voltage and reactive power limits. The
work provided in this thesis is the first work in the literature that addresses the effects of the
emission limits on power system planning from reliability prospective. Although the effects
of voltage and reactive power limits on power system reliability have been considered, the
sensitivity of the load shedding with respect to these limits is the first time introduced in
this thesis.

2

1.1

Motivation

In response to the rapid growth of the global economics, large power grids are interconnected
in order to provide the required energy. As a result of these interconnections, the number of
failures also increases and may spread through these interconnections and cause cascading
failures. Also, due to the integration of renewable energy sources, voltage and reactive power
limits should be included in the planning and operating reliability assessments. Further, after
the introduction of the regulations and amendments to reduce the emissions from the power
generation sector, emission constraints should also be included in the reliability assessments.
The time and computation effort to evaluate the robustness of power systems against
the triggering events are of concern in both planning and operation prospective. Developing
tools that can reduce the computation burden will help in applying strategies to harden
power systems against catastrophic failures.

1.2

Challenges

The challenges that were addressed in this thesis are summarized as follows:
A. Determining the most vulnerable components in power systems that can lead to power
interruptions is a challenging task. Several indices have been introduced in the literature
of power system reliability evaluation but they lack the ability of identifying which
component that has the highest effect on system deficiency.
B. The task of power system reliability evaluation is becoming more computationally demanding due to the increase in the complexity of the power industry infrastructure.
The time and computation effort to evaluate power system reliability indices and their

3

sensitivities with respect to component parameters and system operating limits are of
concern in both planning and operation prospective. As the number of system components increases, the number of states that need to be tested to assess power system
reliability exponentially increases (2n in case of state space enumeration where n is the
number of system components).
C. Most of the existing power system reliability evaluation methods do not account for
the effects of the operating constraints such as voltage and reactive power limits. Also,
new regulations such as emission reduction regulations should be included in performing
power system reliability evaluation.
D. Mitigation of cascading failures in power systems has become of necessity especially after
the recent blackouts. Cascading failures usually develop in a sequence of events within
a time horizon of seconds to hours. If the potential causes of the cascading failures are
incorporated with the load shedding procedure, cascading failure could be mitigated.

1.3

Approaches

In addressing the aforementioned challenges, the following approaches were used:
A. Sensitivity analysis of the reliability indices are used to identify the most vulnerable
components in a power system. The sensitivity analysis can identify which component
that contributes in most and most sever load interruptions.
B. A state space reduction technique is proposed and used in reducing the computational
effort. This technique classifies the state space into success, failure and unclassified subspace by using heuristic devices. Controlled population-based intelligent search methods,
4

in particular, dynamically directed particle swarm optimization search method, are used
to search for the desired states in the unclassified subspace.
C. The effects of the voltage and reactive power constraints on power system availability
are considered by utilizing the AC power flow model and a linearized power flow model.
Also, the constraints of the emission limits are accommodated by piecewise linearizing
the quadratic heat rate equations.
D. The sensitivity of the risk of load curtailments with respect to voltage and reactive power
limits are used in mitigating cascading failures in power systems. After identifying the
buses that hit the limits, reactive power compensation and/or generation rescheduling
are suggested to prevent voltage collapse.

1.4

Contributions

The contributions of the presented work can be summarized as follows:
A. Developing and proposing new sensitivity analysis expressions that can relate the expected load curtailment to the component parameters as well as operating limits.
B. Proposing a new state space reduction technique that can classify the state space of
a given system into success, failure and unclassified subspaces. The validity and effectiveness of the proposed technique are tested through comparing the results with those
obtained using Monte Carlo simulation.
C. Proposing a dynamically directed particle swarm optimization search method to search
for failure states in the unclassified subspace or in any desired subspace.

5

D. Introducing a complementary concept that can be used with the population-based intelligent search methods to estimate the reliability of power systems without the need
of testing the entire state space.
E. Inclusion of emission constraints in power system reliability evaluation.
F. Proposing new indices to include the effects of voltage and reactive power limits on
power system reliability.
G. Developing a method to mitigate cascading failures in power systems through the use
of the sensitivity analysis concept.

1.5

Organization of the Thesis

The remainder of this thesis is organized as follows: Chapter 2 presents the power flow models that are used in this thesis and the incorporation of minimum load curtailments in an
optimization framework. Chapter 3 discusses the developments of the sensitivity analysis
expressions that are used to estimate the sensitivity of power reliability indices with respect
to component parameters and the operating limits. Chapter 4 introduces the proposed state
space reduction technique and the dynamically directed particle swarm optimization search
method. Chapter 5 describes the use of the shadow price concept (Lagrange multipliers)
in evaluating the sensitivity of power system reliability indices with respect to component
parameters. Also, chapter 5 provides benchmark results that can be used to test the effectiveness of the proposed approaches. Chapter 6 shows the applications of the proposed
approaches on performing sensitivity analysis of the reliability indices. Chapter 7 shows the
applications of the proposed methods with considering emission constraints. Chapter 8 con-

6

siders the effects of voltage and reactive power limits on power system reliability. Chapter 9
shows the applications of sensitivity analysis in hardening power systems against cascading
failures. Chapter 10 provides concluding remarks and possible future work.

7

Chapter 2
Power Flow Modeling and
Optimization Framework

2.1

Introduction

Power flow in power system networks plays a significant part of power system planning
and operating studies. Incorporation of a power flow model and redispatch optimization
problem with minimizing operation costs and load curtailment is commonly used in planning
and operation decisions. Two types of power flow models have been commonly used in the
literature in modeling power system networks which are AC power flow and DC power flow.
The DC power flow ignores the effect of line resistances and the reactive power and assumes
voltages magnitudes are 1 p.u. These assumptions make the DC power flow less accurate
because in some scenarios the voltage collapse due to voltage limits and insufficient reactive
power support can lead to under-voltage load shedding. The full AC power flow model is
the most accurate power flow model but it is computationally expensive especially for the
applications that require repetitive runs of power flow. In this thesis, in addition to the AC
and DC power flow models, we have used another power flow model which is a linearized
version of the AC power flow that closely represents power system networks. This model
was developed by the authors in [1]. This model can be considered as a modified version of
the DC power flow model or a linearized form of the full AC load flow model.
8

2.2

Power Flow Models

In evaluating the reliability indices of composite systems, an optimal power flow with an
objective of minimum load curtailment is performed. This section provides a brief overview
of the AC, DC and the linearized power flow models.

2.2.1

AC Power Flow Model

AC power flow model has been amply explained in the literature. This section provides a
brief introduction about the formulation of the AC power flow. The AC power flow equations
at bus k can be presented as follows (in the polar form),

Pk = Vk

Vm (Gkm cos δkm + Bkm sin δkm )

(2.1)

Vm (Gkm sin δkm − Bkm cos δkm )

(2.2)

m∈k

Qk = Vk
m∈k

where m ∈ k means the set of the buses connected to bus k, Pk and Qk are the real and
reactive power injected at bus k, Vk and Vm are the voltage magnitudes at buses k and m,
Gkm is the conductance between buses k and m, Bkm is the susceptance between buses k
and m and δkm is the angle difference between voltages of buses k and m (δkm = δk − δm ).
In solving the AC power flow equations, each bus has four variables which are Pk , Qk ,
Vk and δkm . There are three types of buses, load buses, generation buses and a slack or
reference bus. For load buses, which are sometimes called PQ buses, the Pk and Qk are
pre-specified for constant power load model. For generation buses, which are commonly
known as PV buses, Pk and Vk are pre-specified. For the slack bus (one slack bus), Pk and
δkm are pre-specified.
9

2.2.2

DC Power Flow model

DC power flow model has been amply explained in the literature. This section provides a
brief introduction about the formulation of the DC power flow. The DC power flow model
is an approximated version of the AC power flow model with the following assumptions:
1. Bus voltage magnitudes are assumed 1.0 p.u.,
2. Voltage angles, δkm , between buses are small such that cos(δkm ) ≈ 1 and sin(δkm ) ≈
δkm = δk − δm ,
3. Line resistances are much smaller than the line reactances so that line impedances can
be approximated as Zkm = jxkm (line susceptance Bkm = −j/xkm ),
4. Susceptances between buses and the ground are neglected.
The real power flow between buses can be expressed as,

Pkm =

|Vk ||Vm |
sin δkm
xkm

(2.3)

From the above assumptions, the real power flow between buses can be approximated as,

δ − δm
Pkm = k
xkm

(2.4)

Then, the real power at buses can expressed as follows,

Pkm .

Pk =
m∈k

10

(2.5)

2.2.3

Linearized Power Flow model

The development of the linearized power flow model starts with the well-known AC real and
reactive power injection equations at bus k (equations (2.1) and (2.2)). If the voltage at bus
k is assumed 1 p.u. for calculating the power injection equations, then equations (2.1) and
(2.2) can be approximated as follows,

Pk ≈

Vm (Gkm cos δkm + Bkm sin δkm )

(2.6)

Vm (Gkm sin δkm − Bkm cos δkm ) .

(2.7)

m∈k

Qk ≈
m∈k

It should be noted that, this assumption will not prevent the voltage at bus k from being
calculated, rather it approximates the power injection at bus k. In other words, this is only
an approximation that enables the linearization; it is not an assumption that the voltage
magnitude equals 1.0 p.u.
Differences in bus voltage angles usually are very small and equations (2.6) and (2.7)
can be further approximated by applying the cosine and sine rules. When δ is very small,
cos(δ) ≈ 1 and sin(δ) ≈ δ. Therefore,

Pk ≈

(Vm Gkm + Vm Bkm δkm )

(2.8)

(Vm Gkm δkm − Vm Bkm ) .

(2.9)

m∈k

Qk ≈
m∈k

With further manipulation, the above equations can be rewritten as shown in equation

11

(2.10).






 PK   −B

=
QK
−G





G 

−B

δK 

VK

(2.10)

where B and G are the susceptance and conductance sub matrices of the Ybus and B and
G are the susceptance and conductance sub matrices of the Ybus without including the
susceptances and conductances of the lines. The diagonal elements of B and G can be
expressed as,
Bkk = Bkk − bkk
Gkk = Gkk − gkk
where bkk is the sum of susceptances of the shunt elements that are connected at bus k and
gkk is the sum of conductances of the shunt elements that are connected at bus k.
This model has been applied on several known systems and the results show its validity
and it closely approximates the full AC power flow model. For further details about this
model, the readers are suggested to refer to [1].
The developed linearized power flow model does not include the losses in the transmission
lines. In large power systems with thousands of buses, real power losses can be substantial
which is not advisable to be neglected.
2

δ
If we keep the term of − km
2 of the Tylor expansion of the cosine function, while neglecting

other higher order terms, then we obtain,

Pk ≈

(Vm Gkm + Bkm δkm ) −
m∈Ωk

1
2

2
Gkm δkm

(2.11)

m∈Ωk

Note that the first part of (2.11) describes the real power flow that other buses receive

12

from bus k, while the second part indicates the real power loss in the power transition. Let
us define,
k =−
Ploss

1
2

2
Gkm δkm

(2.12)

m∈Ωk

Similarly, the reactive power losses can be expressed as,

Qkloss =

1
2

2
Bkm δkm

(2.13)

m∈Ωk

Equations (2.12) and (2.13) explicitly indicate active and reactive power losses are quadratic
k is a non-negative
forms of δkm . Gkm , the real part of Ykm , is always non-positive, hence Ploss

number. Ref. [2–4] described a piecewise method to linearize transmission line losses.
Note that (2.12) and (2.13) only include the losses produced by the series impedances
of transmission lines. Losses produced from the shunt part of the transmission lines have
already been considered in the original formulation expressed by (2.10).

2.3

Optimization Problem for Minimum Load Curtailment

In composite system reliability studies, power flow analyses are usually carried out in solving
optimization problems of minimum load curtailment. Optimization with minimum load
curtailment has been extensively used in calculating the reliability indices of composite power
systems. If, under any scenario, the curtailment is unavoidable, the optimization problem
tries to minimize the amount of load to be shed.

13

2.3.1

Network Modeling using Nonlinear Programming and the
AC Power Flow Model

This section describes the formulation and incorporation of the objective function of minimum load curtailment in the nonlinear programming problem with the AC power flow model.
This objective function is subjected to equality and inequality constraints of the power system operating limits. The equality constraints include the power balance at each bus and the
inequality constraints are the capacity limits of generating units, power carrying capabilities
of transmission lines, voltage limits at the nodes and reactive power capability limits. The
minimization problem is formulated as follows [5],




Nb

Loss of Load = min 

Ci 

(2.14)

i=1

Subject to
P (V, δ) − PD + C = 0
Q(V, δ) − QD = 0
min ≤ P (V, δ) ≤ P max
PG
G
max
Qmin
G ≤ Q(V, δ) ≤ QG

V min

≤V ≤

(2.15)

V max

S(V, δ) ≤ S max
0 ≤ C ≤ PD
δ unrestricted.
In (2.14) and (2.15), Ci is the load curtailment at bus i, C is the vector of load curtailments (Nd × 1), V is the vector of bus voltage magnitudes (Nb × 1), δ is the vector of

14

bus voltage angles (Nb × 1), PD and QD are the vectors of real and reactive power loads
min , P max , Qmin and Qmax are the vectors of real and reactive power limits of
(Nd × 1), PG
G
G
G

the generators Ng × 1 , V max and V min are the vectors of maximum and minimum allowed
voltage magnitudes (Nb × 1), S(V, δ) is the vector of power flows in the lines (Nt × 1), S max
is the vector of power rating limits of the transmission lines (Nt × 1) and P (V, δ) and Q(V, δ)
are the vectors of real and reactive power injections (Nb × 1). Also, Nb is the number of
buses, Nd is the number of load buses, Nt is the number of transmission lines and Ng is the
number of generators.
In the standard minimization problem given by (2.14) and (2.15), all generation and
network constraints have been taken into consideration. Also, it has been assumed that one
of the bus angles is zero in the constraints (2.15) to work as a reference bus.

2.3.2

Network Modeling using Linear Programming and the DC
Power Flow Model

The DC power flow model has been widely utilized in reliability assessment of power systems
due to its simplicity of formulation and implementation [6–8,10–17]. Moreover, the DC power
flow model has the advantage of being suitable for studies that require extensive computational burden such as composite system reliability and security assessment. Using DC power
flow, there are three main constraints which are power balance equation, generation capacity
limits and transmission lines power carrying capabilities. The linear programming model of
the network with DC power flow is given in equation (2.16) which is adapted from [7],


Nb

Loss of Load = min 

Ci 
i=1

15


(2.16)

subject to

Bδ + G + C = D
G ≤ Gmax
C ≤ D
bAδ ≤ Ffmax

(2.17)

−bAδ ≤ Frmax
G, C ≥ 0
δ

unrestricted

where Nb is number of buses, Nt is number of transmission lines, B is the augmented node
susceptance matrix (Nb × Nb ), b is the transmission line susceptance matrix (Nt × Nt ), A is
the element-node incidence matrix (Nt × Nb ), δ is the vector of node voltage angles (Nb × 1),
C is the vector of bus load curtailments (Nb × 1), D is the vector of bus demand (Nb × 1),
Gmax is the vector of maximum available generation (Nb × 1), Ffmax is the vector of forward
flow capacities of lines (Nt × 1), Frmax is the vector of reverse flow capacities of lines (Nt × 1),
and G is the solution vector of the generation at buses (Nb × 1).
In the standard minimization problem given by (2.16) and (2.17), all generation and
network constraints have been taken into consideration. Moreover, it has been assumed that
one of the bus angles is zero in the constraints (2.17) to work as a reference bus.

16

2.3.3

Network Modeling using Linear Programming and the Linearized Power Flow Model

In this section, the linear programming optimization problem is incorporated with the linearized power flow to minimize load curtailments. Using the linearized power flow model,
in addition to the power balance equation, generation capacity limits and transmission lines
power carrying capabilities, more constraints can be added. These constraints are voltage
allowable limits and generation reactive power limits. This section provides explanations for
the incorporation of theses constraints to the optimization problem.

2.3.3.1

Power Balance Constraints

Power balance equation is an equality equation represents the sum of the complex power at
a bus. The balance equation can be derived from equation (2.10) as follows:

B δ − GV + PG = PD
(2.18)
G δ + BV + QG = QD

where δ is the vector of bus voltage angles, V is the vector of bus voltage magnitudes, PG
and QG are the vectors of the real and reactive power of the generators respectively and PD
and QD are the vectors of the real and reactive power of the bus loads respectively.
Since the problem is to minimize the curtailment, equation (2.18) is augmented by fictitious generators that are equivalent to the required curtailment. Therefore, equation (2.18)

17

becomes,
B δ − GV + PG + PC = PD
(2.19)
G δ + BV + QG + QC = QD
where PC and QC are the vectors of the real and reactive power of the load curtailments
respectively.

2.3.3.2

Real and Reactive Power Constraints of the Generators

Real power limits of the generators are bounded by maximum available capacity and minimum power where the latter is considered zero. Reactive power limits are the maximum
reactive power provided by a generator and the minimum reactive power assigned for each
generator. Therefore, the constraints of real and reactive power of the generators can be
expressed as,
max
0 ≤ PG ≤ PG

Qmin
G

≤ QG ≤

(2.20)

Qmax
G

max is the maximum available capacity for each generator, Qmin is minimum reactive
where PG
G

power can be absorbed by a generator and Qmax
is the maximum reactive power can be
G
produced by a generator.

2.3.3.3

Voltage Limits Constraints

Voltage constraints are limited according to the allowed voltage fluctuations. The maximum
voltage limit, V max , and the minimum voltage limit, V min , which are assumed throughout
this work as 1.05 p.u. and 0.95 p.u. respectively, are expressed as follows,

V min ≤ V ≤ V max .
18

(2.21)

2.3.3.4

Line Current Limit

The current flowing through a line connecting buses k and m can be calculated from the
voltage difference and branch impedance as follows,

V − Vm
Ikm = k
Zkm

(2.22)

Equation (2.22) can be manipulated to produce a linearized expression for the line current.
This approximation is based on the assumption that the differences between bus voltage
angles are very small such that cos δkm ≈ 1 and sin δkm ≈ δkm and line resistances is much
smaller than line reactances, i.e., Rik

| Ikm | ≈

Xik . Therefore,

(δk − δm )2 + (Vk − Vm )2 / | Zkm |
(2.23)

=

p

Ikm

2

q

+ Ikm

2

where
p

(δk − δm )
| Zkm |

q

(Vk − Vm )
| Zkm |

Ikm =
Ikm =
p

q

where Ikm and Ikm are the real and imaginary components of the line current respectively.
Therefore, the line current can be expressed by two linear components but the magnitude of
these two components is non-linear.
With a constant maximum line current, the locus of the real and imaginary components
forms a circle as shown in Figure 2.1. To linearize the current, the ratio of the reactive power
to the real power flowing through a line needs to be known in advance which is not possible in
19

I ikq

I max

I max sin( m )

I max sin( 3 )

I max sin( 2 )

I max sin(1 )

m

I

l

max

3
2
1

I max

cos( m )

I

max

cos( 3 )

I

max

cos( 2 )

I ikp
I

max

cos(1 )

Figure 2.1: Line Current Linearization
solving linear programming problem. However, this ratio also can be linearly approximated
around the operating point. This can be done by running unconstrained linear programming
problem for several load scenarios to evaluate the average and standard deviation of these
ratios. From the average and standard deviation of these ratios, the two components of the
line current can be approximated by n linear segments according to the required accuracy.
Therefore, the curve around the average ratio can be linearized as mentioned above. The
angle θl in Figure 2.1 is the tangent inverse of the ratio that corresponds to segment l.
From the analysis of Figure 2.1 the linear equations of the linear segments can be expressed as,
p

q

max = 0
Ikml cos θl + Ikml sin θl − Ikm

(2.24)

where l = 1, 2, ..., n
Therefore, by choosing the number of linear segments, n, the current limit constraints will
be increased with order of n times. The current limit constraints in the linear programming

20

problem can be expressed as:

bA δ + bA V ≤ Ifmax

(2.25)

−bA δ − bA V ≤ Irmax
where b is a diagonal matrix of the transmission line admittances (Nt × Nt ), A is the elementnode incidence matrix (Nt × Na ) and

bA = bA cos θl
bA = bA sin θl
As mentioned above, as the number of linear segments increases, the number of line
constraints increases accordingly which it may slow down the computation speed. However,
in most practical power systems it is not economical to transmit reactive power through
transmission lines; therefore, the ratio of reactive power flowing in a line to the real power
is usually small. Hence, linear segments can be limited to one or two segments around the
maximum real power limits.

2.3.3.5

Linear Programming Formulation

The Linear Programming is used here to minimize load curtailments similar to the formulation of the DC power flow model [7] with incorporating voltage and reactive power
constraints. From the above constraints, the minimization problem can be formulated as,


Na

Loss of Load = min 

Ci 
i=1

21


(2.26)

Subject to:

B δ − GV + PG + PC = PD
G δ + BV + QG + QC = QD
max
0 ≤ PG ≤ PG
max
Qmin
G ≤ QG ≤ QG

0 ≤ P C ≤ PD

(2.27)

0 ≤ QC ≤ QD
V min ≤ V ≤ V max
bA δ + bA V ≤ Ifmax
−bA δ − bA V ≤ Irmax
δ

unrestricted

where Na is number of buses, Nt is number of transmission lines, B , B, G and G are
as given in section II with dimensions of (Na × Na ), δ is the vector of node voltage angles
(Na × 1), V is the vector of bus voltage magnitudes (Na × 1), PG is the vector of real power
generation (Na × 1), QG is the vector of reactive power generation (Na × 1), PC is the
vector of real load curtailments (Na × 1), QC is the vector of reactive load curtailments
(Na × 1), PD is the vector of bus real loads (Na × 1), QD is the vector of bus reactive loads
max is the vector of maximum available real power generation (N × 1), Qmax
(Na × 1), PG
a
G

is the vector of maximum available reactive power generation (Na × 1), Qmin
is the vector
G
of minimum available reactive power generation (Na × 1), V max is the vector of maximum
allowable voltages (Na × 1), V min is the vector of minimum allowable voltages (Na × 1),

22

max
Ffmax
−real is the vector of real forward flow capacities of lines (Nt × 1), Ff −imaginary is the
max is the vector of real
vector of imaginary forward flow capacities of lines (Nt × 1), Fr−real
max
reverse flow capacities of lines (Nt × 1) and Fr−imaginary
is the vector of imaginary reverse

flow capacities of lines (Nt × 1).
In (2.26) and (2.27) all generation availability and network constrains have been taken
into considerations. Also, in order to get a feasible solution for the standard problem, it has
been assumed that one of the bus angles is zero in the constraints of equation (2.26).
Formulation of the optimization problem in this manner may produce a solution with
different curtailments for the active and reactive powers. In other words, if a load is assumed
to be curtailed in case of insufficient real power generation, this formulation may shed the
real part of the load and leaves the reactive part. To insure that the curtailment solves for
real and reactive powers, an additional constraint has to be added. This constraint can be
expressed as,
PC − (PD /QD ) QC = 0

(2.28)

The factor (PD /QD ) is the ratio between the real load to the reactive load at each bus.
Therefore, this constraint will insure that both the active and reactive power will be curtailed
if there is no enough real power as well as if there is no enough reactive power.

2.3.4

Other Constraints

Other constraints can be incorporated into the optimization problem such as emission constraints, transient stability constraints, etc. In this thesis, the emission caps and limits have
been included in performing sensitivity analysis for the reliability indices.
The emission constraints are used in this work to determine the point at which the

23

action of buying or selling the emission allowances can be decided. Also, the shadow prices
of emission constraints at each power plant are used in estimating the sensitivity of the
objective function with respect to the emission limits. The emission constraints are derived
as follows,
max
ECO2j ≤ ECO2j

∀j.

(2.29)

In (2.29), ECO2j is the total emission of CO2 in tonne/t at the power plant j and
max is the cap on the CO emission at the power plant j. More details about the
ECO2j
2

incorporation of the emission constraints are given in chapter 7.

24

Chapter 3
Power System Reliability Indices and
Sensitivity Analysis

3.1

Introduction

Due to the ongoing changes in generation portfolios and environmental concerns several
reliability indices for composite systems have been introduced. Even though these indices are
of a great importance in both planning and operational power system reliability evaluation,
they lack the ability of identifying the influence of each area or equipment on the system
reliability. Moreover, these indices do not give a full picture of the most-effective way to
revamp the system in order to keep it reliable. Due to these reasons, significant efforts
have been devoted to this area in recent years to evaluate what the system’s reliability
justifications are and where the best location to invest is. This chapter provides explanations
for the reliability indices and the developments of the sensitivity of these indices with respect
to the component parameters.

3.2

Calculation of the Reliability Indices

In this work, we have evaluated the well know composite power system reliability indices,
namely loss of load probability (LOLP), expected demand not supplied (EDNS), loss of

25

energy expectation (LOEE), loss of load expectation (LOLE), loss of load frequency (LOLF)
and expected duration of load curtailment (EDLC).

3.2.1

Calculation of System Probability Indices

Failure probability indices evaluate the probability of failure of the system to meet the
demand. Another index was introduced which is the probability of load availability. This
index represents the probability of the system to succeed to meet the demand. Through the
search process, if the state under consideration is a failure state, the probability of this state
is added to the failure probability index, q; otherwise, the probability of this state is added
to the success probability index, p. The failure probability index and the success probability
index are bounded between 0 and 1.0. If the algorithm were to visit the entire state space,
the sum of q and p would equal to 1.0. In other words, q and 1 − p would reach each other
and would evaluate the same index, loss of load probability.
The probability of system failure to meet the demand is given by
nf

P xi : xi ∈ X f

q=

(3.1)

i=1

where X is the set of all states, Xf is the set of failure states (Xf ⊂ X) and nf is the
number of failure states.
The probability of system success to meet the demand is given by
ns

P {xi : xi ∈ Xs }

p=
i=1

where Xs is the set of success states (Xs ⊂ X) and ns is the number of success states.

26

(3.2)

The estimated q, (ˆ
q ), index can be calculated based on the complementary concept (the
complementary concept is introduced in chapter 4) as follows,

qˆ =

q
q+p

(3.3)

The estimated loss of load expectation can be calculated in the same manner. Let Lq
ˆ q ), index can be calculated
denote loss of load expectation index. The estimated Lq , (L
directly from the qˆ index as follows,

ˆq =
L

q
×T
q+p

(3.4)

where T is the period of study in hours.

3.2.2

Calculation of System Energy Indices

The well-known reliability energy indices are expected demand not supplied and expected
energy not supplied. Let dn denote the expected demand not supplied index and en denote
the expected energy not supplied index. In this work, another index was introduced which
is expected demand supplied and denoted as ds . This index estimates the amount of power
supplied to system loads. For every tested state, if the state is a success state, the product of
the probability of this state and the peak load is added to the ds index. On the other hand,
if the state is a failure state, the product of the probability of this state and the amount of
load curtailment is added to the dn index and the product of the probability of this state and
the amount of supplied load (the supplied load is the peak load minus the curtailed load)
is added to the ds index. The expected demand not supplied and the expected demand

27

supplied indices are bounded between 0 and the peak load. If the algorithm were to visit the
entire state space, the sum of dn and ds would equal to the peak load. In other words, dn
and (Peak Load − ds ) would reach each other and would evaluate the same index, expected
demand not supplied.
The expected demand not supplied is given by,
nf

P xi : xi ∈ Xf × Lc xi : xi ∈ Xf

dn =

(3.5)

i=1

where Lc is amount of load curtailment of state xi .
The expected demand supplied is given by,
nf

P xi : xi ∈ Xf × Ls {xi : xi ∈ X}

ds =

(3.6)

i=1

where Ls is amount of load supplied of state xi and can be expressed as,




Peak Load,

x i ∈ Xs
Ls {xi : xi ∈ X} =



Peak Load − Lc, xi ∈ Xf

(3.7)

The estimated value of dn , (dˆn ), index can be evaluated as follows,

dˆn =

dn
× Peak Load
dn + ds

(3.8)

The en index can be calculated from the estimated dn index. The estimated value of en ,

28

(ˆ
en ), can be calculated as follows,

eˆn =

3.2.3

dn
× Peak Load × T .
dn + ds

(3.9)

Calculation of System Frequency and Duration Indices

Calculation of frequency and duration indices is generally more difficult than calculation
of probability and energy indices. As it has been mentioned, the probability indices are
bounded between 0 and 1 and energy indices are bounded between 0 and the peak load.
On the other hand, frequency and duration indices do not have this property. However, as
a result of the frequency balance property, Ff = Fs , failure frequency can be calculated
either from the success states that have downward transitions crossing the boundary between
the success and failure states, Fs , or from the failure states that have upward transitions
crossing the boundary between the success and failure states, Ff . If the algorithm were to
visit the entire state space, these two frequencies would converge to the same value. Since
the search method is not intended to visit the entire state space (as shown in chapter 4) and
the frequency is related to the probabilities of the states, frequency of failure index can be
estimated by taking the average of these frequencies, Fs and Ff , over the visited subspace.
If we denote the frequency of the failure index as F , then the estimated value of F , (Fˆ ), can
be expressed as follows,
Fˆ =

Ff + Fs
2 (q + p)

(3.10)

After estimating the probability of failure index, qˆ, and the frequency of failure index,

29

Fˆ , it is trivial to estimate the expected duration of load curtailment, Tˆc , index as follows,

2q
qˆ
=
Tˆc =
Ff + Fs
Fˆ

(3.11)

The approach used in this work to calculate frequency and duration indices was adapted
from [8, 18, 19]. The work presented in [8, 18] was based on concepts described in details
in [20, 21]. We will not reproduce the rigorous derivation here; rather, the expressions will
be presented for convenience.
Every failure state comprises of functional components and failed components. Failure
states can transit to higher states by upward transitions, and the resulted states can be
success or failure states. These transitions are said to be crossed the boundary between
the success and failure states if the failure states transit to success states. Also, the same
procedure can be applied on the success states that transit downwards to failure states.
The frequency of encountering success states (or failure states) is the sum of the individual
frequencies associated with those transitions of the failure states (or success states) which
cross the boundary.
The sum of the frequencies of the upward transitions from the failed states can be expressed as,

nf

F (+) =





 P xi : xi ∈ X f
i=1

µj 

(3.12)

j∈Fi

where Fi is the set of failed components in state i, and µj is the repair rate of component j.
The sum of the frequencies of the downward transitions from the failed states can be
expressed as,

nf

F (−) =




P xi : xi ∈ Xf

i=1

λk 
k∈Si

30

(3.13)

where Si is the set of functional components in state i, and λk is the failure rate of component
k.
Some of F (+) transitions will cross the boundary and the others will not. Also, by the
assumption that the system is coherent, none of F (−) will cross the boundary. As a result
of frequency balance property, all those transitions included in F (+) which do not cross the
boundary are balanced by corresponding transitions included in F (−) . Consequently, the
frequency of system success can be calculated as,

Fs = F (+) − F (−)

(3.14)

In the steady state, the frequency of system success and frequency of system failure will
balance, that is Ff = Fs . Therefore, the frequency of system failure can be calculated from
failure states as follows,
nf







P xi : xi ∈ Xf 

Ff =
i=1

µj −
j∈Fi

λk 

(3.15)

k∈Si

Following the similar approach, frequency of system failure can be calculated from the
success states as follows,
ns







P {xi : xi ∈ Xs } 

Fs =
i=1

λk −
k∈Si

µj 
j∈Fi

where Xs is the set of success states and ns is the number of success states.

31

(3.16)

3.2.4

Calculation of Bus Indices

Calculation of bus indices utilizing optimization problems produces multiple solutions and
hence bus indices will not be unique [22]. Depending on the manner in which the program
scans the vertices of the feasible polytopes, it can favor curtailing power at certain buses depending on how the buses are numbered. In the literature, load priority philosophy technique
is usually used to overcome this problem. Such philosophies can be the priority of each load
or part of each load. Also, load can be curtailed according to closeness to the fault location.
In this method, the loads are divided into several parts with different weighting factors. Each
part is assumed to represent a percentage of the total load. However, if the amount of load
curtailment at each bus is less than the load priority level, the multi-optimum problem will
occur. This bias can be removed by dynamic numbering of the buses, i.e., altering the bus
numbers after the occurrence of each event in the simulation.
In this work we have combined both approaches by adapting the load priority philosophy
technique and dynamic numbering of the buses. Loads were divided into three parts with
different weighting factors for each part, e.g., w1 , w2 and w3 respectively. The first two parts
were assumed to represent 25% of the total load and the third part was assumed to represent
50% of the total load. The weighting factor of the first part is assumed to be less than that
of the second part and the weighting factor of the second part is assumed to be less than
that of the third part or in other words, w1 < w2 < w3 .

3.2.5

Modeling of System Components

There several methods to model system components in power system reliability studies. This
section presents the modeling of generators, loads and transmission lines.

32

3.2.5.1

Modeling of Available Generation

Most buses have several generators which may be similar or different. The unit addition
algorithm [23] is used to construct a discrete probability distribution function for each bus
which is known as capacity outage probability and frequency table, COPAFT. This table is
constructed based on the capacity states and forced outage rates of units at each bus.

3.2.5.2

Modeling of System Load

Loads at the buses is modeled based on the cluster load model technique [23–26]. From the
chronological loads, clusters are constructed according to the load level and its probability.
These clusters are used for each load bus as a percentage of the peak load of the given bus.

3.2.5.3

Modeling of Transmission Lines

A discrete probability density function is constructed for every transmission line. If a line is
tripped for some system state, the line is removed from the bus admittance matrix and its
capacity is set to zero.

3.2.6

A Stopping Criterion

A convergence criterion should be applied to stop the algorithm if there is not much change in
the reliability indices. In power system reliability analysis using Monte Carlo simulation, it
was found to be that energy indices are the slowest indices from convergence view point [22].
In this work, we have applied the stopping criterion on the EDNS index.
The stopping criterion considering the EDNS index can be expressed as,

σ=

V ar(EDNS)
E [EDNS]
33

(3.17)

where E [.] is the expectation operator and V ar (.) is the variance function.
After few iterations, the amount of the change in σ is calculated, if this amount is less
than or equal to the specified tolerance, the algorithm is to be terminated; otherwise, the
simulation will continue. This stopping criterion can be used for both the Monte Carlo
simulation method and the population-based intelligent search methods.

3.3

Concept of the Sensitivity Analysis

Lagrange multipliers based method is proposed for performing sensitivity analysis in composite system reliability. Lagrange multipliers are defined as the sensitivity of the value of
the objective function to the change in the right hand side of the linear/non-linear programming problem constraints [27–30]. Lagrange multipliers are used in evaluating the sensitivity
analysis of the reliability indices with respect to component parameters. The sensitivity of
the reliability indices with respect to the component parameters are used as a decision making tool in identifying the system’s component or area that has the highest impact on system
reliability; and which area or component need to be reinforced to enhance the overall system reliability. Lagrange multipliers can be determined by solving for the dual solution of
the optimization problem with an objective function of minimizing buses load curtailments.
The proposed method relies on component availability data, generation capacity, transmission line capability and load scenarios. In most of the practical applications reported in the
literature, the sensitivity evaluation techniques include approximations in generating capacity model, load model and the evaluation techniques. The proposed method uses Monte
Carlo simulation and population-based intelligent search methods that do not necessitate
such approximations.

34

3.4

Sensitivity Analysis of the Reliability Indices

Sensitivity analyses of reliability indices reported in the literature have been based on calculating the amount of change of these indices with respect to component parameters such as
availability/unavailability, capacity, failure rate and repair rate. These analyses have been
conducted either by examining every single state of the system or by Monte Carlo simulation. The analyses have been performed by casting the dispatch operation as an optimization
problem, through minimizing load curtailment or maximizing load supplying capability. The
optimization has been constrained by generator and transmission line capacities. The LOLP
and LOLF indices from their definitions are based on reliability parameters and cannot be
directly related to the operating limits. On the other hand, the EDNS index is based on
failure rates, repair rates and unit capacities. In this section, we provide a brief description
of the sensitivity analysis of the reliability indices and the detailed derivations are given in
the appendix.

3.4.1

Sensitivity Analysis of the LOLP and LOLF Indices with
Respect to Component Parameters

Sensitivity studies of LOLP and LOLF to component reliability have been amply described in
the literature. Sensitivity analysis can be conducted analytically by enumerating all system
states or by simulation. Refs. [16, 31–33] provide relationships that are suitable for use in
state space enumeration.
The sensitivity of LOLP with respect to unavailability ui of component i can be calculated
as follows,
If (x)P (x)[(1/ui ) − Si /(ai ui )]

∂LOLP/∂ui =
x∈X

35

(3.18)

where X is the set of all states, x is the state of the system, P (x) is the probability of
occurrence of state x, ui is the probability of failure of component i and ai is the probability
of success of component i. Si is the state indicator of component i, i.e., Si = 0 if component
i is in the down state (failure state) and Si = 1 if component i is in the up state (success
state) and If (x) is the system state indicator function which can be expressed as follows,

If (x) =



 1 if x is failure state

(3.19)


 0 if x is success state
The sensitivity of LOLP with respect to component failure rate λi is given by

If (x)P (x)[ai /λi − Si /λi ]

∂LOLP/∂λi =

(3.20)

x∈X

The sensitivity of LOLP with respect to component repair rate µi is given by

∂LOLP/∂µi =

If (x)P (x)[−ai /µi + Si /µi ]

(3.21)

x∈X

The sensitivity of LOLF with respect to unavailability ui of component i can be calculated
as follows.

−If (x)P (x)Si µi /a2i

∂LOLF/∂ui =
x∈X

+F (x)P (x) ((1/ui ) − Si /(ai ui ))]

(3.22)

where F (x) is the sum of the repair rates of a failure state x that crosses the boundary and

36

can be expressed as:
m

λin
i (x)

F (x) = If (x)
i=1

where m is the number of components, and

λin
i (x) = (1 − Si ) µi − Si µi ui /ai

The sensitivity of LOLF with respect to component failure rate λi is given by

−If (x)P (x)Si

∂LOLF/∂λi =
x∈X

+F (x)P (x) (ai /λi − Si /λi )]

(3.23)

The sensitivity of LOLF with respect to component repair rate µi is given by

If (x) (1 − Sk ) P (x)

∂LOLF/∂µi =
x∈X

+F (x)P (x) (−ai /µi + Si /µi )]

3.4.2

(3.24)

Sensitivity Analysis of the EDNS index with Respect to
Component Capacities

The expected demand not supplied (EDNS) is an important index because it inherently
reflects the severity of the events. The sensitivities of this index with respect to ui , λi and
µi are the same as for the LOLP index except that they are multiplied by the amount of
load curtailments. The sensitivity analysis of the EDNS with respect component capacities

37

was derived as follows [16, 32, 33],

∂EDN S/∂Ci =

If (x) P (x) ∂Lc (x) /∂Ci

(3.25)

x∈X

where Ci is the capacity of component i and Lc (x) is the total load curtailment when the
system is at state x.
The derivative of the total load curtailment with respect to component capacity can be
expressed as


 πg,i if component i is a generator
∂Lc (x) /∂Ci =

 π
t,i if component i is a circuit element
where πg,i and πt,ij are the Lagrange multipliers or shadow prices of generation capacity
constraints and transmission lines carrying capability constraints respectively. πg,i can be
calculated directly from the optimization problem. However, πt,ij depends on circuit parameters which are circuit capacity and susceptance. These two parameters are dependent
variables and cannot be treated separately. Pereira and Pinto [34], combined the effect of
circuit capacity and susceptance on circuit sensitivity and developed the following expression.

πt,ij = πd,i − πd,j

θj − θi

(3.26)

where πd,i and πd,j are Lagrange multipliers or shadow prices of load constraints of buses i
and j respectively and θj and θi are voltage angles of buses j and i respectively.

38

Chapter 4
State Space Reduction and
Population-based Intelligent Search
Methods

4.1

Introduction

Composite system reliability evaluation aims at determining the reliability of the given power
system taking into consideration both transmission and generation systems. Numerous techniques have been proposed in the literature to assess composite system reliability. In this
context, analytical methods [7, 35, 36] and Monte Carlo simulation [37] have been used for
composite system reliability evaluation. In evaluating the reliability indices of composite
power systems, a power flow or optimal power flow with an objective of minimum load
curtailment is usually required to test whether the state under consideration is a failure or
success state. Performing optimal power flow for huge number of scenarios can be computationally demanding. Consequently, numerous techniques have been proposed in the literature
to reduce the computational burden and the time spent in evaluating the reliability indices
of composite systems.
This chapter introduces a heuristic technique, which classifies the search space into fail-

39

ure, success, and unclassified subspaces, and introduces an algorithm, which is developed
based on a directed binary particle swarm optimization search method, to search for failure
states in the unclassified subspace. The proposed heuristic technique is developed based on
calculating the maximum capacity flow of the transmission lines and the available generation.
A key element in using particle swarm optimization to search for failure states in the unclassified subspace lies in selecting the weighting factors associated with the objective function.
Appropriate values of these weighting factors should be carefully selected in order to prevent
the swarm from being trapped to one corner of the state space. The work presented in this
chapter proposes an intelligent particle swarm optimization based search method to adjust
these weighting factors in a dynamic fashion. The effectiveness of the proposed method is
demonstrated on several test systems which are given in the next chapters. The results show
not only that the reliability indices obtained using the proposed method correspond closely
with those obtained using Monte Carlo simulation, but also the computational burden can
be significantly reduced using the proposed method.
The use of binary particle swarm optimization in power system reliability is found to be an
efficient tool in evaluating the reliability indices. However, the task of choosing the weighting
factors associated with the objective function tends to be tedious and time consuming as
trial and error approaches are often times being used to prevent the particles from being
flying in one direction. Consequently, the development of an automated approach to adjust
these weighting factors could save some effort in solving the task of reliability evaluation
efficiently.
This work proposes a technique to search for failure states in the unclassified subspace,
which is developed based on an intelligent bounded and directed binary particle swarm optimization. This technique dynamically adjusts the weighting factors associated with the
40

objective function of the binary particle swarm optimization so that the swarm has always
forced to fly on the unclassified subspace. A dynamically directed binary particle swarm
optimization search method was developed by the authors in [11] to search for failure states
in the entire state space. In this work, instead of using the dynamically directed binary
particle swarm optimization search method to search for failure states in the entire state
space, we use this technique to search for failure states in the unclassified subspace. Another intelligence factor is added to the binary particle swarm optimization that reverses the
direction of the particles to search in the unclassified subspace if the algorithm discovers that
some particles have entered the success and/or failure subspaces. The proposed algorithm
tracks the behavior of the particles and then adjusts the objective function coefficients correspondingly. The behavior of the swarm can be evaluated by examining the probabilities
of the visited states.

4.2

State Space Reduction Technique

Searching for failure states in the entire state space is time consuming; and hence, a reduction technique is required to reduce the computational time. Several methods have been
introduced in the power system literature to reduce the search space and the computational
effort. Singh and Mitra [6] have introduced the concept of the state space pruning. In [6],
an arbitrary set of coherent acceptable subspaces were pruned out from the state space using the concept of partitioning vectors. Then, Monte Carlo simulation was performed for
the remaining subspaces. The number of pruned subspaces is system dependent and it is
a tradeoff between the required time to prune the acceptable subspaces and the required
time to perform Monte Carlo simulation for the rest of the subspace. The application of

41

the state space pruning to calculate the frequency and duration indices of composite power
systems was carried out by Mitra and Singh in [8]. State space partitioning has been introduced in [9]. Some other researchers have developed the concept of state space pruning using
population-based intelligent search methods [38–44]. The genetic algorithm and the modified genetic algorithm have been used in [38–40] to prune the state space and Monte Carlo
simulation has been then used for the remaining part of the state space. Binary particle
swarm optimization technique has been utilized in [42, 44] to prune the state space. Green
et al. [43] have applied the artificial immune systems optimization technique in pruning the
state space. A comparative study of using population-based intelligent search methods, in
particular, genetic algorithms, repulsive binary particle swarm optimization, and binary ant
colony optimization has been held in [41]. It is worth to point out here that in the work presented in [38–44], population-based intelligent search methods have been utilized as heuristic
techniques and the selection of the acceleration factors was a key element in guaranteeing
the best performance of these methods.
In this section, a heuristic technique is introduced to prune the state space; and thereby
reduces the computational effort. This heuristic technique is developed based on calculating
the line flow capacity limits of the transmission lines. This technique models power system networks based on network configuration, available generation, loading conditions and
transmission line availabilities, and power carrying capabilities.
By examining the optimization problems presented in chapter 2, it can be noted that for
any sampled state if at least one of the following two conditions is satisfied, we can conclude
that this state is a failure state: (1) total generation is less than the total load and (2) sum
of the capacity of the lines connected to a bus is less than the loading of that bus as shown
in Figure 4.1. Given these conditions, we can construct a heuristic technique that reduces
42

the search space significantly. The procedures can be described as follows,
Bus K
Line i-k

Line k-j

Generator

Load

Figure 4.1: A bus connecting two transmission lines. If the total generation at this bus and
the sum of line capacities are less than the load, then this state is a failure state

1. Subtract the loads connected to the generation buses from the sampled amount of
generation at these buses.
2. Construct a line capacity capability flow matrix (LCM) for the given sampled configuration. This matrix can be constructed in the same manner as constructing a “Ybus ”
from the branch-node incidence matrix (A) except that instead of using “−1” for the
branch that is assumed to enter a node, we use “1”. In this case, the incidence matrix
is denoted A. Further, instead of using the diagonal susceptance matrix, we use a
diagonal capacity matrix (K) so that the diagonal entries of this matrix are the line
capacities. The LCM can be expressed as follows,

LCM = AT K A

(4.1)

3. At any bus, if the total generation is larger than the corresponding diagonal element of
the LCM, adjust the generation at that bus downward to the corresponding diagonal
element.
4. For a state x, if the value of any of the diagonal elements of the LCM is less than the
43

absolute value of the power injection at the corresponding bus or the total generation
is less than the total demand, this state is guaranteed to be a failure state.
After reading system data, system states are sampled. For every sampled state, the above
search space reduction technique is applied. Then, failure states are passed to the optimization problem and the reliability indices are updated. These processes are repeated for every
sampled system state.
As an example, consider the hypothetical system shown in Figure 4.2. This system
consists of five buses and 6 transmission lines. Line capacities are shown on the figure and
denoted as C1 , C2 , ... and C6 .
The LCM is shown in (4.2) and (4.3). The diagonal entries of (4.3) represents the sum
of available line capacities at each bus.




 1 1


 1 0



LCM =  0 1



 0 0


0 0


1 0 0 1 




0 0 0 0 




0 1 0 0 



1 0 1 0 



0 1 1 1 

C1
0
0
0
0
0





0 



C2 0
0
0
0 




0 C3 0
0
0 




0
0 C4 0
0 



0
0
0 C5 0  


0
0
0
0 C6

1 1 0 0 0 


1 0 1 0 0 



1 0 0 1 0 



0 0 1 0 1 



0 0 0 1 1 


1 0 0 0 1

0

0

0

0



(4.2)



0
0
0
0
 C1 + C2 + C3



0
C1 + C6
0
0
0



LCM = 
0
0
C2 + C4
0
0




0
0
0
C3 + C5
0


0
0
0
0
C4 + C5 + C6
44
















(4.3)

1

4
C3

G1

C2
3
C4
C5

C1
5

2
G2

C6

Figure 4.2: A hypothetical system to demonstrate the application of the state space reduction
technique

This heuristic technique was first developed by the authors in [10]. The heuristic technique is different from the pruning techniques, which are already presented in the literature,
in the sense that it classifies the state space into success, failure, and unclassified subspaces.
The method developed in [10] has given satisfactory results in evaluating energy indices due
to the fact that it sums up the probabilities of the visited states, and not the estimated
values of the indices. In section 4.5, we introduce a new method, which is developed based
on the complementary concept, to estimate the energy indices with less computational effort.

4.3

Power System Reliability Evaluation using Binary
Particle Swarm Optimization Search Method

Population-based intelligent search methods have been amply used in power system studies [45–47]. Genetic algorithm and bacterial foraging have been applied in power system
reliability evaluation in [48, 49]. Applications of binary particle swarm optimization in

45

power system studies have reviewed in [50]. Binary particle swarm optimization has been
used as a searching tool for success or failure states in power system reliability evaluation
in [14, 38–43, 51–55]. Binary particle swarm optimization has been used in the reliability
evaluation of distribution systems and microgrids in [56–59]. The objective function of the
binary particle swarm optimization can be single objective or multi-objective. The two objective functions were intended to prevent the particles from being trapped to one corner
of the state space. The weighting factors associated with the objective functions have been
carefully selected in order to keep the particles to search in the entire state space [12, 13].
Even with two objective functions, if the weighting factors are not carefully chosen, the
particles might be trapped to one corner of the search space. To circumvent this difficulty,
an effective approach that adjusts these weighting factors and using an appropriate fitness
function are presented in this work.

4.4

Non-Sequential Monte Carlo Simulation vs. Binary Particle Swarm Optimization Search Method

One of the differences between using binary particle swarm optimization in power system reliability evaluation and non-sequential Monte Carlo simulation is that binary particle swarm
optimization sums up the values of the reliability indices of the visited states, whereas in
non-sequential Monte Carlo simulation, the reliability indices are evaluated by dividing the
number of encountered failure states to the number of samples. It is well known that nonsequential Monte Carlo simulation takes long time to calculate reliability indices for the
designated accuracy; however, it is obvious that by using non-sequential Monte Carlo simulation, no need to visit all system states since it calculates the estimated values not the
46

exact values. Therefore, through the sampling process, Monte Carlo simulation reaches the
designated index from both sides, under estimation and over estimation. On the other hand,
intelligent search methods such as binary particle swarm optimization accumulate the values
of the designated indices as the search progresses. Therefore, all failure states or success
states need to be visited to calculate reliability indices which are computationally expensive
if it is not impossible for large systems. Hence, the values of these indices will be always
less than the exact values unless the search method visits all the failure states or all the
success states or the entire state space. Therefore, another way to overcome this necessity
is required.
New indices have been introduced in [12, 13] which are normalized LOLP, (LOLPnorm )
and projected LOLP, (LOLPproj ) to eliminate the necessity of visiting all the failure and
success states. The LOLPnorm is the ratio between the sum of the probabilities of the
failure states (LOLP) to the sum of probabilities of the states encountered so far (1+LOLPDLOLP); where DLOLP is the dual value of the LOLP which is one minus the sum of
probabilities of the success states encountered so far (1-sum of the probabilities of the visited
success states). The LOLPproj is the crossings between the tangents of the LOLP and
DLOLP. The normalized LOLP can be expressed as follows,

LOLPnorm =

LOLP
1 + LOLP − DLOLP

(4.4)

Again, LOLPnorm and LOLPproj are the same as the weighted LOLP index used in
Monte Carlo simulation. However, LOLPproj converges to the approximate value of the
LOLP faster than LOLPnorm [13]. It can be noted that even though binary particle swarm
optimization and Monte Carlo simulation both evaluate the weighted index, the number of

47

states to be visited using binary particle swarm optimization is significantly less than the
number of states to be visited using Monte Carlo simulation [13].
The binary particle swarm optimization is used in this work to search intelligently for
the most relevant states rather than randomly searching for failure states. Also, we present
comparisons between using binary particle swarm optimization and the conventional Monte
Carlo simulation.
The assumptions made in [12,13] are only valid for calculating loss of load probability not
for the energy, frequency and duration indices. In this work, we have used the complementary
concept (section 4.5) to evaluate all the well-known composite systems reliability indices.

4.5

The complementary Concept

The concept of complementary has been applied in several disciplines to calculate/estimate
a parameter from its complementary value. The complementary concept in composite power
system reliability evaluation can be defined as “every reliability index that can be evaluated
from the failure states, has a complementary value that can be evaluated from the success
states as long as the boundaries of an index and its complementary value are known such as
probability and energy indices or the index can be equally likely determined from the success
or failure states such failure frequency indices”. For example, the boundaries (minimum and
maximum values) of the probability of system failure and the probability of system success
are between 0 and 1. Also, the boundaries of expected load curtailment and the expected
load supplied are between 0 and the peak load. Frequency indices can be evaluated either
from the failure states that transit upwards to success states or from success states that
transit downwards to failure states. The complementary concept estimates reliability indices

48

from visited states; rather than the entire state space.

4.5.1

Theoretical Background

The complementary concept starts with the fact that the sum of the probabilities of the state
space is one. Therefore, the sum of the probabilities of the success states and the failure
states is also one. The number of states in a state space is 2n where n is the number of
the components. Most power system components (generation and transmission) are highly
reliable. Therefore, the state that all system components are in the up state has the highest
probability. Not only that, the probability of a state with only one component being in the
down state usually less than the probability of the state with all the components are in the
up states of order of 1/qi where qi is the unavailability of component i. This probability
order can be generalized for the rest of the states, that is the probability of a state that
has two components in the down state (i and j) is less than the probability of a state
with all the components are in the up state of order of 1/(qi qj ). Therefore, as the number
of the components in the down state increases, the probability of a state decreases with
reciprocal of the multiplication of the un-availabilities. Therefore, if an intelligent search
method discards the states with very small probabilities and captures the states that have
most effect on system reliability, the number of states that need to be evaluated can be
decreased significantly.
For example, consider a system with ten identical components with availability of 0.95.
The probability of all the components are in the up state is 0.59874 and the probability of
one component being in the down state is 0.03151. Therefore, the probability order between
the two states is 19. The probability of a state with two components are in the down state
is 0.00166. The probability order between this state and the state with all components are
49

up state is 361. The probability of the state with all components are in the down state is
9.76563×10−14 . Therefore, even with ignoring the states that have small probabilities, we
can accurately evaluate the behavior of the system.
For this system, if we consider at most one component being in the down state, we need
to evaluate 11 states out of 1024 states and the sum of the probabilities of these states is
0.91386. If we consider at most two components are in the down state, we need to evaluate
56 states out of 1024 states and the sum of the probabilities of these states is 0.98850. This
behavior is trivial for small systems with identical components; however, for large systems
with non-identical components, it is difficult to trace the behavior of the system analytically.
Therefore, if we allow some error tolerance, an intelligent search method can capture the
states that have most effect on system reliability. Having said that, the reliability of the
system can be approximately evaluated based on these states. The indices evaluated from
these states, however, will not be equal to the exact indices. Therefore, a method to estimate
these indices from the specified states is of necessity. In this work we introduced the concept
of complementary to estimate reliability indices.
As the algorithm searches through the state space, the reliability indices and their complementary values become close to each other. For example, the difference between the
probability of system failure calculated from the failure states, q, and the probability of system failure calculated from the success states, (1 − p), goes gradually to zeros as the search
progresses. Since, the search method is not intended to visit the entire state space and the
reliability indices can approximately evaluated from the relevant states, the reliability indices
can be referred to the relevant subspace. Since the sum of probabilities of the entire state
space is 1 and the sum of probabilities of the visited states is q + p, then a reliability index

50

can estimated over the evaluated subspace as follows,

fˆi (x) = fi (x) ×

1
q+p

(4.5)

where, fi (x) is the reliability index under consideration and fˆi (x) is its estimated value.

4.5.2

Mathematical Justification

If we evaluate the entire state space and if we denote the loss of load probability as q and
the probability of load availability as p, then the following relationship holds,

p+q =1

On the other hand, if we evaluate part of the state space, then q + p < 1. In this case the
estimated value of q, (ˆ
q ), and estimated value of p, (ˆ
p), will be less than or equal to the exact
values, that is qˆ ≤ q and pˆ ≤ p.
If the search tool were able to visit all the failure states without the need to visit the entire
state space, then qˆ = q and pˆ < p. On the other hand, if the search tool were able to visit
all the success states without the need to visit the entire state space, then qˆ < q and pˆ = p.
Searching for all the failure states or all the success states is computationally expensive
and equivalent to state space enumeration. If the search tool intelligently searches uniformly
for both success and failure states and if the algorithm does not evaluate the entire state
space, then qˆ < q and pˆ < p.
If the algorithm searches for the states that have the most effect on system reliability,
then a high cumulative probability of the state space can be evaluated through a small

51

fraction of the state space. Therefore, q and p can be approximated as,

q+p≈1

Through the search process, the q index represents the sum of the probabilities of the
failed states visited so far and p index represents the sum of the probabilities of the success
states visited so far. The sum of the probabilities of the states visited so far is q + p.
Therefore, the estimated value of the probabilities of failure over the visited subspace can
be expressed as follows,
qˆ =

q
q+p

and the estimated value of the probabilities of the success over the visited subspace can be
expressed as follows,
pˆ =

p
q+p

Same relations can be applied for the other indices that can be calculated from the failure
states as well as from the success states as shown in chapter 3.

4.6

Particle Swarm Optimization

Particle swarm optimization (PSO) is a population-based intelligent search method, which
has been proposed by Kennedy and Eberhart in [60]. Later, Kennedy and Eberhart proposed
a discrete binary particle swarm optimization to solve combinatorial optimization problems
[61]. PSO has been proven to be an effective optimization technique and has, therefore, been
used in the presented work. We use the binary PSO to search for success and failure states
rather than looking for optimal solution. Further, generation and transmission lines states
52

are represented as up states or down states so that up states are represented by 1’s and down
states are represented by 0’s. Moreover, we apply bi-objective function that searches through
the entire state space for success and failure states. Such objective function is also capable
to prevent particles from being trapped to one area of the search space. It is worth noting
here that the size of the state space is 2n , where n is the number of system components.
Therefore, even for medium systems, the number of states is very large and it is impractical
to test the entire state space. Hence, the concept of the complementary is used in this
work so that the condition of searching for all failure states over the entire state space is
eliminated.
The binary PSO technique searches through the state space and tests the visited states
for a possibility of load curtailment. In performing the search process, for every iteration,
particle velocity vi or the direction of movement of particle i from position xi can be governed
by an objective function (sections 4.6.1 and 4.6.2).
The change in particles positions can be defined by a sigmoid limiting transformation
function and a uniformly distributed random number in [0,1] as following,

xid =





1, rand() < S(vid )

(4.6)




0, otherwise
where, xid is the dth component of particle i and S(vid ) is the sigmoid function of dth ’s
component of particle i which can be expressed as,

S(vid ) =

1
1 + e(−vid )

53

(4.7)

X2
- Large probabilities
- Small load curtailments

Two conflicting forces

- Small probabilities
- Large load curtailments
X1

Figure 4.3: State space representation with two conflicting objective functions
The objective function of the BPSO can be single objective (section 4.6.2) or multiobjective (section 4.6.1). In [38–43], single objective function has been used which is minimum load curtailment to search for success states to be pruned from the state space. BPSO
with two objective functions has been proposed in [12, 13] to calculate reliability indices.
In [12, 13], one objective function has been used to maximize load curtailments and the
other has been used to maximize states probabilities. These two conflicting objective functions have been intended to prevent the particles from being trapped to one corner of the
state space. The representation of the state space with these two objective functions is
depicted in Figure 4.3 where x1 and x2 represent the possible states for two system components. State space pruning using BPSO with two objective functions has been introduced
in [44]. The BPSO has been used with two objective functions to prune the success states
from the state space and then to apply Monte Carlo Simulation on the remaining part of the
state space [44]. Even with two objective functions, if the weighting factors are not chosen
carefully, the particles might be trapped to one corner of the search space. Therefore, using
multi-objective function does not guarantee that the particles will fly on the entire state
space unless the right weighting factors are chosen as shown in section 4.7.2.
54

4.6.1

Multiple Objective Functions

The two objective function can be described as follows [12, 13, 61]:

vi = vi + rand() × c1 × (P pbesti − xi )
+ rand() × c2 × (P gbest − xi )

(4.8)

+ rand() × c3 × (Cpbesti − xi )
+ rand() × c4 × (Cgbest − xi )
where rand() is a uniformly distributed random number between [0,1], P pbesti is the particle
best position from the probability of a state prospective particle i has ever encountered,
P gbest is the group of particles best position from the probability of a state prospective the
group has ever encountered, Cpbesti is the particle best position from the load curtailment
of a state prospective particle i has ever encountered, Cgbest is the group of particles best
position from load curtailment of a state prospective the group has ever encountered and c1 ,
c2 , c3 and c4 are acceleration factors.

4.6.2

Single Objective Function

To circumvent the difficulty of choosing appropriate weighting factors, a fitness function
based on the severity of the states should be used. In this work, we propose a fitness function
that is based on the risk index for each state which is the product of state probability and
amount of load curtailment. This fitness function eliminates the need to use multi-objective
function and guides the swarm to search around the states that have most effect on system
reliability and prevents the swarm from being trapped to one corner of the search space.

55

The proposed fitness function can be expressed as follows,

ρ(x) = p(x) × Lc(x)

(4.9)

where p(x) is the probability of state x and Lc(x) is the amount of load curtailment of state
x. The objective function becomes,

vi = vi + rand() × c1 × ρpbest−i − xi

(4.10)

+ rand() × c2 × ρgbest − xi
where ρpbest−i is the particle best position from the risk of a state prospective particle i has
ever encountered and ρgbest is the group of particles best position from the risk of a state
prospective the group has ever encountered.

4.7

Dynamically-directed Binary Particle Swarm Optimization Search Method

In performing the heuristic technique, some failure events may not be captured. In this
context, three failure scenarios can exist. In the first scenario, an area of the system is isolated
from the generation areas (area with loads only) and the transmission lines connecting the
loads of this area are sufficient to carry the loads. In the second scenario, an area of only load
buses is connected to the rest of the system through a tie-line that cannot withstand the total
load of this area, but the transmission lines can do. In the third scenario, line impedances
force the power to flow through some capacity limited lines and the power flow will not
converge without load curtailment. These scenarios cannot be detected by considering the
56

power carrying capabilities of the transmission lines. Therefore, a search tool should be
used to capture these events. However, the search space can be bounded by definite failure
subspace and definite success subspace, which makes the search method converges to these
events with lower computational effort.

4.7.1

State Space Description

Figure 4.4 shows the state space of a three dimension system. The above reduction technique
can detect the definite failure subspace, but cannot determine the boundary of the definite
success subspace. However, since the reliability indices and their sensitivities are evaluated
based on the failure states, analysis of the entire success subspace is not important. If the
boundary of the success subspace is known, then an algorithm can be developed to search
for failure states outside the success subspace. A large part of the success subspace can
be truncated by considering all the transmission lines in the up state and determining the
vector of minimum generation that satisfies the load. Any state in the state space that has
a generation vector larger than or equal to the minimum generation vector is guaranteed to
be a success state. Therefore, the boundary of the success subspace can be used to control
the binary particle swarm optimization search method. Also, most of the unclassified failure
states lay around the boundary of the definite failure subspace. Therefore, an intelligent
search method that is bounded between the boundaries of the success and failure subspaces
can efficiently converge to the failure states in the unclassified subspace. In this work, we
use a method based on a dynamically directed binary particle swarm optimization to search
for these states.

57

X3
Failure Events with High
Probabilities and Low Load
Curtailments
Definite Failure
Events Subspace
Boundary of Definite
Success Events
Subspace

X1

X2

Figure 4.4: State space representation of three dimension system (three components with
multiple outcomes)

4.7.2

Directed BPSO

By examining the state space, all possible outcomes of the combination of generation and
transmission, the probabilities of the upper part states are larger than the probabilities of
the lower part states. Further, load curtailments of the lower part states are larger than
load curtailments of the upper part states. If the swarm was to fly on the upper part, it is
unlikely to encounter failure states and the opposite is true for the lower part. As mentioned
above, using multi-objective function does not guarantee that the particles will fly on the
desired search space unless the weighting factors are chosen carefully. Figure 4.5 and Figure
4.6 show two situations for which the weighting factors of the objective function were badly
chosen. Figure 4.5-(a) and Figure 4.6-(a), show the behavior of the particles and Figure
4.5-(b), and Figure 4.6-(b) show the expected profile of the reliability indices. In Figure 4.5
the swarm is trapped to the upper corner of the state space. In this case, the normalized

58

LOLP, (LOLPnorm ) as well as the projected LOLP, (LOLPproj ), will be under estimated.
As the particles continue to fly on the upper corner, most of the visited states will be success
states and more likely no failure states will be captured. Therefore, no much change in the
LOLP and the 1-DLOLP is sharply increasing. On the other hand, if the swarm is trapped
into the lower corner of the state space as shown in Figure 4.6 most of the visited states will
be failure states and unlikely to encounter success states. In this case, LOLP, (LOLPnorm )
as well as the projected LOLP, (LOLPproj ), will be overestimated. The LOLP index will
continue to increase whereas the 1-DLOLP index remains almost constant.
x2

(k  1 )th  n th iterations

k th iteration
(k  1 )th iteration

(k  2 )th iteration
(a)

(k  3 )th iteration

x1

LOLP
1  LOLP - DLOLP
Reliability indices

0.5
1  DLOLP

(b)

Exact LOLP
LOLP

Number of iterations

Figure 4.5: Reliability indices profile in case the particles are trapped to the upper corner of
the state space

59

x2
(k  3 )th iteration

(k  2 )th iteration

(k  1 )th iteration

(a)

k th iteration

(k 1 )th  to  n th iterations

x1

LOLP
1  LOLP  DLOLP
Reliability indices

0.5
(b)

Exact LOLP
NLOLP

LOLP
Number of iterations

Figure 4.6: Reliability indices profile in case the particles are trapped to the lower corner of
the state space
To direct the BPSO to fly on the unclassified subspace and overcome the difficulty of
choosing the weighting factors as well as the change in probability limit, we propose a
technique that tracks both the normalized LOLP and the number of the new visited states
or the mutation rate. This technique was developed based on the normalized loss of load
probability index (LOLPnorm ). The LOLPnorm has been introduced in [12, 13] to estimate
the LOLP index without the need of searching the entire state space. After performing few
iterations, if LOLPnorm continues to increase that means the particles are trapped in the
lower corner of the search space. On the other hand, if LOLPnorm continues to decrease,

60

that means the particles are trapped in the upper corner of the search space. Also, if no
new states are visited and the designated reliability index does not converge to the specified
accuracy, then Vmax is limiting the probability of bit change. Also, if the algorithm discovers
that the swarm entered the classified success or failure subspaces, the weighting factor are
adjusted accordingly as shown in Figure 4.7. The procedures of preventing the particles from
being trapped into one area are described in detail in the solution algorithm.
X2
Definite Success
Subspace

Definite Failure
Subspace

X1

Figure 4.7: Redirecting the swarm in case it enters the success or failure subspaces

61

Chapter 5
Sensitivity Analysis in Composite
System Reliability Using Weighted
Shadow Prices

5.1

Introduction

In this chapter, a shadow-price based method is proposed for performing sensitivity analysis
in composite system reliability. Shadow prices (or Lagrange multipliers) are used as a decision
making tool in identifying the system’s component or area that has the highest impact
on system reliability; and which area or component need to be reinforced to enhance the
overall system reliability. In determining components’ shadow prices, power system grids
are represented by the DC power flow model (chapter 2) and solving for the dual solution
with an objective function of minimizing load curtailments. The proposed method relies
on component availability data, generation capacity, transmission line capability and load
scenarios. In most of the practical applications, the sensitivity evaluation techniques include
approximations in generating capacity model, load model and the evaluation techniques.
The proposed method uses a Monte Carlo simulation method that does not necessitate such
approximations.

62

Monte Carlo simulation has been amply used in the literature in evaluating power system
reliability. The use of Monte Carlo simulations also serves as a benchmark for testing the
validity of the proposed state space reduction technique and the dynamically directed particle
swarm optimization search method in evaluating power system reliability. In this chapter,
we use Monte Carlo simulation without using any device to reduce the computation time.
Therefore, the results of using the proposed techniques will be compared with results reported
in this chapter.

5.2

Shadow Price Theory and Description

Shadow price principle was first proposed in 1930s by Conte Petrovic as an application of
linear programming optimization technique [62]. Shadow price is defined as the sensitivity
of the value of the objective function to the change in the right hand side of the linear
programming problem constraints [27]. The objective is to minimize load curtailments and
the constraints are power system network components’ limitations such as generation capacities and transmission lines’ carrying capabilities. From generation condition or status,
total available power generation can be determined. Also, from transmission line availability
and carrying capacity, the constraints can be set. From this solution, shadow prices of the
generating units and transmission lines can be calculated.
In general, suppose we have m constraints with n variables, the standard maximum linear
programming problem can be formulated as follows:
n

max(Z) =

cj x j
j=1

63

(5.1)

Subject to the following constraints,
n

aij xj ≤ bi

(5.2)

xj ≥ 0

(5.3)

j=1

The matrix of the coefficients of the constraints is,




 a11 a12


 a21 a22

A=
 ..
..
 .
.


am1 am2

. . . a1n 


. . . a2n 


.. 
..
.
. 


. . . amn

(5.4)

The right hand side vector b of the constraints consists of m constants,
T

b=

(5.5)

b1 , b 2 , · · · , b m

The row vector of the objective function c consists of n coefficients
T

c=

c1 , c 2 , · · · , c n

(5.6)

The dual of this standard maximum problem is the standard minimum problem, that is,
m

min(W ) =

bi yi
i=1

64

(5.7)

Subject to the following constraints,
n

aij yj ≥ cj

(5.8)

yi ≥ 0

(5.9)

j=1

If the slack variables of equations (5.2) and (5.3) are imported, the standard linear programming will have the following form,

max(Z) = CX

(5.10)

AX = b

(5.11)

X≥0

(5.12)

With the equality constraints,

The optimal and feasible solution of (5.7) can be expressed as,

Xb = B −1 b

(5.13)

where B is the optimal basis at the optimal and feasible solution and Xb is the basic variable
sub-vector. In the standard problem of (5.2), the significance of the dual variables yi is
described as the ith shadow price. In other words, the yi of the dual problem are interpreted
here as the shadow prices of the linear constraints bi .

65

5.3

Solution Algorithm

This section describes the approach used for calculating shadow prices of various system
components. Shadow prices can be calculated for one load scenario as well as for multiple
loading scenarios. However, to consider the shadow price as a reliability index, system
states are to be sampled till the shadow prices of all system components converge to the
pre-specified tolerance.
Since the calculation of shadow prices requires only system state and not the state duration, system state sampling approach is used in determining system states. Components
outages are assumed independent events and their behavior is categorized by a uniform
distribution [0, 1].
For two state components, if the state of the ith component is denoted by Si and component unavailability is denoted by qi , system state vector S will be,

S = [S1 , S2 , · · · , Si , · · · , Sn−1 , Sn ]

where, n is the number of components in the system.
A uniformly distributed random variable Ui under [0, 1], using improved prime number
multiplicative congruential generator method, is drawn for every component such that,

Si =




1

(upstate)

if Ui ≥ qi



0

(downstate)

if 0 ≤ Ui < qi

For three state components such as generating units with derated states, let P di denote

66

the probability of the derated state of the ith generating unit such that,

Si =





1




2






0

(up state)

if Ui ≥ (qi + P di )

(derated state)

if qi ≤ Ui < (qi + P di )

(down state)

if 0 ≤ Ui < qi

System load is modeled by multistep model (cluster technique). Probabilities of load
steps P Li are put successively between [0, 1]. By drawing a uniformly distributed random
number, the load step can be determined according to the position of the random number
of the cluster discrete probability distribution function.
After sampling system states, system bus and line data are modified according to the
new state. Using the DC power flow model and linear programming, components shadow
prices can be calculated for each state. The procedures are illustrated in Figure 5.1.

5.4

Case Studies

The proposed method is applied on the IEEE RTS [63], the Modified IEEE RTS, and the
Saskatchewan Power Corporation in Canada (SPC) and a hypothetical six bus system [64].
System descriptions are given in the appendix. Since IEEE RTS has two 400 MW and one
350 MW units which have high capacity compared to the rest of the system, it is desirable to
consider the unit derated states which is given in the modified test system [63]. Accordingly,
the proposed method is applied on the base case and the modified IEEE RTS test systems.
Also, compared to the generating units, IEEE RTS transmission system is too reliable or over
designed [63]. Due to these reasons, two additional case studies were conducted to enforce
the transmission lines to become binding constraints instead of being always redundant
67

Start

Read System Data

Sample System States

Is there any
Change?

No

Yes
Modify Bus Data and Line Data

Perform Linear Programming
and Calculate the Shadow Prices

No

Is there a Feasible
Solution?
Yes

Calculate the Average Shadow Prices

ε < Tolerance?

No

Yes
Stop

Figure 5.1: Flow chart of the proposed method

68

constraints. Furthermore, for the purpose of comparison, sensitivity analysis of the IEEE
RTS is compared with those given in [32]. Another case study was conducted to determine
the reliability indices using Monte Carlo simulation. The results of using Monte Carlo
simulation will be compared with the results of using the state space reduction technique
and the dynamically directed particle swarm optimization search method.

5.4.1

Case Study I: IEEE RTS Base Case

For the base case, generating units were represented by two states (up state and down state).
The load levels were modeled using 20 clusters load levels. To model bus loads, the RTS peak
load of 2850 MW was taken as 1 p.u. and bus loads sampled to 20 clusters as a percentage
of 1 p.u. value.
Shadow prices values are given in column 2 of Table 5.1 for the generating units, column
3 of Table 5.2 for forward direction of transmission lines and column 3 of Table 5.3 for reverse
direction of the transmission lines. Also, it should be noted that, transformers were modeled
as transmission lines with length zero and they are placed between buses 3-24, 9-11, 9-12,
10-11 and 10-12.

5.4.2

Case Study I: Modifed IEEE RTS

As mentioned previously, the two 400 MW and one 350 MW generating units have considerably large capacity compared to the other components of the system. These three units
were assumed to have derated states as it is given in [63]. In other words, every unit has
three states (full capacity state, derated state and down state). In simulating these states,
a cumulative distribution function was formed for every generating unit. For example, 400

69

MW has a cumulative probability with three steps; down, derated and full capacity levels.
Shadow prices values are given in column 3 of Table 5.1, column 4 of Table 5.2 and column
4 of Table 5.3.
Table 5.1: Generators Shadow Prices

5.4.3

Bus

Case Study

Case Study

Case Study

Case Study

No.
1
2
7
13
15
16
18
21
22
23

I
0.40310
0.40378
0.15612
0.36767
0.36832
1.54987
3.69131
4.12508
0.36371
0.36551

II
0.43786
0.43881
0.17088
0.40017
0.39590
2.07930
3.18961
3.60711
0.42293
0.38899

III
0.18527
0.18739
0.00018
0.16363
0.13240
1.16986
3.76693
3.62967
0.11530
0.15000

IV
0.48389
0.48390
0.03040
0.45123
0.44642
1.50237
3.33811
16.9437
0.44142
0.44435

Case Study III: Increasing Generation Capacity and Loading
Level

Since the transmission lines are too reliable and have high capacity compared to the generating units, they are forced to become binding constraints by increasing generating units’
capacities by two folds and increasing buses’ loads by a factor of 1.8. These modifications
have forced the lines to carry power near to their rated capacities. Shadow prices values of
the system components are given in column 4 of Table 5.1, column 5 of Table 5.2 and column
5 of Table 5.3.

70

Table 5.2: Transmission Lines Shadow Prices (Forward)
From

To

Case Study

Case Study

Case Study

Case Study

Bus
1
1
1
2
2
3
3
4
5
6
7
8
8
9
9
10
10
11
11
12
12
13
14
15
15
15
15
16
16
17
17
18
18
19
19
20
20
21

Bus
2
3
5
4
6
9
24
9
10
10
8
9
10
11
12
11
12
13
14
13
23
23
16
16
21
21
24
17
19
18
22
21
21
20
20
23
23
22

I
0.01910
0.00513
0.03250
0.01127
0.10014
0.02112
0.13659
0.06473
0.31012
0.29032
0.00945
0.08640
0.26387
0.28179
0.13933
0.19194
0.20482
0.05855
0.00888
0.04368
0.06323
0.07747
0.03233
0.05662
0.00897
0.03370
0.02765
0.01830
0.00999
0.08521
0.03194
0.00389
0.01249
0.01843
0.02044
0.02704
0.01535
0.06845

II
0.02797
0.03031
0.01738
0.02497
0.06714
0.00654
0.33778
0.11148
0.18250
0.14138
0.00413
0.04409
0.08410
0.34020
0.13406
0.21197
0.16501
0.01416
0.00847
0.06997
0.08933
0.03390
0.04106
0.03764
0.04632
0.05385
0.01104
0.02639
0.02918
0.02969
0.24006
0.09870
0.03610
0.02143
0.05901
0.07893
0.02283
0.01126

III
0.02478
0.05207
0.01347
0.04716
0.03747
0.00682
0.33118
0.03406
0.01780
0.20795
0.00594
0.04403
0.05354
0.06319
0.21986
0.04182
0.20676
0.09561
0.03080
0.04625
0.10376
0.03663
0.24144
0.04132
0.02449
0.13175
0.14698
0.14053
0.02474
0.03175
0.16942
0.02135
0.01660
0.03792
0.03024
0.00847
0.02537
0.02745

IV
0.02521
0.03879
0.01494
0.02783
0.07840
0.01670
0.29264
0.03459
0.00904
0.08561
0.02115
0.06676
0.02620
0.27834
0.22375
0.26880
0.33482
0.02822
0.02819
0.02875
0.05803
0.01391
0.07660
0.01650
0.00793
0.01835
0.06034
0.02727
0.01844
0.02398
0.08765
0.02919
0.01697
0.05032
0.16131
0.02087
0.01261
0.07253

71

Table 5.3: Transmission Lines Shadow Prices (Reverse)
From

To

Case Study

Case Study

Case Study

Case Study

Bus
2
3
5
4
6
9
24
9
10
10
8
9
10
11
12
11
12
13
14
13
23
23
16
16
21
21
24
17
16
17
17
18
18
19
19
20
20
21

Bus
1
1
1
2
2
3
3
4
5
6
7
8
8
9
9
10
10
11
11
12
12
13
14
15
15
15
15
16
19
18
22
21
21
20
20
23
23
22

I
0.02465
0.00641
0.03412
0.01545
0.10526
0.01914
0.10015
0.05499
0.07976
0.05761
0.55643
0.08052
0.02616
0.21403
0.07665
0.23862
0.52747
0.04484
0.00890
0.02962
0.03859
0.07103
0.02340
0.06181
0.00620
0.02223
0.04418
0.01508
0.01000
0.08367
0.01489
0.00389
0.00854
0.01594
0.02044
0.01843
0.01104
0.03180

II
0.03037
0.03088
0.01854
0.02762
0.06962
0.00831
0.28250
0.10741
0.11993
0.06277
0.55998
0.03544
0.01965
0.23538
0.07778
0.18203
0.22230
0.00789
0.00833
0.04041
0.05952
0.02514
0.03051
0.03872
0.02904
0.03742
0.01751
0.02060
0.02766
0.02786
0.06379
0.09616
0.02988
0.01601
0.05406
0.06626
0.01526
0.01126

III
0.02492
0.05251
0.03744
0.06835
0.03698
0.00774
0.24665
0.02991
0.01688
0.05999
0.66622
0.02852
0.03037
0.04538
0.13748
0.03000
0.13929
0.03609
0.03069
0.02747
0.05860
0.03419
0.02585
0.04110
0.01828
0.10214
0.16761
0.04628
0.04628
0.02223
0.03157
0.12308
0.02135
0.01660
0.02937
0.02618
0.00462
0.01694

IV
0.02345
0.03763
0.02722
0.03738
0.09066
0.01517
0.24444
0.02632
0.00853
0.06797
0.93444
0.05715
0.02211
0.22259
0.16371
0.20555
0.28966
0.02184
0.02696
0.01715
0.03715
0.01051
0.05528
0.01688
0.00520
0.01022
0.07624
0.01816
0.01816
0.01563
0.02138
0.02764
0.02533
0.01185
0.04566
0.15679
0.01408
0.00780

72

5.4.4

Case Study IV: Reducing Lines Capacity

In this case, the carrying capacities of the transmission lines are decreased to 75% of their
rated capacity. This modification was assumed to examine the effect of the lines capacities
on the system reliability. The shadow prices values of the system components are given in
column 5 of Table 5.1, column 6 of Table 5.2 and column 6 of Table 5.3.

5.4.5

Case Study V: Calculation of the Reliability Indices of Four
Test Systems

In this case study, Monte Carlo simulation method was used to estimate the well-known annualized reliability indices of the IEEE RTS [63], the Modified IEEE RTS, and the Saskatchewan
Power Corporation in Canada (SPC) and a hypothetical six bus system [64]. The annualized
reliability indices are evaluated at the system peak load. The annualized indices are shown
in Table 5.4.
Table 5.4: Annualized Indices of the Tested Systems.
qˆ

dˆn

Fˆ

ˆq
L

eˆn

Tˆc

Number of

MW

occ./yr

h/yr

MWh/yr

hour

sampled states

The Hypothetical Six Bus System
0.002400 0.034630

1.34832

21.02400

303.3588

15.59274

21,981

41.07354

73,925

The IEEE RTS System
0.08539

14.19236 18.16174

745.967

123,984

The Modified IEEE RTS System
0.069052 11.88175 17.14523

603.238

103,799

35.18403

67,884

Saskatchewan Power Corporation (SPC) System
0.001081 0.039098

3.23125

9.44362

341.560

2.92259

122,941

The results of the annualized indices using Monte Carlo simulation will be used in the
73

next chapters as a benchmark to validate the proposed methods. Also, the number of samples
will be used to test the efficiency of the proposed state space reduction technique.

5.5

Results and Discussion

From Table 5.1 it is clear that the generating unit (400 MW) at bus 21 has the highest
shadow price followed by the shadow prices of the generating units of buses 18 and 16 for the
four study cases. Consequently, these generating units have the largest impact on system
reliability. In other words, these buses are at the highest risk and the most beneficial way
to make the system more reliable is to increase the availability of these units by adding
new units, decreasing repair time and/or using sophisticated methods to decrease units’
failure rate. When the derated states were introduced, these shadow prices decreased due
the increase in the availability of the two 400 MW and one 350 MW generating units.
In comparison to the case study I, the shadow prices of the generating units of buses 21,
18 and 16 of the case study III have slightly changed due to the unequal increase in the total
installed capacity and the peak load. Such change in the installed capacity and the peak
load has forced the transmission lines to have larger shadow prices than it was in case study
I instead of forcing the generating units.
For the case study IV, the shadow price of the generating unit at bus 21 has increased
for more than four times than it was in case study I whereas the generating units at buses 18
and 16 stayed almost the same. In addition, while buses 18 and 21 have identical generating
units, they have different shadow prices; which indicates that not only the availability of
components has an impact on systems’ reliability, but also, components’ location in the
grid. By looking back to the network configuration, bus 21 is connected to bus 15 through

74

two lines with normal capacities of 500 MW and to bus 18 also with two lines with normal
capacities of 500 MW and they have the highest loadings in the network.
From Table 5.2 and Table 5.3, the overhead transmission line between buses 7 and 8 has
the highest shadow price followed by the transformer between buses 10 and 12 and the line
that connects buses 5 and 10. From the network configuration, line 7 - 8 is the only line
that connects bus 7 to the rest of the system. Also, the nominal real power injection at bus
7 is 175 MW (300 MW - 125 MW) and the normal capacity of the line 7 - 8 is 175 MW.
Therefore, line 7 - 8 is almost always carrying its rated power. In addition, the transformer
between buses 10 and 12 connects bus 23 which have the highest generation to bus 10 which
in turn connected to load buses (buses 5, 6, 8 and 14 through bus 11).

5.6

Summary

This chapter has introduced the use of shadow price concept in estimating sensitivity of
the reliability indices with respect to component parameters. The shadow price can be
obtained from the dual solution of the standard linear programming problem. Four study
cases have been conducted on the IEEE RTS test system. It has been found that shadow price
analysis is an efficient tool for composite generation and transmission planning, expansion
and operation studies when compared to the sensitivity analysis technique. Moreover, it
gives essential information about which component is in unsafe state from reliability point
of view and the most effective way to enhance system reliability.

75

Chapter 6
Power System Reliability Evaluation
and Sensitivity Analysis Using the
State Space Reduction Technique and
Population-based Intelligent Search
Methods
The task of power system reliability evaluation is becoming complicated due to the increase
in the complexity of the power industry infrastructure. The time and computation efforts
to evaluate power system reliability indices and their sensitivities with respect to component parameters and system operating limits are of concern in both planning and operation
prospectives. In this chapter, we use the heuristic, which is introduced in 4 to prune the
state space and to reduce the computational efforts in evaluating the well-known reliability
indices and their sensitivities with respect to component parameters and system operating
limits. In this chapter we use the DC power flow model which is given in chapter 2.
This chapter serves as a validation for the proposed state space reduction technique and
the use of population-based intelligent search methods in power system reliability evaluation.

76

The results of these methods will be compared with results of using Monte Carlo simulation
which is presented in chapter 5. The method presented in this chapter is based on particle
swarm optimization search method and it is applicable to the other population-based intelligent search methods. In addition to the validation of the results, this chapter also shows
the reduction in the computational time and effort. The proposed method is applied on the
IEEE RTS, the Modified IEEE RTS and PSC.

6.1

Discussion on the Reduction in the Computational
Time

In calculating reliability indices, an optimization problem usually needs to be performed for
every system state with different loading scenarios. The process of preparing and modeling
power system data to be tested for reliability analysis is performed one time. However,
running an optimization problem is required for every system state or for any change in
system parameters. From a study of wide variety of system sizes and configurations, the
optimization problem usually takes most of the calculation time. For instance, on average,
the time required to solve an optimization problem for the systems we have tested took
more than 95% of the calculation time for each system state. By performing the reduction
technique, not all the states need to be passed to the optimization problem. In other words,
if a state is guaranteed to be a success state, no useful results can be shown by solving an
optimization problem. On the other hand, failure states should be passed to the solver to
calculate bus indices. It should be pointed out here that if only system indices are to be
evaluated and not bus indices and the algorithm was able to capture most of the failure states,
the optimization problem need not to be performed since system indices can be evaluated
77

from the failure states directly.
Some failure states may not be captured by using the state space reduction technique
which denoted as unclassified subspace. We have used binary particle swarm optimization
search method to search for such states. The problem of deciding whether we need to perform
binary particle swarm optimization search method for the unclassified subspace can be solved
by running the binary particle swarm optimization search method for few iterations. If the
reliability index under consideration does not change, the algorithm is terminated and the
indices calculated by the line capacity model will be enough to evaluate the reliability of the
system.

6.1.1

Reduction in Computational Time in Case of Ignoring the
Unclassified Subspace

If the transmission lines of the system under consideration are very reliable and have high
power limits in comparison with system loading conditions which is the case for most power
systems, the proposed state space reduction technique can capture most of the failure states.
In this case, performing the binary particle swarm optimization search method is not necessary and ignoring this subspace will not cause significant errors.
If we assume the total time required to evaluate system indices by passing every single
state to the solver (optimization problem) is T , the total time required to evaluate system
indices by passing only the failure states to the solver is Tr , the time required to evaluate
every system state is t, the time required by the optimization problem to evaluate every
system state is tr , number of samples is N and the loss of load probability of the system
is LOLP, then by ignoring the unclassified subspace, the percentage of reduction in the

78

computational time η can be calculated as,

η=

T − Tr
× 100 %
T

(6.1)

where T = N × t and

Tr ≈ N × [LOLP × t + (1 − LOLP ) (t − tr )]

here Tr is approximated because the time of performing the state space reduction technique
to identify the failure states is very small in comparison with the time of preparing system
data for every sampled state and the time required by solving the optimization problem.
After some manipulations, equation (6.1) can be rewritten as,

η ≈ (1 − LOLP )

tr
× 100 %
t

(6.2)

As an example, if we assume 100,000 states were sampled and evaluated in order to reach
the desired tolerance and every state takes about 0.01 seconds and if we assume the LOLP
of the system is 0.02 and the optimization solver takes 95% of the calculation time of each
state, then the time required with solving an optimization problem for “every state” is 1000
seconds (100000 × 0.01) and the time required by passing only the failure states to the solver
is 69 seconds,

Tr ≈ 100000 × [0.02 × 0.01 + (1 − 0.02) (0.01 − 0.95 × 0.01)] = 69 (seconds)

79

Therefore, the percentage of time reduction of using this technique is 93.10%,

η ≈ (1 − 0.02)

6.1.2

0.95 × 0.01
0.01

× 100 % = 93.1%

Reduction in Computational Time in Case of Including the
Unclassified Subspace

If ignoring the unclassified subspace causes considerable error in evaluating system indices,
the proposed binary particle swarm optimization search method should be used. In this case,
calculation of the percentage of time reduction is not straight forward and system dependent.
Therefore, calculation of this time cannot be determined beforehand.
To calculate the time reduction, another time which is the time spent by performing the
binary particle swarm optimization search method, TBP SO , should be added to the time
spent to perform the line capacity model, Tr . Now, equation (6.1) can be rewritten as,

η=

T − (Tr + TBP SO )
× 100 %
T

(6.3)

The computational efforts required by the algorithm proposed in [13] for the entire state
space is found to be only 19.31% of performing Monte Carlo Simulation with State Space
Pruning [8] on the modified IEEE RTS. Since we are applying a dynamically directed binary
particle swarm optimization search method on a subspace of the entire state space, TBP SO
in our case should be less than 19.31% of T . Also, the number of states to be visited is very
small, hence, the required memory will be small too.

80

6.2

Solution Algorithm

The solution algorithm explains the flow of the procedures of evaluating the reliability indices
and their sensitivities with respect to component parameters and the operating conditions
as well as the technique used to direct and prevent the particles from being trapped in a
search subspace. The steps can be explained as follows:
1. Initialize the positions and velocities of the particles, xi and vi respectively. Particles positions are initialized by using the forced outage rates of system components;
components that are in the up state are represented by 1’s and components in the
down state are represented in 0’s. The length of a particle string equals the number of
components. Velocities are initialized to have Vmax equals 2.2 (probability limits are
between 0.9 and 0.1).
2. Check if there are identical particles. If so, discard the identical ones and save the
rest of the particles in a temporary array vector by converting the binary numbers to
decimal numbers. Check if there are particles already exist in the database, if so set
probabilities and load curtailments of the existing particles to zeros to decrease the
chance of visiting these states again. Then, save the rest in the database and go to the
next step.
3. Compute the exact probability of each particle which represents the probability of a
system state. If the probability of a particle is less than a threshold ε discard this
particle, otherwise go to the next step. In this work, ε is chosen to be 10−15 .
4. Set system parameters and update generation and transmission lines status for every
particle. Perform the linear programming optimization problem to check if there are
81

load curtailments. Update the reliability indices.
5. Determine and update personal best states’ probabilities and best load curtailments.
Also, determine and update best global probability and best global load curtailment.
Update particles velocities using (2.16) and update particles positions using (4.7) and
(4.10).
6. Check the normalized LOLP (LOLPnorm ) whether it is increasing or decreasing and
adjust c1 , c2 , c3 and c4 accordingly. If it increasing, increase the values of c1 and c2 and
decrease the values of c3 and c4 . Perform the opposite if LOLPnorm is decreasing. In
this work, we adjust the amount of increasing or decreasing these coefficients according
to the amount of change of the LOLPnorm . Any strategy can be used to adjust these
factors. However, the optimum adjustment of these factors is out of the scope of this
work. We relate the amount of change of these coefficients to the amount of change in
LOLPnorm , (∆LOLPnorm ), in an exponential form as follows,

c(1,2)−new = c(1,2)−old × e(100×∆LOLPnorm )

c(3,4)−new = c(3,4)−old × e(−100×∆LOLPnorm )

7. Check if any of the particles has entered the success and/or failure subspaces. If a
particle has entered the success subspace, redirect the particle by using the objective
of maximum load curtailment only (ignoring the objective of maximum probability)
by setting c1 and c2 to zeros. If a particle has entered the failure subspace, redirect the
particle by using the objective of maximum probability only (ignoring the objective of
maximum load curtailment) by setting c3 and c4 to zeros.
82

8. Check for convergence. If the stopping criterion is met, stop; otherwise go to the next
step.
9. Determine the number of newly visited states since the previous few iterations. If the
number of the new states is less than a threshold, adjust Vmax and then go to step 2;
otherwise go directly to step 2. The threshold is chosen to be 80% of the total number
of particles for the first 10 iterations, 50% for the next 10 iterations and 20% thereafter.
Vmax is adjusted according to the direction of the LOLPnorm index. If the particles
are searching around the success subspace, adjust Vmax to a large number. On the
other hand, if the particles are searching around the failure subspace, adjust Vmax to
a small number. In this work we have chosen 4.0 to be the large number and 2.0 to
the small number. Again, the optimum value of Vmax to better improve the behavior
of the search method is out of the scope of this work.

6.3

Case Studies

The proposed formulation was applied on four systems, the IEEE RTS [63], the Modified
IEEE RTS, and the Saskatchewan Power Corporation in Canada (SPC) and a hypothetical
six bus system [64]. More details about these system are given in the appendix. The
case studies are divided into two cases: (I) calculating the reliability indices, validating
the results and estimating the reduction in the computational time and (II) performing
sensitivity analyses for the reliability indices.

83

6.4

Case Study I: Calculation of the Reliability Indices

The model was applied on the IEEE RTS and the Modified IEEE RTS with 50 particles.
Analysis of the annual as well as the annualized indices was conducted on the two systems.
From Table 6.1, Table 6.2, Table 6.3 and Table 6.4, it has been found that in calculating
the annual indices of both systems and the annualized indices of the original IEEE RTS,
the heuristic technique was able to evaluate system indices and no BPSO was performed
(terminated after few iterations). On the other hand, in calculating the annualized indices
of the modified IEEE RTS, the heuristic technique was not able to evaluate system indices
and therefore the proposed BPSO was used to search for the unclassified subspace. Also, from
Table 6.4, the unclassified subspace has a significant effect on the LOLP and LOEE indices
and very small effect on the EDNS and LOLE indices. These effects can be interpreted as the
unclassified subspace has failure states with relatively high probability but with relatively
small load curtailments. This also can be attributed to the fact that failure states that are
located far from the origin (all components are in the down state) have high probability and
small load curtailments. Therefore, if we were looking for just the EDNS and LOLE indices,
the heuristic technique would be enough for this case.
Table 6.5 and Table 6.6, show the reduction in the computational time by using the
proposed algorithm as a percentage of the time of evaluating the same indices using Monte
Carlo simulation. These times account for calculating the bus indices for all cases. The
average base time using Monte Carlo simulation was found to be 1175 seconds for both
systems.
Figure 6.1 shows the profiles of the LOLP, DLOLP and the projected LOLP as the
algorithm progresses. As the number of visited states increases, the estimated LOLP and

84

Table 6.1: Annual Indices of the IEEE RTS
LOLP
Algorithm

LOEE

LOLE

MW/yr MWh/yr
Heuristic 0.00117

Proposed

EDNS

BPSO
Sum

Conventional

h/yr

0.16374

1430.465

10.22112

–

–

–

–

0.00117

0.16374

1430.465

10.22112

0.00119

0.16412

1433.752

10.39584

Table 6.2: Annualized Indices of the IEEE RTS
LOLP
Algorithm

EDNS

LOEE

LOLE

MW/yr

MWh/yr

h/yr

Heuristic 0.08429 14.67784 128225.60 736.360
Proposed

BPSO
Sum

Conventional

–

–

–

–

0.08429 14.67784 128225.60 736.360
0.08451 14.59761 127524.72 738.279

Table 6.3: Annual Indices of the Modified IEEE RTS
LOLP
Algorithm
Heuristic 0.00106
Proposed

BPSO
Sum

Monte Carlo

EDNS

LOEE

LOLE

MW/yr

MWh/yr

h/yr

0.07395

645.90125 9.26016

–

–

–

–

0.00106

0.07395

645.90125 9.26016

0.00105

0.07436

649.61330 9.17280

Table 6.4: Annualized Indices of the Modified IEEE RTS
LOLP
Algorithm

EDNS

LOEE

LOLE

MW/yr

MWh/yr

h/yr

Heuristic 0.06039 10.31755 90134.083 517.259
Proposed

BPSO

0.00976

Sum

0.07015 10.43926 91197.375 612.830

Monte Carlo

0.12171

1063.256

85.263

0.07017 10.44237 91224.549 613.005
85

Table 6.5: Reduction in the computational time of the IEEE RTS
Time Spent by
Category

Heuristic tech.
(seconds)

Total

Time

BPSO
Time
Reduction
(seconds) (seconds)
(%)

Annual

28.66

0

28.66

97.56

Annualised

144.28

0

144.28

87.72

Table 6.6: Reduction in the computational time of the Modified IEEE RTS
Time Spent by
Category

Heuristic tech.
(seconds)

Total

Time

BPSO
Time
Reduction
(seconds) (seconds)
(%)

Annual

55.60

0

55.60

95.27

Annualised

336.15

409.75

745.90

36.52

the estimated DLOLP becomes close to each other. If the simulation were to continue to
make the particles visit the entire state space which is impractical, these two indices (LOLP
and DLOLP) would converge to the same value. Also, Figure 6.1 shows some changes in
the LOLPproj index specially at the early stages of the searching process. These changes
are due to the adjustments of the weighting factors of the objective function as well as the
adjustments in the velocity limit, Vmax .

6.5

Case Study II: Sensitivity Analyses

The proposed formulation was applied on the IEEE RTS [63]. From Table 6.8, it is clear that
the generating unit (400 MW) at bus 18 has the highest effect on all the studied reliability
indices followed by the generating units of buses 21 and 24 for all indices. Consequently,
these generating units have the largest impact on system reliability. In other words, these
buses are at the highest risk and the most beneficial way to improve the reliability of the

86

Table 6.7: Bus Annual Indices for the IEEE RTS
Bus
No.

LOLP

EDNS
MW

LOEE
MWh/yr

LOLE
hr/yr

1

0.00032 0.00555

48.50943

2.75184

2

0.00002 0.00034

2.929904

0.17472

3

0.00006 0.00188

16.41914

0.52416

4

0.00044 0.00530

46.26799

3.80016

5

0.00124 0.01473 128.67821 10.78896

6

0.00042 0.00423

36.93749

3.66912

7

0.00021 0.00398

34.78341

1.79088

8

0.00013 0.00308

26.94283

1.09200

9

0.00005 0.00125

10.90315

0.43680

10

0.00008 0.00238

20.78260

0.69888

13

0.00028 0.01001

87.42760

2.40240

14

0.00096 0.02881 251.65359

8.34288

15

0.00070 0.03126 273.12752

6.11520

16

0.00110 0.01813 158.37621

9.60960

18

0.00018 0.00862

75.30948

1.52880

19

0.00012 0.00329

28.73840

1.04832

20

0.00037 0.00744

64.97050

3.18864

system is to increase the availability of these units by adding new units, decreasing repair
time and/or using sophisticated methods to decrease units’ failure rate.

6.6

Summary

The heuristic technique, which is presented in chapter 4, was applied to reduce the search
space. The heuristic technique is conservative in that not all the failure states can be captured
and no success state was classified as failure state. The heuristic technique produces three
subspaces, definite failure subspace, definite success subspace and unclassified subspace. An
87

0.8
DLOLP
LOLP
Projected LOLP

0.6
0.4

Probability

0.2
0.0
1.0
Due to the adjustments
in particles direction

0.8
0.6
0.4
0.2
0.0
0

2

4
6
8
Number of Visited States ( x103)

10

12

Figure 6.1: Reliability indices profile of the IEEE RTS System

Table 6.8: Sensitivity Analysis of the LOLP index
Bus No. & ∂LOLP/∂ui

∂LOLP/∂λi

∂LOLP/∂µi

Gen. Cap.

×10−2

×10−4

×10−4

#18G400

33.29905

44.15546

-6.02120

#21G400

33.18920

44.00979

-6.00134

#23G350

26.16603

25.28188

-2.19842

#13G197

20.02105

10.31336

-0.54281

#23G155

11.34010

4.77216

-0.19884

#15G155

11.01719

4.63627

-0.19318

#16G155

10.77760

4.53545

-0.18898

#22G100

4.94948

2.60356

-0.10848

#2G76

4.48367

1.96626

-0.04013

#1G76

4.36122

1.91257

-0.03903

88

Table 6.9: Sensitivity Analysis of the LOLF index
Bus No. &

∂LOLF/∂ui

∂LOLF/∂λi

∂LOLF/∂µi

×10−2

×10−4

Gen. Cap.
#18G400

58.64305

77.76229

-660.80623

#21G400

57.52415

76.27860

-641.89222

#23G350

49.01718

47.36089

-202.50560

#13G197

43.04602

22.17411

-16.60056

#23G155

41.51697

17.47125

-27.43647

#15G155

39.63948

16.68116

-25.43609

#16G155

37.98464

15.98477

-23.49279

#22G100

16.69121

8.78003

-16.78555

#2G76

17.06241

7.48253

-6.30312

#1G76

16.41083

7.19679

-5.96487

Table 6.10: Sensitivity Analysis of the EDNS index
Bus No. & ∂EDN S/∂ui

∂EDN S/∂λi

∂EDN S/∂µi

×10−2

×10−4

Gen. Cap.
#18G400

73.95452

98.06572

-1337.25977

#21G400

73.78475

97.84059

-1334.18991

#23G350

68.91659

66.58790

-579.02520

#13G197

29.80471

15.35317

-80.80615

#23G155

20.02147

8.42547

-35.10614

#15G155

18.68813

7.86438

-32.76824

#16G155

18.83241

7.92509

-33.02120

#22G100

11.23545

5.91016

-24.62565

#2G76

9.13458

4.00587

-8.17524

#1G76

8.76626

3.84434

-7.84560

89

intelligent search method based on binary particle swarm optimization was used to search
for failure states in the unclassified subspace. From the failure states which were captured
by the heuristic technique and the binary particle swarm optimization, the reliability indices
and their sensitivities with respect component parameters were calculated.
An element key in using binary particle swarm optimization in composite system reliability evaluation is the selection of weighting factors of the objective function. These
weighting factors are system dependent and, therefore, their appropriate values should be
carefully selected, in order to prevent the swarm from being trapped to one corner of the
state space. An effective particle swarm optimization technique has been proposed and used
for reliability assessment of composite systems. The proposed method adjusts the weighting
factors associated with the objective function in a dynamic fashion so that the swarm would
always fly on the entire search space rather of being trapped to one corner of the search
space. The effectiveness of the proposed method was demonstrated on four well-known test
system. The results show that the reliability indices obtained by the proposed method correspond closely with those obtained using Monte-Carlo simulation, while requiring lower
computational burden.

90

Chapter 7
Reliability and Sensitivity Analysis of
Composite Power Systems Under
Emission Constraints

7.1

Introduction

In view of the increasing role of environmental considerations in power generation, there is
an evolving body of minimum cost and minimum emission methods that are cognizant of
these emerging factors. This chapter presents a method of power system expansion planning,
considering emission constraints. The method presented here uses a linearized power flow
representation in the form of a linear programming (LP) model similar to those frequently
used in power system studies to minimize the total operation cost and is applied to both
generation and transmission expansion planning. The cost of the fuel, the cost of purchasing
additional emission allowances in case a power plant reaches its maximum emission allowance
and the gain from selling the excess emission credits where the emissions are less than the
limit are combined to develop a piecewise linear objective function. Emissions are represented
as a cost in case of buying extra emission allowances or as benefits in case of selling the excess
emission allowances. It also utilizes the concept of the shadow price to perform sensitivity

91

analysis with respect to the operation constraints to identify the best cost effective generation
and transmission expansion planning. The dual solution of the LP provides the shadow
prices that are used for determining the sensitivity of the minimum cost with respect to
power generation.
Following the introduction of the Clean Air Act Amendments of 1990, utilities have been
required to reduce and monitor emissions from their power generating units starting from
1995. This act enforced allowances on SO2 emissions that could be used or traded within the
specified banking period [65]. The most significant reductions have occurred very recently.
On November 14, 2012, California had its first auction and trading on greenhouse emission
allowances [66]. Starting January 2013, California’s emission cap-and-trade program came
into effect, and greenhouse gas emissions from large electric power plants are required to
be reduced by 16% between 2013 and 2020 [66]. The Regional Greenhouse Gas Initiative
(RGGI), which consists of nine North-Eastern and Mid-Atlantic states, agreed on reduction
of CO2 budget by 45% by 2014 and the reduction will decrease by 2.5% yearly from 2015 to
2020 [67]. The American Clean Energy and Security Act of 2009 provides for the reduction of
carbon emissions by 17% by 2020. The European Union Emissions Trading System reduced
greenhouse gas emissions more than 8% between 2008 and 2012 and is expected to reduce
greenhouse gas emissions below 20% of the 1990 levels by 2020, starting from 2013 [68].
Several methods have been proposed to reduce emissions from power generating units;
these include installing post-combustion cleaning equipment, switching to fuel with less pollutants, increasing the penetration of the renewable energy, and modifying existing dispatch
strategies to include emissions [69]. With these reduction methods, utilities are to operate
and expand their generation within the emission allowances.
Some of the earliest work on minimum emission generation dispatch was performed by
92

Gent and Lamont [70]. Nada, Hari and Kothari [71] solved economic emission load dispatch with line capacity constraints. Subsequently, considerable research was reported on
variations and extensions of these methods [72–78]. These variations use different means
of accommodating emissions within the optimization problem. Some of the proposed formulations also take into account the different emission reduction equipment and methods
described above, and thereby provide for strategies that may allow utilities to forgo, defer,
or minimize additional capital costs. For example, [79] shows how a utility may avoid installation of new emission equipment by changing its commitment and dispatch schedules,
or by switching to fuels with low emission (and high cost), or both.
Reducing emissions not only affects power dispatch and cost-effective planning, but also
influences system reliability requirements. Imposing emission limits on some specific areas may change reliability-based expansion planning procedures. Planning decisions may
comprise a trade-off between the cost of buying extra emission allowances and the cost of
improving system reliability. Evaluating system reliability under the condition of reducing
the emission became necessary and reliability indices need to be adapted to this new constraint. Therefore, there is an emerging need to include emission considerations in reliability,
security and economy studies.
Sensitivity analyses of conventional reliability indices incorporate remedial action through
load curtailment minimization under the constraints of power balance, generation capacity
limits and transmission capabilities to rank system components according to their risk. However, in accommodating regulatory measures, emission allowances have been considered as
additional constraints. While planning studies are concerned with reliability and cost, it
is often useful to quantify the sensitivity of the reliability indices of interest to component
parameters and operating constraints. In the emerging context, it is reasonable to include
93

emission constraints in the sensitivity analysis of the reliability indices. In other words,
sensitivity of reliability indices with respect to the emission limits needs to be added to expansion planning studies. Therefore, generation expansion planning is no longer based only
on reliability enhancement or cost minimization.
A considerable amount of research has been conducted on economic/emission dispatching
strategies. However, incorporating emission allowances in power system reliability studies
is rarely found in the literature. In [80], SO2 , N Ox and CO2 emissions were considered as
constraints in their generation expansion model which is based on multi-area reliability exponential analytic model. Minimization of emissions while maintaining reliability requirements
in microgrids is presented in [81]. Most of the prior research has focused on evolutionary
techniques and multiobjective minimization function to minimize emissions and generation
costs in generation expansion planning [72–76, 82].
The work reported in this chapter calculates the reliability indices of composite power
systems and utilizes the concept of the shadow price in ranking the generating units from the
reliability point of view while considering the emission allowances as additional constraints.
Therefore, this work adds another dimension to system planning processes through the inclusion of emission constraints. Sensitivity analyses of the well-known reliability indices, LOLP,
LOLF and EDNS, are presented. Also, a new expression for calculating the sensitivity of
the EDNS with respect to the emission limits is developed. The generating units are ranked
based on their sensitivities to the emission allowances. Based on this ranking, the emission
allowances of generating stations are evaluated to determine a need for expansion. For instance, if a generating plant is ranked high in terms of risk and is viewed as a candidate
for reconditioning in order to increase the overall system reliability, the emission allowance
should be applied. If the emission limit is reached, the planner should choose between buy94

ing emission allowances and investing in another power station that has available emission
allowances. The method is applicable to any source that produces undesirable byproducts
that are capped by regulatory and other policies, as long as these byproducts are a function
of the output power. The method is demonstrated on the IEEE Reliability Test System
(IEEE-RTS).

7.2

Emissions from Electric Generation

A significant amount of energy supply is used in producing electricity worldwide, and in
most industrialized countries the bulk of electric generation comes from fossil fueled plants.
According to the Environmental Protection Agency (EPA), about 16 million tons of SO2
and 7 million tons of N Ox were emitted by electric utilities in the United States in 1985 [83].
Further, EPA reports that in 1995 about 23 million tons of SO2 and 19 million tons of N Ox
were emitted from all resources. In the UE-27, 40% of energy resources is used in generating
electricity [84]. Also, in the UE-27 [84], 55% of the generated electricity in 2005 was supplied
by fossil fuel resources. In 2008, 601.32 GW was generated from coal in China which is about
75.87% of the total electric energy generated [85]. Most of these emissions were produced
by fossil fuel plants. In 2006, 12 million tons of Sulphur-dioxide emissions were produced in
China [80].
Some regulations were introduced and categorized according to some definitions such
as, “new source performance”, “modified source”, “reconstruction”, “degree of emission
limitation”, etc. For instance, “new source” is defined as any source that is commissioned
after the regulation is proposed and is subjected to the NSPS (New Source Performance
Standards) requirements even if the regulation is not final [83]. Also, a “modified source”

95

is any source subjected to a physical or operation change that may increase air pollutant or
produce a new air pollutant and is subject to the NSPS (CAA).
The goal of reducing the amount of emissions of SO2 and N Ox is implemented in two
stages: Phase I and Phase II. Phase I began in January 1, 1995 and Phase II began in
January 1, 2000. For example, in phase I, all units with 100 MW or greater and had SO2
emission greater than 2.5 lb/mmBtu in 1985 are subjected to emission reduction according
to table A of section 404, title IV of the CAA.

7.3

CO2, SO2 and N Ox Emission Constraints

This section describes the formulation and incorporation of the emission in the load curtailment minimization problem (chapter 2). In this section, the emission allowances are
considered as additional constraints.
Several methods have been proposed to model the emission function of the thermal
generating units. Some of these models are based on the heat rate functions of the generating
units with modifying the heat rate coefficients. Some other model emissions based on the
heat rates with multiplying the heat rates with appropriate emission factors. In addition,
some models use more detailed models such as including the valve loading effect. The heat
rates themselves can be modeled in a quadratic, cubic or exponential function. Emission
rates were modeled using a straight linear form as in [77,80,86], combination of a straight line
and two exponentials as in [70], a quadratic form as in [71–75, 78, 79, 82, 87, 88], a quadratic
and an exponential form as in [76, 89, 90], and multiplying the heat rate function with an
emission factor as in [91]. In this work, for simplicity, the emissions are modeled based on
the heat rates, while ignoring the valve loading effects, and linearizing the quadratic form.

96

The quadratic heat rate for a generating unit i can be expressed as:

HR = ai + bi Pi (t) + ci Pi2 (t) (MMBtu/MWh)

(7.1)

where HR is the heat rate, ai , bi and ci are the heat rate coefficients and Pi (t) is the real
power generation of unit i at time t. Another widely used unit for heat rates is Btu/kWh.
MMBtu/MWh has been used so as to be consistent with units of the emission rates. To
convert from MMBtu/MWh to Btu/kWh, the given heat rate can be multiplied by 1,000.
Since the heat rates are functions of the output power and the emissions are functions
of the heat rates, the emissions of a generating unit i can be expressed as a function of the
output power as follows,

Eji = αji + βji Pi (t) + γji Pi2 (t) (lbs/MMBtu)

(7.2)

where α, β and γ are the coefficients of the emission rates, j denotes emission type (1 for
ECO2 , 2 for ESO2 and 3 for EN Ox ), ECO2i is the CO2i emission of unit i, ESO2i is the
SO2i emission of unit i and EN Oxi is the N Oxi emission of unit i.
Assuming a power generating plant with n thermal units at time t, the total CO2 emission
is expressed as
n

ECO2 =

ECO2i (Pi (t))

(ton/h)

(7.3)

(ton/h)

(7.4)

i=1

Similarly, the total SO2 emission is expressed as
n

ESO2 =

ESO2i (Pi (t))
i=1

97

Also, the total N Ox emission is expressed as
n

EN Ox =

EN Oxi (Pi (t))

(ton/h)

(7.5)

i=1

These emissions are nonlinear functions in the unit output power. In order to include
these emissions as constraints in the linear programming problem, these equations are linearized around the operating points.

7.4

Sensitivity Analysis of EDNS Index with respect
to the Emission Limits

The optimization is constrained by generator and transmission line capacities. However,
after the introduction of the Clean Air Act, it is necessary for such analyses to accommodate
emission limits. This work includes the effect of emission limits on the sensitivity analysis
of the reliability indices. The most commonly used indices in composite power system are:
LOLP, LOLF and EDNS. The sensitivity of the LOLP, LOLF and EDNS with respect to
component parameters has been introduced in chapter 3. In this section, an analytical
expression is developed to calculate the sensitivity of the EDNS with respect to emission
limits. The development of this expression is based on the linearized relationship between
the amount of emission and the output power. The derivation starts with utilizing the partial
derivative of the EDNS to the component capacity.
The sensitivity analysis of the EDNS with respect component capacities was derived as
follows [32, 33].
∂EDN S/∂Ci =

If (x) P (x) ∂Lc (x) /∂Ci
x∈X

98

(7.6)

where Ci is the capacity of component i and Lc (x) is the total load curtailment when the
system is at state x.
The derivative of the total load curtailment with respect to component capacity is given
in section 3.4.2. Sensitivity of EDNS with respect to circuit capacity or susceptance is
inherently not affected by the emission limits; however, with generation redispatch under
the condition of emission limits reduction they may become active constraints.
In section 7.3 it was stated that the emission functions are linearized about the operating
points. Therefore, the effect of the emission limits on EDNS sensitivity can be evaluated by
developing a relationship between shadow prices of generating units with its emission. From
(2.16), the Lagrange multiplier of a generating unit can be expressed as

πg,i = ∂Lc (x)/∂gi

(7.7)

where gi is the capacity of generator i.
Now a new term is introduced which is the partial derivative of total load curtailment
with respect to emission limit. This term can be thought of it as a new constraint to (2.16).
This new constraint will produce a new Lagrange multiplier which relates the emission limit
to total load curtailment. The partial derivative of total load curtailment with respect to
emission limits can be expressed as

πE,ij = ∂Lc (x)/∂Eij

(7.8)

where πE,ij and Eij are the Lagrange multiplier or shadow price of the emission limit and

99

the emission limit of unit i for emission type j. ∂Lc (x)/∂Eij can be expressed as

∂Lc (x)/∂Eij = ∂Lc (x)/∂gi · ∂gi /∂Eij

(7.9)

The first term of the right hand side of (7.9) is the shadow price of the generating capacity.
The second term can be calculated using the relationship of generation output power and
emission. This relationship can be given as,

Eij = Ec,ij · gi

(7.10)

where Ec,ij is the linearized emission coefficient of unit i for emission type j.
Therefore, from (7.10), the partial derivative of the output power with respect to the
amount of emission can be calculated as follows.

∂gi /∂Eij = 1/Ec,ij

(7.11)

Now, the partial derivative of the total load curtailment with respect to the emission
allowances can be given by

∂Lc (x)/∂Eij = ∂Lc (x)/∂gi · 1/Ec,ij

(7.12)

This can be expressed in terms of shadow price or Lagrange multiplier as

πE,ij = πg,i /Ec,ij .

100

(7.13)

Therefore, the sensitivity of the EDNS with respect to emission allowances now can be
defined as,
∂EDN S/∂Eij =

If (x)P (x)∂Lc (x)/∂Eij

(7.14)

x∈X

Expressed in terms of shadow price or Lagrange multiplier,

∂EDN S/∂Eij =

If (x)P (x) πg,i /Ec,ij

(7.15)

x∈X

Equations (7.14) and (7.15) calculate the sensitivity of the EDNS with respect to emission
allowances. Therefore, by solving the optimization problem of (2.16) and getting the shadow
prices of the generating units, the shadow prices of the emission limits are calculated by
dividing the shadow prices of the generating units by the linearized emission coefficients.

7.5

Solution Algorithm

This section describes the approach used for calculating the sensitivity analysis of LOLP,
LOLF and EDNS with respect to component reliability and emission allowance limits of
various system components. As mentioned above, emission allowances are considered as
additional constraints. The sensitivity analysis was carried out using the heuristic technique
and the directed binary particle swarm optimization search method as described in chapter
4.
The stopping criterion used in this method consists of comparing the coefficient of variation of a selected index (e.g., the EDNS) to a pre-specified tolerance. The simulation steps
are summarized in Figure 7.1.

101

Start
Read system data

Sample system state
Solve the LP problem

Curtailment = 0 ?
No
Determine shadow prices
and update indices

Simulation
converged?

No

Yes
Output results and stop
Figure 7.1: Flowchart of Estimating the Sensitivity Analysis

7.6

Case Studies

This section illustrates the method on the IEEE RTS test system [63]. The data of this
system is given in the appendix. The IEEE RTS has two nuclear steam generating units
each of which 400 MW and six hydro driven generating units each of which 50 MW. These
units do not have SO2 , N Ox or CO2 emissions. The emission rates were calculated from
the heat rates. The emission coefficients were calculated from the emission rates based on
2nd order polynomials.

102

For every state, the partial derivative of the total load curtailment with respect to generation capacity is calculated from the optimization problem of (2.16). The emission curves
were linearized about the operating points for two linear segments. After the simulation
converges to the specified accuracy, the sensitivity of the EDNS to the emission limits is
summed up as shown in (7.15).
Two case studies have been performed: a base case, and one in which the emission
allowances are reduced by 5%. The annualized reliability indices of the IEEE RTS are
shown in Table 7.1. The results for two more reliability indices, Loss of Load Expectation
(LOLE) and Loss of Energy Expectation (LOEE, also called expected un-served energy or
EUE), are also presented. LOLE and LOEE can be calculated directly from LOLP and
EDNS respectively by multiplying by the study period in hours.
The sensitivity results of the LOLP with respect to reliability parameters for the base
case and reduced emission case are shown in Table 7.2 and Table 7.3, respectively. The
sensitivity results of the LOLF with respect to reliability parameters for the base case and
reduced emission case are shown in Table 7.4 and Table 7.5, respectively. The sensitivity
results of the EDNS with respect to reliability parameters for the base case and reduced
emission case are shown in Table 7.6 and Table 7.7, respectively. Table 7.8 shows the
sensitivity of the EDNS with respect to ESO2 , EN Ox and ECO2 for the base case. Table
7.9, shows the sensitivity of the EDNS with respect to ESO2 , EN Ox and ECO2 for the
reduced emission case.

103

7.7

Results and Discussion

Is is clear that the generating unit (350 MW) at bus 23 has the highest effect on the LOLP
followed by the generating units (400 MW) of buses 18 and 21 for the two study cases.
Consequently, these generating units have the largest impact on system reliability. In other
words, these buses are at the highest risk and the most beneficial way to make the system
more reliable is to increase the availability of these stations by adding new units, or by
modernizing existing units to improve reliability. When the emission allowances were reduced
(case study 2), these sensitivities increased accordingly.
In comparison to the case study 1, the sensitivities of LOLP of the generating units of
the case study 2 have increased due to the decrease in the emission limits. Such decrease in
the emission limits has forced the generating units to have larger shadow prices values than
it was in case study 1. By decreasing the emission limits, the emission limits constraints
became binding constraints and generating capacities became redundant constraints.
Also, it is clear that the generating unit (350 MW) at bus 23 has the highest effect on
the EDNS for both study cases. Also, after reducing the emissions, the sensitivities of the
EDNS with respect to the emission limits have increased dramatically. For instance, for the
SO2 air pollutant, sensitivity of the 350 MW unit with respect to the emission has increased
by about 2.6 times and for the 197 MW unit has increased by about 4.9 times. The factor
of increase in the sensitivity of the EDNS with respect to emission limit due reducing the
emissions are shown in Table 7.10.
From the above discussion it is obvious that reducing the emission allowance limits leads
to a huge change in system reliability indices. From these studies, planners can decide
where to invest by trading off between reliability requirements and the cost of reducing the

104

emissions or buying emission allowances for the affected areas.
Table 7.1: Annualized Indices of the IEEE Reliability Test System
LOLP

EDNS

LOEE

LOLE

MW

MWh/yr

hr/yr

Base Case

0.08496

14.76607

128996

742.21056

Reduced Emission

0.14093

26.49218

231435

1231.16448

Table 7.2: Sensitivity Analysis of LOLP with Respect to Reliability Parameters for the Base
Case
Bus No. & ∂LOLP/∂ui

7.8

∂LOLP/∂λi

∂LOLP/∂µi

Gen. Cap.

×10−4

×10−2

×10−4

#15G12

0.23468

0.13523

−0.27598

#22G50

1.57580

0.30889

−0.31201

#16G155

14.87768

5.48451

−22.85212

#23G155

13.65057

5.03215

−20.96728

#15G155

14.49772

5.34444

−22.26849

#13G197

15.36175

6.93199

−36.48416

#23G350

42.33958

35.83622

−311.61931

#18G400

30.31500

35.21391

−480.18962

#21G400

30.42285

35.33918

−481.89795

Summary

Most of the emission studies in the literature focused on minimizing emission with economic
dispatch and introducing new concepts to solve multi-objective functions. The goal was
to minimize emissions and fuel costs to accommodate for the requirements of the CAA of
1990. Minimizing the emission is one solution among many methods such as utilizing fuels
with less air pollutant, exploring new techniques to reduce the amount of emissions, etc.
105

Table 7.3: Sensitivity Analysis of LOLP with Respect to Reliability Parameters for the
Reduced Emission Case
Bus No. & ∂LOLP/∂ui

∂LOLP/∂λi

∂LOLP/∂µi

Gen. Cap.

×10−4

×10−2

×10−4

#15G12

1.39303

0.80272

−1.63820

#22G50

4.97198

0.97461

−0.98445

#16G155

17.80651

6.56419

−27.35080

#23G155

16.73514

6.16924

−25.70517

#15G155

18.46648

6.80748

−28.36451

#13G197

19.05290

8.59762

−45.25063

#23G350

46.66679

39.49877

−343.46757

#18G400

33.29836

38.67937

−527.44598

#21G400

32.84112

38.14824

−520.20328

Table 7.4: Sensitivity Analysis of LOLF with Respect to Reliability Parameters for the Base
Case
Bus No. &

∂LOLF/∂ui

∂LOLF/∂λi

∂LOLF/∂µi

Gen. Cap.

×10−4

×10−2

×10−4

#15G12

0.23468

0.13523

−0.27598

#22G50

1.57580

0.30889

−0.31201

#16G155

14.87768

5.48451

−22.85212

#23G155

13.65057

5.03215

−20.96728

#15G155

14.49772

5.34444

−22.26849

#13G197

15.36175

6.93199

−36.48416

#23G350

42.33958

35.83622

−311.61931

#18G400

30.31500

35.21391

−480.18962

#21G400

30.42285

35.33918

−481.89795

106

Table 7.5: Sensitivity Analysis of LOLF with Respect to Reliability Parameters for the
Reduced Emission Case
Bus No. &

∂LOLF/∂ui

∂LOLF/∂λi

∂LOLF/∂µi

Gen. Cap.

×10−4

×10−2

×10−4

#15G12

1.39303

0.80272

−1.63820

#22G50

4.97198

0.97461

−0.98445

#16G155

17.80651

6.56419

−27.35080

#23G155

16.73514

6.16924

−25.70517

#15G155

18.46648

6.80748

−28.36451

#13G197

19.05290

8.59762

−45.25063

#23G350

46.66679

39.49877

−343.46757

#18G400

33.29836

38.67937

−527.44598

#21G400

32.84112

38.14824

−520.20328

Table 7.6: Sensitivity Analysis of EDNS with Respect to Reliability Parameters for the Base
Case
Bus No. & ∂EDN S/∂ui

∂EDN S/∂λi

∂EDN S/∂µi

Gen. Cap.

×10−4

×10−2

×10−4

#15G12

0.23468

0.13523

−0.27598

#22G50

1.57580

0.30889

−0.31201

#16G155

14.87768

5.48451

−22.85212

#23G155

13.65057

5.03215

−20.96728

#15G155

14.49772

5.34444

−22.26849

#13G197

15.36175

6.93199

−36.48416

#23G350

42.33958

35.83622

−311.61931

#18G400

30.31500

35.21391

−480.18962

#21G400

30.42285

35.33918

−481.89795

107

Table 7.7: Sensitivity Analysis of EDNS with Respect to Reliability Parameters for the
Reduced Emission Case
Bus No. & ∂EDN S/∂ui

∂EDN S/∂λi

∂EDN S/∂µi

Gen. Cap.

×10−4

×10−2

×10−4

#15G12

1.39303

0.80272

−1.63820

#22G50

4.97198

0.97461

−0.98445

#16G155

17.80651

6.56419

−27.35080

#23G155

16.73514

6.16924

−25.70517

#15G155

18.46648

6.80748

−28.36451

#13G197

19.05290

8.59762

−45.25063

#23G350

46.66679

39.49877

−343.46757

#18G400

33.29836

38.67937

−527.44598

#21G400

32.84112

38.14824

−520.20328

Table 7.8: Sensitivity Analysis of EDNS with Respect to Emission Allowance Limits for the
Base Case
Bus No. &
Gen. Cap.

∂EDN S/
∂ESO2

∂EDN S/
∂EN Ox

∂EDN S/
∂ECO2

×10−4

×10−4

×10−4

#15G12

−34.95128

−13.98051

−0.04112

#16G155

−94.96082

−37.98433

−0.09044

#23G155

−47.96955

−19.18782

−0.04569

#15G155

−94.63840

−37.85536

−0.09013

#13G197

−47.58625

−19.03450

−0.05598

#23G350

−267.88743 −107.15497

−0.25513

108

Table 7.9: Sensitivity Analysis of EDNS with Respect to Emission Allowance Limits for the
Reduced Emission Case
Bus No. &
Gen. Cap.

∂EDN S/
∂ESO2

∂EDN S/
∂EN Ox

∂EDN S/
∂ECO2

×10−4

×10−4

×10−4

#15G12

−39.70206

−15.88082

−0.04671

#16G155

−274.08317 −109.63327

−0.26103

#23G155

−57.60845

−23.04338

−0.05487

#15G155

−229.14951

−91.65980

−0.21824

#13G197

−233.52268

−93.40907

−0.27473

#23G350

−702.68535 −281.07414

−0.66922

Table 7.10: Percentages of Increase in the Sensitivities of EDNS with Respect to Emission
Allowance Limits Due to Reduced the Emissions
Bus No. & Ratio of Increase
Gen. Cap.
#15G12

1.14

#16G155

2.89

#23G155

1.20

#15G155

2.42

#13G197

4.90

#23G350

2.62

Minimizing the emissions can help in redispatching the existing generators. The need to
establish strategies to include the CAA 1990 in future planning is necessary. Inclusion of
emission constraints in planning studies as described in this chapter helps guide investment
strategies with respect to reconstruction or upgrades.
The work presented here responds to the emerging need for reliability methods that
are cognizant of environmental factors. This work introduced a framework for evaluating
the effect of the emission rates and allowances on system reliability and determining the
sensitivity of commonly used reliability indices to equipment and emission constraints.
109

The Expected Demand Not Supplied (EDNS) is used as a key reliability index. The
sensitivity of the EDNS with respect to emission allowances is evaluated and used to study the
effect of emission allowances on system reliability. From the sensitivity results, components
can be ranked from reliability point of view. Therefore, the component that has the highest
effect on system reliability can be identified. Also, if the reliability of the system needs to
be improved, the sensitivity analyses can identify which component to increase its emission
limit.

110

Chapter 8
Reliability and Sensitivity Analysis of
Composite Power Systems
Considering Voltage and Reactive
Power Constraints

8.1

Introduction

There is a growing need to include the effects of voltage and reactive power constraints in
power system reliability evaluation. Most of the present methods that evaluate the reliability
of composite systems represent the network using the DC power flow model, which ignores
the effects of the voltage and reactive power constraints on the reliability indices. Due to
the integration of the renewable energy resources, the voltage and reactive power limits are
expected to have more effect on the system reliability. The time and computational efforts to
evaluate the reliability indices and their sensitivities with respect to the operating conditions
using models that incorporate voltage and reactive power limits such as AC power flow model
are of concern in both planning and operation prospectives.
Traditionally, composite system reliability methods have used transportation model and

111

dc power flow models to represent system flows, largely due to the complexity of AC flowbased redispatch algorithms. Evaluation of the effects of the voltage and reactive power
limits has been introduced in [92]. The method presented in [92] is based on using the DC
power flow model in two steps and linearizing the relationship between the reactive power
injections and the voltages at the nodes. Several methods have been recently introduced
in the literature of power system reliability to investigate the effects of the lack of reactive
power support on the reliability indices. For example, using the AC power flow model in
two steps to study the aspects of the reactive power on the composite system reliability has
been introduced in [93–96].
Several methods have been introduced in the literature to decrease the time and computational efforts in evaluating power system reliability indices as it was mentioned in chapter
4. This chapter investigates the effects of the voltage and reactive power constraints on composite system reliability indices through performing sensitivity analyses of the Severity Index
with respect to the voltage and reactive power constraints. In this chapter, expressions are
developed to evaluate the sensitivity of the Severity Index with respect to these constraints.
Further, four indices are used to evaluate the contributions of the minimum and maximum
voltage limits at the nodes and minimum and maximum reactive power capabilities in the
loss of load probability index. The full, non-linear AC power flow model, which is given
in (2.14) and (2.15), is used to incorporate these constraints in the reliability evaluation.
The state space pruning technique, which is introduced in chapter 3 is utilized to reduce
the computational time. Also, in this chapter, composite system reliability indices using AC
power flow model are provided to serve as a benchmark for future work involving power flow
models. A comparison between the use of AC and DC power flow models is also provided.
We have used the state space reduction technique of chapter 4.
112

8.1.1

Calculation of Voltage and Reactive Power Limit Indices

Calculations of the LOLP, LOLF, EDNS, LOEE and LOLE are similar to those presented in
chapter 3. This section develops expressions to calculate the voltage and reactive power limit
indices. In this chapter, we use four indices to represent the contributions of violations of the
voltage and reactive power constraints on the loss of load probability index. These indices
are defined as follows: (a) υimin represents the contributions of the minimum voltage level
constraints, (b) υimax represents the contributions of the maximum voltage level constraints,
(c) imin represents the contributions of the minimum reactive power level constraints and
(d) imax represents the contributions of the maximum reactive power level constraints. For
every tested state, if any of these indices is involved in the load curtailment, this index is
updated. These indices are related to the failure subspace not to the system state space;
these indices reflect the contributions of the violations of the voltage and reactive power
constraints on the loss of load probability index. Therefore, the state space of these indices
is Xf .
The estimated values of these indices can be calculated as follows,

fk = E [qk ]

(8.1)

where fk is the index under consideration, (f1 is for the υimin , f2 is for the υimax , f3 is for
the imin and f4 is for the imax ), and qk can be evaluated as follows,
1
qk =
Tf

N

ξki

(8.2)

i=1

where Tf is the sum of the durations of tested failure states and ξki is an indicator function

113

for state violation that can be expressed as,

ξki =





τi , If xi ∈ Xk



0,

(8.3)

otherwise

where Xk is the set in which the index k has a contribution in the load curtailment and
Xk ⊂ Xf .

8.2

Sensitivity Analysis of the EDNS Index with respect to the Voltage and Reactive Power Limit Indices

Sensitivity analyses of composite system reliability indices with respect to system parameters
have been extensively addressed in the literature. To calculate the sensitivity of the EDNS
index with respect to the voltage and reactive power limits; these limits need to be added
to the optimization problem as shown in (2.15).
The sensitivity analyses of the reliability indices with respect to the voltage and reactive
power constraints are important in identifying the impacts of these constraints on power
system reliability and help in hardening power systems against catastrophic failures such as
voltage collapse. Expressions for calculating the sensitivity of the EDNS index with respect
to voltage and reactive power limits are developed. These expressions were developed from
the Lagrange multipliers of the voltage and reactive power constraints.

114

8.2.1

Sensitivity of the EDNS Index with Respect to Voltage
Limit Constraints

The derivation of the sensitivity of the EDNS index with respect to voltage and reactive
power limits starts with utilizing the expression of (3.5). The sensitivity of the EDNS index
with respect to voltage limits allowed at bus k can be expressed as follows,

p (x) × (∂C (x) /∂υk )

∂EDNS/∂υk =

(8.4)

x∈Xf

where ∂C (x) /∂υk is the Lagrange multiplier (πυk ) of the voltage constraint at bus k while
the system is residing at state x.
Also, the sensitivity of the EDNS index with respect to voltage limit constraints can be
evaluated using Monte Carlo state next event method as follows [97, 98],

1
∂EDNS/∂υk =
T

N

ϕi × πυk .

(8.5)

i=1

This expression can be used for both maximum and minimum voltages allowed at the
buses.

8.2.2

Sensitivity of the EDNS Index with Respect to Reactive
Power Limit Constraints

In the same manner, the sensitivity of the EDNS index with respect to the reactive power
available at bus k can be expressed as follows,

115

p (x) × (∂C (x) /∂ k )

∂ρ/∂ k =

(8.6)

x∈Xf

where ∂C (x) /∂ k is the Lagrange multiplier, π k , of the reactive power available at bus k
while the system is residing at state x.
Similarly, the sensitivity of the EDNS index with respect to reactive power constraints
at bus k can be evaluated using Monte Carlo state next event method as follows,

1
∂ ρˆ/∂ k =
T

N

ϕi × π k .

(8.7)

i=1

Also, this expression can be used for both maximum and minimum reactive power limits
at the buses.

8.3

Case Studies

The proposed formulation was applied on the IEEE RTS [63]. The IEEE RTS system has
been extensively tested for power system reliability analysis.

8.3.1

Reliability Assessment

Reliability assessments using both the AC and DC power flow models have been performed
on the IEEE RTS system at the peak load (Annualized Indices). As a comparison, the results
of both models for the IEEE RTS are shown in Table 8.1. As it is obvious from Table 8.1,
estimates from the DC power flow model produce optimistic results in comparison with the
AC power flow model counterpart.
It is well-known that the transmission lines of the IEEE RTS are very reliable with respect
116

Table 8.1: System Annualized Indices of the IEEE RTS
Power Flow

qˆ

Model

ρˆ

φˆ

τˆ

MW

occ./yr

hour

AC

0.10427 17.37479 22.76247 0.00458

DC

0.08455 14.72996 19.76867 0.00428

to the generation. Also, the power carrying capability limits of the transmission lines are
much higher than the normal loading level even in the case of the peak load. Therefore,
the contributions of the transmission lines on the system reliability of this test system are
very small and can be ignored. For this reason, several studies have suggested to use the
modified version of it which is the same as the original system except that the generation
is multiplied by 2 and the load at the buses is multiplied by 1.8. The results of both power
flow models for the Modified IEEE RTS are shown in Table 8.2. Again, from Table 8.2, the
DC power flow model gives very optimistic results in comparison with the AC power flow
model counterpart.
Table 8.2: System Annualized Indices of the Modified IEEE RTS
Power Flow
Model

qˆ

ρˆ

φˆ

τˆ

MW

occ./yr

hour

AC

0.44337 20.13860 63.82885 0.00695

DC

0.07141 10.52347 17.76162 0.00402

Since there are significant differences in the values of the indices for different power flow
models, two case studies have been performed by gradually stressing the transmission lines
to investigate the effects of the voltage and reactive power limits on the reliability indices
and to compare the results of the AC and DC power flow models.

117

8.3.1.1

Case Study I: Gradually Increasing the Generation and Load Levels

In this case study, the generation and load levels are increased gradually to track the behavior
of system indices. The annualized reliability indices are shown in Table 8.3. The first
two columns of Table 8.3, show the factors by which the load and generation are increased
respectively. From Table 8.3, not much change in the system EDNS index between the IEEE
Table 8.3: System Annualized Indices for Case I: Effect of Gradually Increasing the Generation and Load Levels
Loading Generation

qˆ

ρˆ

φˆ

τˆ

MW

occ./yr

hour

Factor

Factor

1.000

1.000

0.10427 17.37479 22.76247 0.00458

1.080

1.100

0.09245 15.73231 20.50596 0.00451

1.160

1.200

0.07459 13.91987 17.03678 0.00438

1.240

1.300

0.06243 13.07068 16.30765 0.00383

1.320

1.400

0.05668 12.99306 14.51361 0.00391

1.400

1.500

0.05551 12.74350 14.54252 0.00382

1.408

1.600

0.06897 12.51458 18.99008 0.00363

1.560

1.700

0.11062 13.18529 25.88540 0.00427

1.640

1.800

0.13719 14.96639 30.11179 0.00456

1.680

1.850

0.15788 15.37179 36.89292 0.00428

1.720

1.900

0.17134 16.72898 43.63865 0.00438

1.740

1.925

0.19511 16.66352 45.35073 0.00430

1.760

1.950

0.23595 17.61181 51.25988 0.00460

1.780

1.975

0.35426 18.56049 59.52425 0.00595

1.800

2.000

0.44337 20.13860 63.82885 0.00695

RST and Modified IEEE RTS. System loss of load probability and loss of load frequency
indices are smoothly increasing with the increase in the loading and generation levels until
the load level reaches around 170% of the peak load. After this point, these indices show
a sharp increase. This sudden increase can be attributed to the voltage and reactive power
118

Table 8.4: System Annualized Indices for Case II: Effect of Gradually Increasing the Load
Levels and Fixing the Generation
Loading Generation

qˆ

ρˆ

φˆ

τˆ

MW

occ./yr

hour

Factor

Factor

1.000

2.000

0.00139

0.00800

0.37248

0.00373

1.080

2.000

0.00171

0.02570

0.68971

0.00248

1.160

2.000

0.00185

0.04468

0.85359

0.00217

1.240

2.000

0.00202

0.07438

0.85251

0.00237

1.320

2.000

0.00201

0.09948

0.82788

0.00243

1.400

2.000

0.00272

0.20431

1.23875

0.00220

1.480

2.000

0.00557

0.40744

2.55331

0.00218

1.560

2.000

0.01041

1.05772

4.68652

0.00222

1.640

2.000

0.04629

2.97801

16.07442 0.00288

1.680

2.000

0.08679

4.68413

23.47966 0.00412

1.720

2.000

0.12261

7.91592

27.23764 0.00450

1.740

2.000

0.14321 10.02282 33.99955 0.00421

1.760

2.000

0.18729 12.58832 42.44476 0.00441

1.780

2.000

0.22284 15.67655 49.43979 0.00451

1.800

2.000

0.44337 20.13860 63.82885 0.00695

limits as it is explained in section 8.3.3.

8.3.1.2

Case Study II: Gradually Increasing the Load Level and Keeping the
Generation Level Constant

This case study is similar to the case study I except that the generation level is kept constant
by a factor of 2. The reason of performing this case study is that in Case Study I, both
generation and load levels are increasing and the point at which the indices start sharply
increasing due to load effect is difficult to detect. The annualized reliability indices are shown
in Table 8.4.

119

0.45
0.4
0.35

Using AC Power Flow - Case I
Using DC Power Flow - Case I
Using AC Power Flow - Case II
Using DC Power Flow - Case II

System LOLP

0.3
0.25
0.2
0.15
0.1
0.05
0
0
10
20
30
40
50
60
70
80
Increase in the system loading as a percentage of the peak load (%)

Figure 8.1: System loss of load probability (LOLP) index profile against the change in the
system loading level
From Table 8.4, the system EDNS index smoothly increases and system loss of load
probability and loss of load frequency indices are smoothly increasing as the load level
increases until the load level reaches around 145% of the peak load. After this point, these
indices increase sharply. Also, a very sharp increase for these indices happens after around
175% of the peak load. Again, this sudden increase can be attributed to the voltage and
reactive power limits as it is explained in subsection 8.3.3.

8.3.2

Comparison with the DC Power Flow Model

To compare the results of the reliability indices with those found using the DC power flow
model, the analyses of the above case studies were repeated using the DC power flow model.
To better visualize the differences in these values, the results of the system loss of load
probability index are depicted in Figure 8.1.

120

Contributions of voltage limits on system LOLP

0.45
0.4
0.35

Effect of minimum voltage limits - case I
Effect of maximum voltage limits - case I
Effect of minimum voltage limits - case II
Effect of maximum voltage limits - case II

0.3
0.25
0.2
0.15
0.1
0.05
0
0
10
20
30
40
50
60
70
80
Increase in the system loading as a percentage of the peak load (%)

Figure 8.2: Contributions of the voltage limit violations on the system loss of load probability
(LOLP) index against the change in the system loading level
From Figure 8.1, it is obvious that the DC power flow model closely depicts the AC power
flow model, although optimistic, until the load reaches around 145% of the base peak load
level. After this point, the differences between the results of using the AC and DC power
flow models are significant where the results of the DC power flow are very optimistic.

8.3.3

Effects of Voltage and Reactive Power Limits on the Reliability Indices

From the above results, after a certain loading level, the two power flow models produce
significantly different results. Several factors may cause this difference. The effects of the
voltage and reactive power constraints on the loss of load probability index are investigated
through tracking the contributions of these constraints and the estimation of the proposed
indices of the voltage and reactive power limits.
121

Contributions of the reactive power limits in the LOLP

0.45
0.4
0.35

Effect of maximum reactive power limits - case I
Effect of minimum reactive power limits - case I
Effect of maximum reactive power limits - case II
Effect of minimum reactive power limits - case II

0.3
0.25
0.2
0.15
0.1
0.05
0
0
10
20
30
40
50
60
70
80
Increase in the system loading as a percentage of the peak load (%)

Figure 8.3: Contributions of the reactive power limit violations on the system loss of load
probability (LOLP) index against the change in the system loading level
8.3.3.1

Contributions of the Voltage and Reactive Power Constraints

The contributions of the voltage limit violations on the system loss of load probability index
against the change in the system loading level are shown in Figure 8.2. The contributions of
the reactive power limit violations on the system loss of load probability index against the
change in the system loading level are shown in Figure 8.3.
The results totally agree with the observations of the sharp increase in the loss of load
probability and loss of load frequency indices using AC power flow model as it is shown in
sections 8.3.1.1 and 8.3.1.2. The contributions of the voltage and reactive power violations
are very small until the load level reaches about 145% of the original load level for the Case
Study I and 170% for the Case Study II. Also, if the contributions of the voltage and reactive
power violations are subtracted from the loss of load probability index, both models would
produce similar results.
122

8.3.3.2

Indices of the Voltage and Reactive Power Limits

During the simulation, if the voltage or reactive power limits contributes in the load curtailment, the related index is updated according to (8.2). The results of these indices of the
Modified IEEE RTS are shown in Table 8.5. From Table 8.5, it is clear that the maximum
voltage limit at bus 2 has contributed in almost all the failure states. Maximum voltage
limits at buses 7, 13, 21 and 22 have contributed in around 80% of the failure states. One of
the possible solutions is to install condensers at these buses. Minimum voltage limit at bus 6
has contributed in around 83% of the failure states. One possible solution is to add a reactor
at this bus. Maximum reactive power limit at buses 2, 3 and 15 have contributed in around
90% of the failure states. The minimum reactive power limit at bus 6 has contributed in
around 15% of the failure states.

8.3.3.3

Relaxing Voltage Constraints

In this part, the voltage limit constraints are relaxed to investigate the effects of these
constraints on the reliability indices of the Modified IEEE RTS. Figure 8.4 shows the profile
of the loss of load probability index with relaxing the voltage constraints. The results of
the system indices with relaxing voltage constraints are given in Table 8.6. From Figure 8.4
and Table 8.6, it is clear that relaxing voltage limit constraints significantly reduces system
indices. However, even with relaxing these constraints, the values of these indices still high
in comparison with DC power flow counterpart. This can be related to the fact that the
DC power flow model ignores the losses in the lines and in the case of the AC power flow,
voltage constraints have been relaxed without increasing the reactive power support.

123

Table 8.5: Contributions of Voltage and Reactive Power Constraints
Bus No

8.3.4

υimax

υimin

imax

imin

1

0.15438 0.00000 0.00684 0.00004

2

0.99643 0.00000 0.84201 0.00000

3

0.00059 0.00941 0.96751 0.00000

4

0.00002 0.00016

–

–

5

0.00002 0.00007

–

–

6

0.00011 0.83354 0.00309 0.15427

7

0.85299 0.00000 0.03379 0.00000

8

0.00000 0.04804

–

–

9

0.00011 0.00787

–

–

10

0.04770 0.00000

–

–

11

0.00000 0.00034

–

–

12

0.00000 0.00598

–

–

13

0.78751 0.00000 0.12792 0.00000

14

0.00259 0.00435 0.03030 0.00002

15

0.18674 0.00000 0.91008 0.00015

16

0.07582 0.00050 0.19863 0.08940

17

0.00011 0.00000

18

0.20134 0.00000 0.08019 0.00000

19

0.00018 0.00000

–

–

20

0.00050 0.00000

–

–

21

0.82095 0.00000 0.32305 0.00000

22

0.29347 0.00007 0.00003 0.00002

23

0.82510 0.00402 0.02226 0.00007

24

0.00029 0.04100

–

–

–

–

Sensitivity Analysis

The results of the sensitivity analyses of the EDNS index with respect to the voltage and
reactive power constraints are shown in Table 8.7. The following points can be observed

124

0.5
0.45

System LOLP Index

0.4
0.35
0.3
0.25
0.2
0.15
0.1
1.05

1.1
Maximum Voltage Limits (p.u.)

1.15

Figure 8.4: The profile of system loss of load probability (LOLP) index with relaxing voltage
limit constraints
Table 8.6: Reliability Indices with Relaxing Voltage Limit Constraints
qˆ

ρˆ

φˆ

τˆ

MW

occ./yr

hour

Voltage

Voltage

Max.

Min.

1.05

0.95

0.44341 19.63547 62.98836 0.00704

1.06

0.94

0.28647 18.52728 58.25106 0.00492

1.07

0.93

0.23826 18.10026 45.03072 0.00440

1.08

0.92

0.19594 17.13641 44.27397 0.00443

1.09

0.91

0.18810 16.55377 42.44361 0.00443

1.10

0.90

0.16416 16.71187 36.94868 0.00444

1.11

0.89

0.15108 16.40434 33.24840 0.00454

1.12

0.88

0.14926 16.21602 33.35913 0.00447

1.13

0.87

0.14719 15.83632 32.79869 0.00449

1.14

0.86

0.14642 15.60547 32.05549 0.00457

1.15

0.85

0.14687 15.68066 32.06075 0.00458

125

Table 8.7: Sensitivity Analyses of the Severity Index with respect to the Voltage and Reactive
Power Constraints
Bus
No
1

Sensitivity of the Severity Index with respect to:
Vmax
Vmin
Qmax
Qmin
0.32604 0.00000 0.00003
0.00001

2

1.35209 0.00000 0.01412

0.00000

3

0.00002 0.17834 0.14772

0.00000

4

0.00000 0.00431

–

–

5

0.00000 0.00102

–

–

6

0.00005 10.3881 0.02302

0.16863

7

0.11922 0.00000 0.00088

0.00000

8

0.00000 0.10789

–

–

9

0.00000 0.02155

–

–

10

0.02656 0.00000

–

–

11

0.00000 0.00126

–

–

12

0.00000 0.03430

–

–

13

0.31398 0.00000 0.00579

0.00004

14

0.00745 0.00453 0.00120

0.00001

15

0.09631 0.00000 0.02306

0.00000

16

0.06129 0.00095 0.00106

0.00164

17

0.00000 0.00000

18

0.06080 0.00000 0.00071

19

0.00000 0.00000

–

–

20

0.00000 0.00000

–

–

21

0.81786 0.00000 0.00460

0.00000

22

0.14240 0.00000 0.00000

0.00000

23

0.51735 0.00125 0.00000

0.00005

24

0.00056 0.25576

–

–

126

–
0.00047

–

from the sensitivity analyses: (a) bus 2 has the highest effect in terms of maximum voltage
limit, therefore, the planner may consider installing a condenser at bus 2, (b) the minimum
voltage limit at bus 6 has the highest effect on the EDNS index, therefore, the planner may
consider installing a reactor, (c) bus 3 has the highest effect in terms of maximum reactive
power limit, and (d) bus 6 has the highest effect in terms of minimum reactive power limit.
It should be noted that the voltage sensitivity analyses are based on the p.u. values.

8.4

Summary

This chapter has investigated the effects of the voltage and reactive power constraints on
power system reliability through developing expressions and performing sensitivity analyses
of the EDNS index with respect to these constraints. Also, four indices were proposed to
address the contributions of the voltage and reactive power constraints on system load curtailments. Extensive studies on the IEEE RTS and the Modified IEEE RTS were conducted
to investigate the effects of the voltage and reactive power constraints on the power system reliability indices. Several case studies were performed through gradually stressing the
transmission lines to track the profile of the reliability indices with increasing system loading
level. A comparison between the use of the AC and DC power flow models in the composite
system reliability studies are provided to test the effects of ignoring the voltage and reactive
power constraints. The results show that, for the cases where the system is stressed, the
accuracy of the DC power flow in evaluating the reliability indices is deteriorated. Also, the
state space reduction technique has been used to reduce the computation time and burden.

127

Chapter 9
Hardening Power Systems Against
Cascading Failures Based on Risk
Sensitivity Analysis

9.1

Introduction

In response to the rapid growth of the global economics, large power grids are interconnected
in order to provide the required energy. As a result of these interconnections, the number of
failures also increases and may spread through these interconnections and cause cascading
failures. Large blackouts in power system networks are rare in nature, however they usually
start with a small triggering event that by successive sequence can cascade and develop a
complete blackout. Cascading failures in power systems have been recently encountered more
frequently. Examination of catastrophic events such as the August 14, 2003 blackout shows
that many factors were involved in the development of the cascading failure. These factors
can be classified into three categories: physical behavior of the networks, software failures and
operators’ mistakes. Among the physical behavior problems of the networks are overloading,
inadvertent tripping, transient instability, voltage collapse, etc. Such factors are interrelated
and their effect cannot be separated from each other. For instance, overloading may cause

128

inadvertent tripping and tripping of loaded lines may cause voltage collapse and/or transient
instability. An effort to address all of these factors faces many burdens such as modeling
and computational time. Most of current cascading failure mitigation methods can capture
only a subset of the cascading failure phenomena mechanisms. This chapter focuses on using
the sensitivity of the reliability risk index with respect to voltage limits and reactive power
reserve.

9.2

Sensitivity Analysis-based Cascading Failure Mitigation

This chapter introduces a method based on risk sensitivity analysis to determine strategies
to harden power systems against catastrophic failures. This method utilizes risk sensitivity
analysis in mitigating the cascading failures in power systems. In other words, the proposed
method utilizes sensitivity analysis of the risk index with respect to the control variables to
identify vulnerabilities that can be corrected by strategic or tactical measures. For any initiating event, the sensitivity of the risk of load curtailments with respect to the power quality
constraints such as voltage and reactive power is derived, and thereby preventive/correction
actions are determined. In addition to the insertion of compensation devices, rescheduling
of generation is considered to increase the reactive power reserve. In order to evaluate the
sensitivity of the risk index with respect to the power quality constraints, voltage and reactive power constraints are incorporated with the optimization problem with an objective of
minimum load curtailments.
Cascading failures usually occurs as a sequence of events within a time horizon of seconds to minutes and in the best scenarios the time horizon is hours. In order to evaluate
129

the sensitivity of the risk index with respect to the power quality within a short time, the
linearized power flow model that accounts for voltage and reactive power was used (chapter
2). The linearized power flow model was incorporated with a linear programming optimization problem with objective of minimum load curtailments. For any initiating event, the
sensitivity of the risk of load curtailment with respect to the control variables is used to
mitigate the risk of cascading failures. The risk sensitivity was combined with the WORMS
and the Dynamic Decision-Event Tress (DDET) techniques for searching for the chain of the
events and preventing the risk of the cascading failures.
From the sensitivity of the risk index with respect to voltage limits and reactive power
reserve, generation rescheduling can be performed to move the system from being close to
a voltage collapse to a state where the system resides within the safe operating limits. For
example, if one of the expected next events may cause the voltage at a certain bus will
be less than or larger than the operating limits, the remedial action can be reactive power
compensation at that bus. However, if there is no reactive power compensator nearby that
bus, the real power output of generators that are close to that bus have to be reduced to
increase the reactive power reserve. This chapter provides suggestions for the control action
but specifically identifying which corrective action should be performed. This model was
applied on the IEEE RTS and the results are presented.

9.3

Remedial Actions

A remedial action scheme based on the sensitivity analysis with economic dispatch helps
the operator in the control room to choose an action that appropriate for the conditions on
hand. The proposed sensitivity analysis-based remedial action scheme can be combined with

130

the other tools using multi-decision-making techniques to choose the “best” option that has
minimum risk and minimum cost. The remedial action scheme can be run periodically as
well as after every triggering event. From the results of the analyzer (sensitivity analysis of
the risk index with respect to voltage limits and reactive power reserve), a remedial action
scheme produces a list of actions that can alleviate the violation if there is any.
The remedial action scheme is intended to act as a preventive/corrective action control.
The ability of distinguishing between these two actions depends upon the initial conditions
of the system and the triggering events. For example, in case there is no disturbance,
the remedial action scheme examines the designated contingencies for possible voltage limit
violations. In case of possible violations, the scheme lists possible generation rescheduling
alternatives based on the sensitivity of the risk index with respect to control parameters
and the best action can be chosen based on multi-decision-making techniques. Possible
actions from the preventive action scheme are insertion of reactive power compensation and
generation rescheduling. The preventive action scheme can help operators to discard the
unsafe alternatives. Also, the suggested actions can be combined with the other tools and
control parameters such as equipment actuation time, cost, available means, etc. using
multi-decision-making techniques.
In case of corrective action mode, the scheme can be used to prevent the possibility of
cascading failures. The corrective action scheme is less costly than the preventive action;
however, it is the last defense line and if the correction actions suggested by this scheme
cannot be applied, the system more likely to go under cascading failures. Therefore, a
careful selection of the proposed actions needs to be implemented in this scheme. Possible
actions from the corrective action scheme are generation rescheduling and load shedding.
In many situations where loss of one of the critical lines may cause the system to go under
131

cascading failures and a fast action needs to be initiated. Based on the preventive action
scheme, the power flowing in the critical lines can be alleviated by generation rescheduling
but this might be costly prohibitive. In such cases, the critical lines will be allowed to work at
their power transfer limits and perform corrective action control in case of tripping of these
lines. The most commonly control actions is one or combination of generation tripping, load
shedding and braking resistor insertion.

9.4

Case Studies

The proposed formulation was applied on the modified IEEE RTS. System’s data are given
in the appendix.
It is well known that the transmission lines of the IEEE RTS are very reliable and they
have high power carrying capabilities in comparison with the loading level of the system.
The contributions of the failures of these transmission lines and their capacity limits on the
annual and annualized indices are very small and can be ignored [17, 99, 100]. For instance,
the annualized LOLP index of the IEEE RTS using DC load flow model is 0.08344. On the
other hand, if we relax the power carrying capability limits and assume the transmission lines
are perfectly reliable, the LOLP becomes 0.0829. The modified IEEE RTS is the same as the
IEEE RTS except that the generation is multiplied by 2.0 and the loads are multiplied by 1.8
to stress the transmission lines. However, these modifications cause the voltage limit at bus
3 to contribute in every sampled state. We have added a 30 MVar synchronous condenser
at bus 3.
In evaluating the sensitivity of the risk index with respect to the control variables, the
state space reduction technique was not able to entirely classify the state space into success

132

and failure subspaces. This can be interpreted as the unclassified subspace has failure states
with relatively high probability but with relatively small load curtailments. This also can
be attributed to the fact that failure states that are located far from the origin (all components are in the down state) have high probability and small load curtailments. Therefore,
the proposed dynamically directed binary particle swarm optimization search method was
used to search for the unclassified subspace. To investigate the benefits of the state space
reduction technique and the dynamically directed binary particle swarm optimization search
method, we have performed the same a case study on the modified IEEE RTS using Monte
Carlo simulation. The computation time using the state space reduction technique with the
dynamically directed binary particle swarm optimization search method was around 45%
less than the computation time using Monte Carlo simulation.
The sensitivities of the risk index of the modified IEEE RTS with respect voltage and
reactive power limits are shown in Table 9.1, Table 9.2, Table 9.3 and Table 9.4.
Table 9.1: Sensitivity of the risk index with respect to minimum reactive power limits for
the modified IEEE RTS
Power limit
Bus No.
Value

Minimum Reactive Power Limit
22

21

16

18

15

-0.00078 -0.00045 -0.00007 -0.00006 -0.00004

Table 9.2: Sensitivity of the risk index with respect to maximum reactive power limits for
the modified IEEE RTS
Power limit
Bus No.
Value

Maximum Reactive Power Limit
6

14

16

13 & 15

1

-0.00362 -0.00036 -0.00014 -0.00005 -0.00003

133

Table 9.3: Sensitivity of the risk index with respect to under voltage limits for the modified
IEEE RTS system
Voltage limit
Bus No.
Value
Bus No.
Value

Minimum Voltage Limit
3

8

6

24

7

-0.06538 -0.03226 -0.01649 -0.01345 -0.01313
4

18

20

19

13

-0.00795 -0.00218 -0.00032 -0.00027 -0.00032

Table 9.4: Sensitivity of the risk index with respect to over voltage limits for the modified
IEEE RTS system
Voltage limit
Bus No.
Value
Bus No.
Value

Maximum Voltage Limit
1

2

13

14

10

-0.03888 -0.02926 -0.02710 -0.01370 -0.01228
6

23

7

22

15

-0.01205 -0.01157 -0.00586 -0.00335 -0.00075

Form Table 9.1, we can see that among the minimum reactive power limits, bus 22
has the highest effect on the risk index. Also, buses 15, 16, 18, 21 and 22 are generation
buses. Therefore, one of the options to harden this system against catastrophic failures is to
reschedule the generation at these buses to increase the minimum reactive power limits.
Table 9.2 shows the effect of maximum reactive power limits on the risk index. Buses 1,
13, 15 and 16 are generator buses. Again, one of the options to harden this system against
catastrophic failures is to reschedule the generation at these buses to increase the maximum
reactive power limits. Bus 14 has a synchronous condenser. Another condenser may be
added at this bus. At bus 6 there is a fixed reactor with 100 MVar. A condenser may be
134

added at this bus.
Table 9.3 and Table 9.4 show the effect of under and over voltage limits on the risk index.
The buses that hit the limits consist of load and generation buses. At load buses that have
only under voltage problems, condensers may be added. At load buses that have only over
voltage problems, reactors may be added. As in the case of the effect of reactive power
limits, at generation buses that have under or over voltage problems or both, generation
rescheduling to increase the reactive power margin may solve the problem. Buses 6, 7 and
13 appear in both under and over voltage violations. These buses have special property.
Bus 6 is connected to bus 10 through a cable with high line susceptance. The high line
susceptance makes bus 6 sensitive to loading levels and line status. Buses 7 and 13 are
connected with a line with capacity of 175 MVA and it is the only line connects bus 7 to the
rest of the system. Any failure in this line, the system will be separated.

9.5

Summary

Several power system networks operate at their capacity and voltage limits and a tool to
harden power systems against cascading failures is of necessity. A risk sensitivity analysis
with respect to voltage and reactive power limits is introduced in this chapter. To account
for voltage and reactive power constraints, the linearized power flow model was used and
incorporated with a linear programming optimization problem with the objective of minimum
load curtailments to calculate the sensitivity of the reliability risk index with respect to
voltage limits and reactive power reserve. The most vulnerable components as well as the
most buses that have voltage collapse and reactive support problems were identified. The
results obtained from the risk sensitivity analysis are used to harden power systems against

135

cascading failures. In addition, these results can be used in preventive or corrective remedial
actions. It is noteworthy pointing out here that the method presented in this chapter differs
from the existing cascading failure mitigation methods in the sense that it evaluates the
sensitivity of the risk index with respect to voltage and reactive power limits for N-k events,
instead of considering sequence of events.
In performing the sensitivity analyses, large number of states has to be tested. The
state space reduction approach in conjunction with the dynamically directed particle swarm
optimization search method was used to reduce the search space, and thereby speed up the
computations. The proposed method was applied on the modified IEEE RTS. The results of
the sensitivity analysis of the risk index with respect to voltage and reactive power limits have
identified the buses that have most effect on system availability and possibility of cascading
failures.

136

Chapter 10
Conclusions and Future Work
Summaries of the findings have been provided at the end of each chapter. These summaries
discuss specific findings related to each chapter. This chapter provides a general conclusion
about the methods that have been used, the general outcomes of this work and suggestions.
Also, this chapter provides possible future developments that can be built on or added to
the presented work.

10.1

Concluding Remarks

The objective of this work was to develop a comprehensive formulation of the sensitivity
analysis and an efficient solution algorithms and methods to decrease the computation burden. The work presented in this thesis uses sensitivity analysis-based method to develop
strategies to improve power system reliability. Lagrange multipliers were used in estimating
the sensitivity of the reliability indices with respect to component parameters. The Lagrange
multipliers can be obtained from the dual solution of a standard linear/non-linear programming problem. Three power flow models were utilized in modeling power system networks,
namely, DC power flow model, AC power flow model and a linearized power flow model. A
state space reduction technique in conjunction with a dynamically directed particle swarm
optimization search method was used to reduce the computational time and burden.
The state space reduction technique, which is presented in chapter 4, was proposed and

137

applied to reduce the search space. This technique is conservative in that not all the failure
states can be captured and no success state was classified as failure state. The state space
reduction technique classifies the state space into three subspaces, definite failure subspace,
definite success subspace and unclassified subspace. An intelligent search method based on a
dynamically directed binary particle swarm optimization was developed and used to search
for failure states in the unclassified subspace. From the failure states, which were captured
by the state space reduction technique and the dynamically directed binary particle swarm
optimization search method, the reliability indices and their sensitivities with respect to
component parameters were calculated.
An element key in using binary particle swarm optimization in composite system reliability evaluation is the selection of weighting factors of the objective function. These weighting
factors are system dependent and, therefore, their appropriate values should be carefully
selected in order to prevent the swarm from being trapped to one corner of the state space.
An effective particle swarm optimization search method was proposed and used for reliability assessment of composite systems. The proposed method adjusts the weighting factors
associated with the objective function in a dynamic fashion so that the swarm would always
fly on the desired search space rather of being trapped to one corner of the search space.
The effectiveness of the proposed methods was tested on several standard test systems
and the results correspond closely with those obtained using Monte-Carlo simulation, while
requiring lower computational burden. The proposed methods were applied on several practical problems such as the effects of the voltage and reactive power limits on power system
reliability, the effects of imposing the regulations of reducing the emissions from power generation, and hardening power systems against cascading failures.
The effects of the voltage and reactive power constraints on power system reliability
138

were investigated through developing expressions and performing sensitivity analyses of the
expected demand not supplied (EDNS) index with respect to these constraints. Four indices
were proposed to address the contributions of the voltage and reactive power constraints on
system load curtailments. Extensive studies on the IEEE RTS and the Modified IEEE RTS
were conducted to investigate the effects of the voltage and reactive power constraints on
the power system reliability indices. Several case studies were performed through gradually
stressing the transmission lines to track the profile of the reliability indices with increasing
system loading level. A comparison between the use of the AC and DC power flow models
in the composite system reliability studies are provided to test the effects of ignoring the
voltage and reactive power constraints. The results show that, for the cases where the
system is stressed, the accuracy of the DC power flow in evaluating the reliability indices is
deteriorated.
The effects of the inclusion of the emission constraints were investigated through performing sensitivity analysis of the EDNS index with respect to the emission constraints.
Most of the emission studies in the literature have focused on minimizing emissions through
economic dispatch and introducing new concepts to solve multi-objective functions. This
work introduced a framework for evaluating the effect of the emission rates and allowances
on system reliability and determining the sensitivity of commonly used reliability indices to
equipment and emission constraints. From the sensitivity results, components can be ranked
from reliability point of view. Therefore, components that have the highest effect on system
reliability can be identified. Also, if the reliability of the system needs to be improved, the
sensitivity analyses can identify which component to increase its emission limit. Inclusion
of emission constraints in planning studies helps guide investment strategies with respect to
reconstruction or upgrades.
139

Hardening power systems against cascading failures by improving voltage stability was
also introduced. The sensitivity of the reliability EDNS index with respect to voltage limits
and reactive power reserve was used to identify the most vulnerable components as well as
the most buses that have voltage collapse and reactive support problems. In addition, the
results of the sensitivity analysis can be used in establishing preventive or corrective remedial
action schemes.

10.2

Future Work

Although the results presented in this thesis show the powerfulness of the use of the sensitivity
analysis in power system planning and operation and the effectiveness of the use of state
space reduction technique with the population-based intelligent search methods, the work
could be extended and further developed by considering some of the following points:
1. Adding another constraints: Another constraints can be added to the system operating
constraint such as transient stability limits, voltage stability limits, operating cost, etc.
Inclusion of such constraints would increase the complexity of the optimization problem
but it would make the solution for optima more realistic.
2. Using different population-based intelligent search methods: In this thesis, we have
used particle swarm optimization due to its simplicity; the other population-based
intelligent search methods could be used and a comparison between these methods in
terms of efficiency and reduction in the computational time could be established.
3. Combining the proposed state space reduction method with the existing methods: In
this work, we have used minimum generation vector for determining the boundary of

140

the success subspace. The state space pruning could be combined with the proposed
method to decrease the size of the unclassified subspace.
4. Determining the effective region of the Lagrange multipliers: In this work, we have used
the value of the Lagrange multipliers in performing sensitivity analysis, but without
accounting for the effective region. The effective region can be defined as the region in
which the Lagrange multipliers of the corresponding constraints are binding. In other
words, there will be a limit for a Lagrange multiplier after which the corresponding
constraint becomes redundant.

141

APPENDIX

142

Derivation of the Sensitivity Relationships
Here, only the derivations of the sensitivity of the LOLP index with respect to reliability
parameters are provided. The remaining sensitivities for LOLF and EDNS can be easily
derived in the same manner. Some of the derivations presented below have been adapted
from [32]. However, this reference does not provide the complete derivations.

Sensitivity of LOLP to Component Failure Probability

∂LOLP/∂ui
{[∂If (x)/∂ui ]P (x) + If (x)[∂P (x)/∂ui ]}

=
x∈X

If (x){[∂[P (S1 )P (S1 ) · · · P (Si ) · · · P (Sm )]/∂ui ]}

=
x∈X

=

If (x)P (x)[1/P (Si )][∂P (Si )/∂ui ]
x∈X

If Si = 1,

[1/P (Si )][∂P (Si )/∂ui ]
= (1/P (Si = 1))(∂P (Si = 1)/∂ui )
= (1/ai )(∂ai /∂ui ) = −1/ai = −Si /ai

143

If Si = 0,

[1/P (Si )][∂P (Si )/∂ui ]
= (1/P (Si = 0))(∂P (Si = 0)/∂ui )
= (1/ui )(∂ui /∂ui ) = 1/ui = (1 − Si )/ui

Therefore, in general,

[1/P (Si )][∂P (Si )/∂ui ] = (1 − Si )/ui − Si /ai
= (1/ui ) − Si /(ai ui )

Then,
If (x)P (x)[(1/ui ) − Si /(ai ui )]

∂LOLP/∂ui =
x∈X

Sensitivity of LOLP to Component Failure Rate

∂LOLP/∂λi
{[∂If (x)/∂λi ]P (x) + If (x)[∂P (x)/∂λi ]}

=
x∈X

=

If (x){[∂[P (S1 )P (S1 )...P (Si )...P (Sm )]/∂λi ]}
x∈X

=

If (x)P (x)[1/P (Si )][∂P (Si )/∂λi ]
x∈X

144

If Si = 1,

[1/P (Si )][∂P (Si )/∂λi ]
= (1/P (Si = 1))(∂P (Si = 1)/∂λi )
= (1/ai )(∂ai /∂λi ) = −Si /(λi + µi )

If Si = 0,

[1/P (Si )][∂P (Si )/∂λi ]
= (1/P (Si = 0))(∂P (Si = 0)/∂λi )
= (1/ui )(∂ui /∂λi )
= (ai /ui )(1/(λi + µi ))
= (ai /ui )((1 − Si )/(λi + µi ))

Therefore, in general,

[1/P (Si )][∂P (Si )/∂λi ]
= −Si /(λi + µi ) + (ai /ui )((1 − Si )/(λi + µi ))
= ai /λi − Si /λi

Then,
If (x)P (x)[ai /λi − Si /λi ]

∂LOLP/∂λi =
x∈X

145

Sensitivity of LOLP to Component Repair Rate

∂LOLP/∂µi
{[∂If (x)/∂µi ]P (x) + If (x)[∂P (x)/∂µi ]}

=
x∈X

=

If (x){[∂[P (S1 )P (S1 )...P (Si )...P (Sm )]/∂µi ]}
x∈X

=

If (x)P (x)[1/P (Si )][∂P (Si )/∂µi ]
x∈X

If Si = 1,

[1/P (Si )][∂P (Si )/∂µi ] = (1/P (Si = 1))(∂P (Si = 1)/∂µi )
= (1/ai )(∂ai /∂µi )
= (ui /ai )(1/(λi + µi ))
= (ui /ai )(Si /(λi + µi ))

If Si = 0,

[1/P (Si )][∂P (Si )/∂µi ]
= (1/P (Si = 0))(∂P (Si = 0)/∂µi )
= (1/ui )(∂ui /∂µi )
= −1/(λi + µi ) = −(1 − Si )/(λi + µi )

146

Therefore, in general,

[1/P (Si )][∂P (Si )/∂µi ]
= (ui /ai )(Si /(λi + µi )) − (1 − Si )/(λi + µi )
= −ai /µi + Si /µi

Then,
∂LOLP/∂µi =

If (x)P (x)[−ai /µi + Si /µi ]
x∈X

Test Systems
IEEE RTS
The IEEE RTS System has been extensively tested for power system reliability analysis [63].
IEEE RTS consists of 24 buses, 38 transmission lines/transformers (33 transmission lines and
5 transformers) and 32 generating units on 10 buses. The total generation of this system is
3405 MW and total peak load is 2850 MW.

Modified IEEE
The modified IEEE RTS System is the same as the original IEEE RTS System except that
the generation is doubled and the loads are multiplied by a factor of 1.8. The reason of
this modification is because the transmission lines of the original system have high power
carry capabilities in comparison with the generation and loading conditions. Therefore, this
modification will make the transmission lines more stressed. Also, this modification will
allow us to test the robustness of the proposed method.
147

A Six Bus Hypothetical System
The six bus system is a hypothetical system consists of 16 generating units on two generating
buses, five load buses and 9 transmission lines. Line and generation data are given in [64].
The total generation of this system is 240 MW and the total real power load is 185 MW.

The Saskatchewan Power Corporation in Canada
The SPC system is the network of the Saskatchewan Power Corporation in Canada [64]. The
system consists of 45 buses, 29 generating units on 8 buses and 71 transmission lines/transformers.
Four of the 45 buses are used to represent equivalent assistance from the Manitoba Hydro
System. One of these four buses is a fictitious bus that represents the power import from
Manitoba Hydro System which is 300 MW. This bus is connected to the other three buses
and the sum of the 300 MW is represented by three generating units each of which 100 MW.
The total generation is 2530 MW and total load is 1802.5 MW. System data are given in [64].

148

BIBLIOGRAPHY

149

BIBLIOGRAPHY
[1] S. Elsaiah, M. Benidris and J. Mitra, “Analytical approach for placement and sizing
of distributed generation on distribution systems,” IET Generation, Transmission &
Distribution, vol. 8, no. 6, pp. 1039–1049, Dec. 2013.
[2] T. N. Santos and A. L. Diniz, “A dynamic piecewise linear model for DC transmission
losses in optimal scheduling problems,” IEEE Trans. Power Syst., vol. 26, no. 2, pp.
508–519, May 2011.
[3] H. Zhang, V. Vittal, G. T. Heydt and J. Quintero, “A mixed-integer linear programming approach for multi-stage security-constrained transmission expansion planning,”
IEEE Trans. Power Syst., vol. 27, no. 2, pp. 1125–1133, May 2012.
[4] N. Alguacil, A. L. Motto and A. J. Conejo, “Transmission expansion planning: a
mixed-integer LP approach,” IEEE Trans. Power Syst., vol. 18, no. 3, pp. 1070–1077,
Aug. 2003.
[5] Wenyuan Li, Risk Assessment of Power Systems. USA, IEEE Press:
Interscience, 2005.

Wiley-

[6] C. Singh and J. Mitra, “Composite System Reliability Evaluation Using State Space
Pruning,” IEEE Trans. on Power System, vol. 12, no. 1, pp. 471–479, Feb. 1997.
[7] J. Mitra and C. Singh, “Incorporating the DC Load Flow Model in the Decomposition–
Simulation Method of Multi-Area Reliability Evaluation,” IEEE Trans. on Power
Syst., vol. 11, no. 3, pp. 1245–1254, Aug. 1996.
[8] J. Mitra and C. Singh, “Pruning and Simulation for Determination of Frequency and
Duration Indices of Composite Systems,” IEEE Trans. on Power Syst., vol. 14, no. 3,
pp. 899–905, Aug. 1999.
[9] J. He, Y. Sun, D. S. Kirschen, C. Singh and L. Cheng, “State-space partitioning method
for composite power system reliability assessment,” IET Gener. Transm. Distrib., vol.
4, no. 7, pp. 780–792, 2010.

150

[10] M. Benidris and J. Mitra, “Composite Power System Reliability Assessment Using
Maximum Capacity Flow and Directed Binary Particle Swarm Optimization,” Proc.
North American Power Symposium, pp. 1–6, 2013.
[11] M. Benidris, S. Elsaiah and J. Mitra, “Composite System Reliability Assessment Using
Dynamically Directed Particle Swarm Optimization,” Proc. North American Power
Symposium, pp. 1–6, 2013.
[12] S. B. Patra, J. Mitra and R. Earla, “A New Intelligent Search Method for Composite
System Reliability Analysis,” Proc. of the IEEE-PES Transmission and Distribution
Conference and Exposition, Dallas, TX, May 21–24, 2006, pp 803–807.
[13] J. Mitra and X. Xu, “Composite System Reliability Analysis Using Particle Swarm
Optimization,” IEEE 11th International Conference on Probabilistic Methods Applied
to Power Systems (PMAPS), pp. 548–552, 2010.
[14] L. Wang and C. Singh, “Population-Based Intelligent Search in Reliability Evaluation
of Generation Systems With Wind Power Penetration,” IEEE Trans. on Pow. Syst.,
vol. 23, no. 3, pp. 1336–1345, 2008.
[15] J. Mitra and M. R. Vallem, “Determination of storage required to meet reliability
guarantees on island-capable Microgrids with intermittent sources,” IEEE Trans. Pow.
Syst., vol. 27, no. 4, pp. 2360–2367, Nov. 2012.
[16] M. Benidris and J. Mitra, “Reliability and Sensitivity Analysis of Composite Power
Systems Under Emission Constraints,” IEEE Trans. Power Syst., vol. 29, no.1, pp.
404–412, Jan. 2014.
[17] M. Benidris, S. Elsaiah and J. Mitra, “Sensitivity Analysis in Composite System Reliability Using Weighted Shadow Prices,” IEEE PES General Meeting 2011, pp. 1–6,
Detroit, MI, USA.
[18] A. C. G. Melo, M. V. F. Pereira and A. M. Leite da Silva, “A Conditional Probability
Approach to the Calculation of Frequency and Duration Indices in Composite System
Reliability Evaluation,” IEEE Trans. Power Syst., vol. 8, no. 3, pp. 1118–1125, Aug.
1993.
[19] A. C. G. Melo, M. V. F. Pereira and A. M. Leite da Silva, “Frequency and Duration
Calculation in Composite Generation and Transmission Reliability Evaluation,” IEEE
Trans. Power Syst., vol. 7, no. 2, pp. 469–476, May 1992.

151

[20] C. Singh, “Calculating the Time-Specific Frequency of System Failure,” IEEE Trans.
on Reliability, vol. R-28, no. 2, pp. 124–126, June 1979.
[21] C. Singh, “Rules for Calculating the Time-Specific Frequency of System Failure,” IEEE
Trans. on Reliability, vol. 30, no. 4, pp. 364–366, Oct. 1981.
[22] R. Billinton and W. Li, Reliability Assessment of Electric Power Systems Using Monte
Carlo Methods, Plenum Press, New York, USA 1994.
[23] C. Singh and Q. Chen, “Generation System Reliability Evaluation Using a Cluster
Based Load Model,” IEEE Trans. on Power Systems, Vol. 4, No. 1, pp. 102–107, Feb.
1989.
[24] C. Singh and A. Lago-Gonzalez, “Improved Algorithms for Multi-Area Reliability Evaluation Using the Decomposition-Simulation Approach,” IEEE Trans. on Power Systems, vol. 4, no. 1, pp. 321–328, Feb. 1989.
[25] A. Lago-Gonzalez and C. Singh, “The Extended Decomposition–Simulation Approach
for Multi-Area Reliability Calculations,” IEEE Trans. on Power Systems, vol. 5, no.
3, pp. 1024–1031, Aug. 1990.
[26] C. Singh and Z. Deng, “A New Algorithm for Multi-Area Reliability Evaluation?Simultaneous Decomposition–Simulation Approach,” Electric Power Systems Research, vol. 41, no. 2, pp. 129–136, June 1991.
[27] P. Ramanujan, Z. Li, and L. Higham, “Shadow prices versus Vickery prices in multipath
routing,” IEEE conference on Information and Communication, pp. 2956–2960, 2009.
[28] S. Cunda, M. Pereria, L. Pinto and G. Oliveira, “Composite Generation and Transmission Reliability Evaluation in Large Hydroelectric Systems,” IEEE Trans. Power
Appar. and Syst., vol. 104, no. 10, pp. 2657–2663, Oct. 1985.
[29] G.Li, C. Liu, and H. Salazar, “Forecasting transmission congestion using day-ahead
shadow prices,” IEEE Power Syatems Conference, pp. 1705–1709, 2006.
[30] Z. Hui, Y. Pei, and L. Jie, “Study on the estimating model of shadow prices for
national economy evaluation of construction projects,” IEEE, International Conference
on Management and Service Science, pp. 1–4, 2009.
[31] A. C. G. Melo and M. V. F. Pereira, “Sensitivity Analysis of Reliability Indices with
Respect to Equipment Failure and Repair Rates,” IEEE Trans. Power Syst., vol. 10,
no.2, pp. 1014–1021, May 1995.
152

[32] Y. Zhao, N. Zhou, J. Zhou and X. Zhao, “Research on Sensitivity Analysis for Composite Generation and Transmission System Reliability Evaluation,” IEEE International
Conference on Power System Technology, pp. 1–5, 2006.
[33] T. X. Zhu, “A New Methodology of Analytical Formula Deduction and Sensitivity
Analysis of EENS in Bulk Power System Reliability Assessment,” IEEE PES, Power
Systems Conference and Exposition, PSCE, pp. 825–831, 2006.
[34] M. V. F. Pereira and L. M. V. G. Pinto, “Application of Sensitivity Analysis of Load
Supplying Capability to Interactive Transmission Expansion Planning,” IEEE Trans.
Power Appar. and Syst., vol. PAS–104, no.2, pp. 381–389, Feb. 1985.
[35] Z. Deng and C. Singh, “A new approach to reliability evaluation of interconnected
power systems including planned outages and frequency calculations,” IEEE Trans.
on Power Systems, vol. 7, no. 2, pp. 734–743, May 1992.
[36] J. Mitra, M. R. Vallem and S. B. Patra, “A Probabilistic Search Method for Optimal Resource Deployment in a Microgrid,” Proc. of the 9th International Conference
on Probabilistic Methods Applied to Power Systems, Stockholm, Sweden, June 11–15,
2006.
[37] G. C. Oliveira, M. V. F. Pereira, and S. H. F. Cunha, “A technique for reducing
computational effort in Monte Carlo based composite reliability evaluation,” IEEE
Trans. on Power Systems, vol. 4, no. 4, pp. 1309–1315, Nov. 1989.
[38] R. C. Green, Z. Wang, L. Wang, M. Alam, and C. Singh, “Evaluation of Loss of
Load Probability for Power Systems Using Intelligent Search Based State Space Pruning,” IEEE 11th International Conference on Probabilistic Methods Applied to Power
Systems (PMAPS), pp. 319–324, 2010.
[39] R. C. Green, L. Wang and C. Singh, “State Space Pruning for Power System Reliability Evaluation Using Genetic Algorithms,” IEEE Power and Energy Society General
Meeting, pp. 1–6, 2010.
[40] D. Zhao and C. Singh, “Modified Genetic Algorithm in State Space Pruning for Power
System Reliability Evaluation and Its Parameter Determination,” IEEE North American Power Symposium, pp. 1–6, 2010.
[41] R. C. Green, L. Wang, Z. Wang, M. Alam, and C. Singh, “Power System Reliability
Assessment Using Intelligent State Space Pruning Techniques: A Comparative Study,”
Proc. IEEE International Conference on Power System Technology (POWERCON),
pp. 1–8, 2010.
153

[42] R. C. Green, L. Wang, M. Alam and C. Singh, “State Space Pruning for Reliability
Evaluation Using Binary Particle Swarm Optimization,” IEEE/PES Power Systems
Conference and Exposition, pp. 1–7, 2011.
[43] R. C. Green, L. Wang, M. Alam and C. Singh, “An Examination of Artificial Immune
System Optimization in Intelligent State Space Pruning for LOLP Estimation,” IEEE
North American Power Symposium (NAPS), pp. 1–7, 2011.
[44] R. C. Green, L. Wang, M. Alam and C. Singh, “Intelligent State Space Pruning Using
Multi-Objective PSO for Reliability Analysis of Composite Power Systems: Observations, Analyses, and Impacts,” IEEE Power and Energy Society General Meeting, pp.
1–8, 2011.
[45] J. Mitra, “Application of Computational Intelligence in Optimal Expansion of Distribution Systems,” Proc. of the IEEE-PES Annual General Meeting, Pittsburgh, PA,
July 20–24, 2008.
[46] J. Mitra, S. B. Patra, M. R. Vallem and S. J. Ranade “Optimization of Generation
and Distribution Expansion in Microgrid Architectures,” Proc. of the 6th WSEAS
International Conference on Power Systems, Lisbon, Portugal, Sept. 22–24, 2006, pp.
14–21.
[47] S. A. Al-Askari, S. J. Ranade and J. Mitra, “Optimal Allocation of Shunt Capacitors
Placed in a Microgrid Operating in the Islanded Mode,” Proc. of the 37th annual North
American Power Symposium, Ames, IA, Oct. 23–25, 2005, pp. 406–411.
[48] X. Xu and J. Mitra, “Reliability Evaluation of Composite System Using Bacterial Foraging Algorithm,” Proc. of the 17th Power System Computation Conference, Stockholm, Sweden, August 22–26, 2011.
[49] S. J. Ranade and R. Kolluru and J. Mitra, “Identification of Chains of Events Leading
to Catastrophic Failures of Power Systems,” Proc. of the IEEE International Symposium on Circuits and Systems, 2005, Kobe, Japan, May 23–26, 2005, pp. 4187–4190.
[50] Y. del Valle, G. K. Venayagamoorthy, S. Mohagheghi, J.-C. Hernandez and R. G.
Harley, “Particle swarm optimization: basic concepts, variants and applications in
power systems,” IEEE Trans. Evol. Comp., vol. 12, no. 2, pp. 171–195, Apr. 2008.
[51] M. Benidris and J. Mitra, “Use of Intelligent Search Methods in Performing Sensitivity
Analysis of Power System Reliability Indices,” Proc. of the IEEE-PES Annual General
Meeting, Washington, DC, July 27–31, 2014.

154

[52] R. Earla, S. B. Patra and J. Mitra, “A Particle Swarm Based Method for Composite
System Reliability Analysis,” Proc. of the 36th annual North American Power Symposium, Moscow, ID, Aug. 8–10, 2004, pp 294–298.
[53] V. Miranda, L. M. Carvalho, M. A. Rosa, A. M. L. Silva and C. Singh, “Improving
Power System Reliability Calculation Efficiency With EPSO Variants,” IEEE Trans.
on Pow. Syst., vol. 24, no. 4, pp. 1772–1779, 2009.
[54] F. Yang and C. S. Chang, “Optimisation of maintenance schedules and extents for
composite power systems using multi-objective evolutionary algorithm,” IET Gener.
Transm. Distrib., vol. 3, no. 10, pp. 930–940, 2009.
[55] K. Suresh and N. Kumarappan, “Generation maintenance scheduling using improved
binary particle swarm optimisation considering aging failures,” IET Gener. Transm.
Distrib., vol. 7, no. 10, pp. 1072–1086, 2013.
[56] S. Elsaiah, M. Benidris and J. Mitra, Reliability-Constrained Optimal Distribution System Reconfiguration, to appear in Computational Intelligence Applications in Modeling
and Control, Springer, 2015.
[57] J. Mitra, S. B. Patra, S. J. Ranade and M. R. Vallem, “Reliability-Specified Generation and Distribution Expansion in Microgrid Architectures,”WSEAS Transactions on
Power Systems, vol. 1, no. 8, pp 1446–1453, Aug. 2006.
[58] M. R. Vallem and J. Mitra, “A Particle Swarm-based Method for Reliability-driven
DG Deployment and Feeder Augmentation in Microgrids,” Proc. of the 15th National
Power System Conference, Mumbai, India, Dec. 16–18, 2008.
[59] J. Mitra, S. B. Patra and S. J. Ranade, “Reliability Stipulated Microgrid Architecture Using Particle Swarm Optimization,” Proc. of the 9th International Conference
on Probabilistic Methods Applied to Power Systems, Stockholm, Sweden, June 11–15,
2006.
[60] J. Kennedy and R. Eberhart, “Particle Swarm Optimization,” IEEE International
Conference on Neural Networks, pp. 1942–1948, 1995.
[61] J. Kennedy and R. Eberhart, “A Discrete Binary Version of the Particle Swarm Algorithm,” IEEE International Conference on Computational Cybernetics and Simulation,
pp. 4104–4108, 1997.

155

[62] B. Li, J. Qu and Y. Cui, “The resources purchase quantity algorithm based on the
shadow price and its application in the enterprise management,” IEEE International
Conference on Test and Management, pp. 346–349, 2009.
[63] Reliability Test System Task Force of the Application of Probability Methods Subcommittee, “IEEE Reliability Test System,” IEEE Trans. on Power Apparatus and
Systems, Vol. PAS–98, No. 6, pp. 2047–2054, Nov./Dec. 1979.
[64] R. Billinton and S. Kumar, “Effect of Higher–Level Independent Generator Outages
in Composite-System Adequacy Evaluation,” IEE Proceedings on Generation, Transmission and Distribution, vol. 134, no. 1, pp. 17–26, Jan. 1987.
[65] F. Lee, J. Liao and A. Breipohl, “Coordination of SO2 emission allowance trading,
energy and spinning reserve transactions, and consumption of take–or–pay fuels,” IEEE
Trans. on Power Systems, vol. 9, no. 3, pp. 1243–1252, Aug. 1994.
[66] Center for Climate and Energy Solutions. (2013,
January). Califorst
nia Cap–And–Trade Program Summary (1
ed.) [Online]. Available:
http://www.c2es.org/docUploads/calif-cap-trade-01-13.pdf
[67] Regional Greenhouse Gas Initiative, Inc. (2013, February 7). RGGI
States Propose Lowering Regional CO2 Emissions Cap 45%, Implementing a More Flexible Cost–Control Mechanism (1st ed.) [Online] Available:
http://www.rggi.org/docs/PressReleases/PR130207 ModelRule.pdf
[68] L. M. Brown, A. Hanafi and A. Petsonk. (2012). The EU Emissions
Trading System–Results and Lessons Learned (1st ed.) [Online]. Available:
http://www.edf.org/climate/eu-emissions-trading-system-report
[69] J. Kabouris and F. Kanellos, “Impacts of Large–Scale Wind Penetration on Designing
and Operation of Electric Power Systems,” IEEE Trans. on Sustainable Energy, vol.
1, no. 2, pp. 107–114, July 2010.
[70] M. Gent and J. Lamont, “Minimum–Emission Dispatch,” IEEE Trans. on Power Apparatus and Systems, vol. PAS–90, no. 6, pp. 2650–2660, June 1971.
[71] J. Nanda, L. Hari and M. Kothari, “Economic emission load dispatch with line flow
constraints using a classical technique,” IEE Proceedings—Generation, Transmission
and Distribution, vol. 141, no. 1, pp. 1–10, Jan. 1994.

156

[72] M. AlRashidi and M. El-Hawary, “Emission–Economic Dispatch Using a Novel Constraint Handling Particle Swarm Optimization Strategy,” IEEE Canadian Conference
on Electrical and Computer Engineering, CCECE ’06., pp. 664–669, 2006.
[73] H. Rughooputh and R. Ah King, “Environmental/economic dispatch of thermal units
using an elitist multiobjective evolutionary algorithm,” IEEE International Conference
on Industrial Technology, vol. 1, pp. 48–53, 2003.
[74] L. Wang and C. Singh, “Reserve–Constrained Multiarea Environmental/Economic Dispatch Using Enhanced Particle Swarm Optimization,” IEEE, Systems and Information
Engineering Design Symposium, pp. 96–100, 2006.
[75] T. Prasanna and P. Somasundaram, “Fuzzy mutated evolutionary programming based
algorithm for combined economic and emission dispatch,” IEEE Region 10 Conference,
pp. 1–5, 2008.
[76] M. Abido, “Environmental/Economic Power Dispatch Using Multiobjective Evolutionary Algorithms,” IEEE Trans. on Power Systems, vol. 18, no. 4, pp. 1529–1537, Nov.
2003.
[77] J. Choi, S. Jeong, S. Bo, J. Park, A. El-Keib and J. Watada, “CO2, SOx and NOx
Emissions Constrained Multi–Criteria–Best Generation Mix Using Fuzzy Set Theory,”
IEEE Conference on Large Engineering Systems, pp. 118–124, 2007.
[78] D. Yamashita, T. Niimura, R. Yokoyama and M. Marmiroli, “Trade–off Analysis of
CO2 versus Cost by Multi–objective Unit Commitment,” IEEE PES General Meeting,
pp. 1–6, 2010.
[79] W. Spens and F. Lee, “Interactive search approach to emission constrained dispatch,”
IEEE Trans. on Power Systems, vol. 12, no. 2, pp. 811–817, May 1997.
[80] R. Sun, L. Chen, Y. Song, and Y. Sun, “Generation Expansion Planning based on
multi–area reliability exponential analytic model and emission control,” IEEE/PES
Transmission and Distribution Conference and Exposition, pp. 1–6, 2008.
[81] A. Prasai, A. Paquette, Yi Du, R. Harley and D. Divan, “Minimizing emissions in
microgrids while meeting reliability and power quality objectives,” IEEE International
Power Electronics Conference (IPEC), pp. 726–733, 2010.
[82] P. Venkatesh, R. Gnanadass and N. P. Padhy, “Comparison and Application of Evolutionary Programming Techniques to Combined Economic Emission Dispatch With
157

Line Flow Constraints,” IEEE Trans. on Power Systems, vol. 18, no. 2, pp. 688–697,
May 2003.
[83] J. R. Domike and A. C. Zacaroli, The Clean Air Act Handbook, American Bar Association, Chicago, USA, 2011.
[84] W. Graus and E. Worrell, “Trend in efficiency and capacity of fossil power generation
in the EU,” Elsevier, Energy Policy, vol. 37, no. 6, pp. 2147–2160, June 2009. Bar
Association, Chicago
[85] X. Yang, Y. Song, G. Wang and W. Wang, “A Comprehensive Review on the Development of Sustainable Energy Strategy and Implementation in China,” IEEE Trans.
on Sustainable Energy, vol. 1, no. 2, pp. 57–65, July 2010.
[86] Y. Kim and B. Ahn, “Multicriteria generation–expansion planning with global environmental considerations,” IEEE Trans. on Engineering Management, vol. 40, no. 2,
pp. 154–161, Aug. 2002.
[87] V. Mendes, S. Mariano, J. Catalao and L. Ferreira, “Emission constraints on short–
term schedule of thermal units,” Universities Power Engineering Conference, 39th
International, vol. 2, pp. 1068–1072, 2004.
[88] Z. Yong, “Control of CO2 Emissions of China Under Kyoto Protocol,” IEEE International Conference on Energy and Environment Technology, ICEET, pp. 32–35, 2009.
[89] A. Farag, S. Al-Baiyat, and T. C. Cheng, “Economic load dispatch multiobjective
optimization procedures using linear programming techniques,” IEEE Trans. Power
Syst., vol. 10, pp. 731–738, May 1995.
[90] D. B. Das and C. Patvardhan, “New multi-objective stochastic search technique for
economic load dispatch,” Proc. Inst. Elect. Eng.-Gen. Transm. Dist., vol. 145, no. 6,
pp. 747–752, 1998.
[91] H. Chang, C. Kuo, J. Jiang and S. Lu, “Generation Dispatch Strategy under Environment Protection Consideration,” IEEE/ACM International Conference on Green
Computing and Communications, pp. 242–247, 2011.
[92] P. L. Noferi, L. Paris “Effects of Voltage and Reactive Power Constraints on Power
System Reliability,” IEEE Trans. Power Appar. and Syst., vol. PAS-94, no. 2, pp.
482–490, March/April 1975.

158

[93] W. Qin, P. Wang, X. Han and X. Du, “Reactive Power Aspects in Reliability Assessment of Power Systems,” IEEE Trans. Power Syst., vol. 26, no. 1, pp. 85–92, Feb.
2011.
[94] P. Wang, W. Qin, X. Han and X. Du, “Reliability Assessment of Power Systems
Considering Reactive Power Sources,” IEEE Power & Energy Society General Meeting,
pp. 1–7, 2009.
[95] W. Qin, P. Wang, J. Song and Z. Wang, “Reactive Power Impact on Reliability of
220kV Taiyuan Power System,” IEEE 11th International Conference on Probabilistic
Methods Applied to Power Systems (PMAPS), pp. 648–653, 2010.
[96] D. Gaikwad and S. Mehraeen, “Reactive Power Considerations in Reliability Analysis
of Photovoltaic Systems,” IEEE Green Technologies Conference, pp. 1–6, 2012.
[97] C. Singh, T. Pravin, M. P. Bhavaraju and M. Lauby, “A Monte Carlo Tool for Estimating Contingency Statistics,” In Proceedings of the Joint International Power
Conference, vol. 2, pp. 670–674, Athens Power Tech, 1993.
[98] R. Billinton and W. Li, “A System State Transition Sampling Method for Composite
System Reliability Evaluation,” IEEE Trans. Power Syst., vol. 8, no. 3, pp. 761–770,
Aug. 1993.
[99] M. Benidris and J. Mitra, “Composite Power System Reliability Assessment Using
Maximum Capacity Flow and Directed Binary Particle Swarm Optimization,” North
American Power Symposium, pp. 1–6, 2013.
[100] M. Benidris, S. Elsaiah and J. Mitra, “Composite System Reliability Assessment Using
Dynamically Directed Particle Swarm Optimization,” North American Power Symposium, pp. 1–6, 2013.

159