IMPROVING THE CAPACITY VALUE OF WIND POWER WITH ENERGY
STORAGE INTEGRATION
By
Elicia A. Sashington

A THESIS
Submitted to
Michigan State University
In partial fulfillment of the requirements
For the degree of
Electrical Engineering-Master of Science
2014

ABSTRACT
IMPROVING THE CAPACITY VALUE OF WIND POWER WITH ENERGY
STORAGE INTEGRATION
By
Elicia A. Sashington
Wind power generation is a steadily growing renewable energy resource.
However, the integration of wind power to the power grid poses reliability issues because wind is
variable and unpredictable. One of the reliability issues is the mismatch between the system’s
load and generation. A solution to alleviating this issue is to combine wind power with energy
storage systems. This combination will reduce the variability of wind power generation output
and improve the capacity value. Capacity value can be defined as the amount of additional load
that can be served due to the addition of a generator, while maintaining the existing levels of
reliability [1]. For planning and reliability purposes, the capacity contribution of wind power
with energy storage systems must be estimated. This will transform wind power from an energy
resource to a capacity resource.
In this thesis, an analysis of systems with wind power and wind power with energy
storage were performed using reliability methods. In addition, a descriptive analysis and
comparison was completed on different types of energy storage systems. The objectives of this
thesis were to show the improvement of capacity values when energy storage systems are
combined with wind power and energy storage systems that are most suitable for the system
being analyzed.

Copyright by
ELICIA A. SASHINGTON
2014

To my heavenly father, who has blessed me with the knowledge and strength to achieve
this higher goal in my life, thank you for the gift of learning and tapping into my potential. To
my mother, Sandra Burrage; I would like to thank you for supporting me and being the best
parent you could be. To my Aggie friends, Julian Alford and Lauren Brownridge: I would like to
thank you both for encouraging me to continue this path and helping me through the process of
my thesis. To my grandmother, Geraldine Thomas: Thank you for checking in on me and
making home-cooked meals. To the inner city kids of Chicago, I am you and you are me. There
is nothing you can’t do if you put your mind to it and work hard. You can do this too!

iv

ACKNOWLEDGEMENTS

I would like to thank Dr. Mitra for giving me the guidance needed to conduct my
research. Also, I would like to thank my committee, Dr. Dong and Dr. Pierre, for their assistance
in my process. The Sloan Program: Thank you for not only the financial support but the
academic and social support as well. Dr. Percy Pierre, thank you for supporting me well before
entering Michigan State University and throughout. Your counseling and guidance is valuable to
the matriculation of students in the Sloan Program.

v

TABLE OF CONTENTS

LIST OF TABLES...………………………………………………………………………...…viii
LIST OF FIGURES..……………………………………………………………………………ix
Chapter 1 Introduction……………………….…………………………………………………1
Chapter 2 Background…….…………………………………………………………………….4
2.1 Wind Energy…………………………………………………………………….…….4
2.1.1 Characteristics of Wind……………………………………………………...5
2.1.2 Reliability Concerns…………………………………………………………5
2.1.3 Solutions…………………………………………………………………….7
Chapter 3 Different Techniques of Estimating Capacity Value….……………………….…..9
3.1 Capacity Value…………………………………………………………………….......9
3.1.1 Garver Approximation……………...………………………….…………....9
3.1.2 Multi-State Unit Representation……………………………….…………..10
3.1.3 Annual Peak Calculation…………………………………………………..10
3.1.4 Peak-Period Capacity Factor………………………………………………11
3.1.5 Z-Statistic Method………………………………………………………....11
3.2 Preferred Method……………………………………………………….……………12
Chapter 4 Energy Storage Systems….………………………………………………………...14
4.1 Compressed Air Energy Storage………………………………………….………….14
4.2 Flywheel Energy Storage………………………………………………….…………15
4.3 Superconducting Magnetic Energy Storage……………………………….………....16
4.4 Battery Energy Storage………………………………………....………….………...17
4.4.1 Lead-Acid Battery……………………………………….………....17
4.4.2 Sodium-sulfur (NaS) Battery…………………………….………...18
4.4.3 Zinc/Bromine (Za/Br) Battery………………………….………….18
4.4.4 Lithium-ion Battery……………………………………….……….19
4.5 Supercapacitor Energy Storage..…………………………………………….…….....19
4.6 Pumped Hydroelectric Storage…………………………………………….………...20
4.7 Application………………………………………………………………….……….21
4.8 Summary…………………………………………………………………….……….21
Chapter 5 Reliability..……………………………………………………………….………….24
5.1 System Description……………………………………………………….………….24
5.2 Discrete Convolution……………………………………………………….………..25
5.2.1 Generation Model……………………………………………….…………25
5.2.2 Load Model……………………………………………………….………..26
5.2.3 Generation Reserve Model………………………………………….……..28
5.3 Monte Carlo………………………………………………………………….………28
5.3.1 Mathematical Model……………………………………………….………29
vi

5.3.2 Sequential Method………………………………………………….……...29
5.4 Results and Discussion……………………………………………………….……...31
Chapter 6 Conclusions and Recommendations…………………….………………………...39
APPENDICES………………………………………………….……………………….……....41
Appendix A: Matlab Code- Discrete Convolution of Original System…………….........42
Appendix B: Matlab Code- Discrete Convolution of System with Wind Power………..45
Appendix C: Matlab Code- Discrete Convolution Iterations to Find LPI……..………...48
Appendix D: Matlab Code- Discrete Convolution of Wind Power and ESS…………....51
Appendix E: Matlab Code- Discrete Convolution Iterations of Wind Power and ESS to
Find LPI Increase………………………………………………………………………...55
Appendix F: Matlab Code- Monte Carlo of Original System…………………………...60
Appendix G: Matlab Code-Monte Carlo of System with Wind Power………………….62
Appendix H: Matlab Code-Monte Carlo Iterations to Find LPI………………….……...64
Appendix I: Matlab Code-Monte Carlo of Wind Power and ESS..………...…..………..67
Appendix J: Matlab Code- Monte Carlo Iterations of Wind Power with ESS to Find LPI
Increase………………….………………….………………….………………….……..71
Appendix K: Test System Data..…………………………………..……..…….……..….75
REFERENCES………………………………………………………………………….….…...76

vii

LIST OF TABLES

Table 4.1: Advantages and Disadvantages of Energy Storage Systems………………………22
Table 4.2: Comparison of Technical Aspects Energy Storage Systems………………………23
Table 5.1: Initialization of COPAFT……………………………………………………..…...26
Table 5.2: LOLE of Original and Wind Integrated Systems………………………………….31
Table 5.3: Factor Increases and Capacity Values……………………………………………..32
Table 5.4: Discrete Convolution Results for ESSs Discharging during All Peak Hours……..33
Table 5.5: Monte Carlo Results for ESSs during All Peak Hours…………………………….33
Table 5.6: Discrete Convolution Results for ESSs Discharging during a Portion of Peak
Hours………………………………………………………………….……………….34
Table 5.7: Monte Carlo Results for ESSs Discharging during a Portion of Peak Hours….......34
Table A.1: Generating Unit Reliability Data………………………………………………….86

viii

LIST OF FIGURES

Figure 2.1: Primary Energy Consumption with Renewable Energy Breakdown…………….....8
Figure 5.1: HLOLE for ESS Discharging During All Peak Hours………………………….....34
Figure 5.2: HLOLE for ESS Discharging During a Portion of Peak Hours…………………...35
Figure 5.3: Capacity Value Comparison for Discrete Convolution…………………………...36
Figure 5.4: Capacity Value Comparison for Monte Carlo………..…………………………...36
Figure A.1: IEEE RTS Load Data……………………………………………………………..86

ix

Chapter 1
Introduction

As fossil fuel prices and environmental pollution increases, countries have invested many
resources into wind energy. Since the 1990s, the annual growth rate of wind energy has exceeded
26% and many countries have set goals for high penetration level on the grid [2]. Wind energy is
the most mature renewable energy source and competitively priced compared to traditional
power plants. Wind energy does not produce any greenhouse gases and is located mostly in rural
areas or offshore. However, sites for wind farms tend to be in remote locations, resulting in
significant infrastructure development, especially transmission capacity. In addition, there are
concerns about noise, aesthetic aspects, and the harm to birds, due to the rotor blades. These
concerns have been mitigated through improvements of proper siting for wind farms. The major
challenge of wind energy, especially at high level of penetration, is the variability of wind and
load demand. Also, since wind energy is unpredictable, the uncertainty makes it more difficult to
operate the power grid. A solution to this challenge is integrating energy storage systems with
wind energy to alleviate the uncertainty and reliability issues.
One of the challenges imposed by increasing penetration of wind generation is balancing
the system load demand and available generation. When integrating wind generation with energy
storage systems, the mismatch between intermittent wind generation and the time varying load
demand is reduced [3].For capacity planning and reliability purposes, the capacity contribution
of an integrated energy storage system (ESS) with a wind generation system needs to be
estimated. A method to determine the capacity contribution is calculating the capacity value.

1

Capacity value can be defined as the amount of additional load that can be served due to the
addition of the generator, while maintaining the existing levels of reliability [1].
In relative works, techniques have been developed to estimate the capacity value of wind.
In [1, 8], there is a thorough description of different techniques to estimate the capacity of wind
and the preferred method. In addition, there are several papers discussing the use of ESSs to
mitigate the issues of high penetrations of wind power in the bulk power system. In [2], a battery
model is developed to demonstrate that energy storage technology can reduce the fluctuation of
wind power output. From there, researchers have analyzed more specific energy storage options
to combine with wind farms. In [22-23], proposed using superconducting magnetic energy
storage (SMES) to improve the power quality and stability of wind farms. In [21], uses SMES as
a solution to minimize the frequency fluctuation for isolated power system with wind farms. In
[24], using flywheel energy storage system for frequency control of an isolated power system
with wind farms. Now that energy storage technology has been analyzed, how to model the
usage of these systems are critiqued to improving wind power. In [7], an availability model is
developed to take in account the impact of outages during charging periods and the unavailability
during the discharging periods. In addition, it compares the difference between using ESSs for
reduction of variability and reliability improvement.
The motivation for this thesis is to investigate the different possibilities of ESSs and
estimate the capacity value of wind generation with an ESS. Capacity value has been heavily
researched with many methodologies but there is a preferred method. This preferred method is
combined with reliability approaches to produce the capacity value of the system. This
information can be used for the advancement of high penetrations of wind in the future.

2

Chapter 2 discusses wind power, characteristics of wind, reliability concerns, and current
solutions. Chapter 3 describes different methods used to estimate the capacity value and the
preferred method. Chapter 4 discusses multiple types of ESSs and their specifications. Within
chapter 4, there is a descriptive analysis of each energy storage system’s functionality,
efficiency, maturity, advantages, and disadvantages, charging capabilities, and discharging
capabilities. Chapter 5 describes two approaches to achieving the reliability of a system to lead to
the capacity value. Furthermore, a comparison of how the capacity value improves from a system
with wind power and a system with wind power and storage. The final chapter will include
further remarks and suggestions for future research.

3

Chapter 2
Background

The current U.S. energy portfolio relies heavily on fossil fuel resources. Due to the
Environmental Protection Agency (EPA) regulations and the depletion of fossil fuel reserves, the
energy portfolio needs to be diversified to handle the current and increasing energy demands.
Currently, natural gas has been an immediate option because the pricing of Shale gas has
decreased and it is cleaner than fossil fuel. With the advancement of new technologies,
renewable energy sources will be the solution to have an environmental-friendly and independent
energy portfolio. About 10% of energy supply is renewable energy sources. Renewable energy is
an energy source that is naturally replenished or inexhaustible such as wind, solar, hydro,
biomass, and many others. Figure 2.1 shows the distribution of different types of renewable
energy resources in the complete primary energy usage.
Wind and solar are the two most rapidly growing renewable energy sources. Globally,
these are 283 GW of installed capacity for wind and 100 GW of installed capacity for solar [25].
However, solar is growing at a higher rate than wind. For this thesis, wind will be the focused
renewable energy source.

2.1

Wind Energy
Wind is the movement of air across the surface of the Earth [2]. Wind energy is a

renewable source of power that uses turbines to harvest the kinetic energy of wind to be
converted to electricity. The minimum speed to produce power is 11km/hour and the optimum

4

wind speed is 60 km/hour [21]. When wind speeds are 90km/hour or greater, the main computer
shuts down as a safety precaution. Power can be calculated by the following equation:

(2.1)
where
swept area of a wind turbine
density of air
wind speed
power coefficient.

2.1.1 Characteristics of Wind
Wind power brings up many environmental and reliability concerns. From an
environmental view, wind is a clean and resourceful way to produce power versus the
conventional way. However, certain interest groups are concerned with the amount of area
needed to commission a wind farm, the aesthetics and noise of the wind turbines, and the habitat
for birds in the area. From a reliability view, wind energy poses several problems for planning
operations such as variability, frequency fluctuation, increase ramping rate and range, and others.
There will be further discussion of these concerns to follow in the next section.

2.1.2 Reliability Concerns
Since wind power is based on wind speed, there are times when there is no generation
due to low wind speeds. Earlier in the text, it stated the wind speed needs to be a minimum of 11
km/hour to start generating power and cannot exceed 90 km/hour for safety reasons. Wind is not

5

predictable and is variable because wind does not blow all the time or meet those specifications
to generate power. This is an issue planning operations have because a plan is develop by the
load and generation to figure out how to meet demands. If power generation is changing
spastically then a proper plan cannot be developed or executed. In addition, wind generates the
most power during the off peak hours of the day, which is completely opposite to the daily load
curve.
The variability of wind power leads to other reliability concerns such as frequency and
power fluctuation. The inconsistent wind speeds create output power fluctuation, which causes
network frequency and voltage fluctuation in the system. Therefore, the power quality of the
system declines. The system’s inertia plays an important role in determining the rate of
frequency change after an event. Renewable energy sources such as wind, contribute little to no
inertia to the system. If wind power continues to increase, frequency control will be difficult in
the future. For isolated power systems, which benefit from having wind power, are lacking in
capabilities of power regulation. The utilization of wind power is more effective in an isolated
power system because the main power is driven by a diesel engine which is bad for the
environment [24].
Other reliability concerns are ramping rate and range. A ramping event is the change of
power at every time interval. If the power change is positive, it is a ramp-up event and if
negative, a ramp-down event. The ramp rate is the power difference from time interval to time
interval. When integrating wind with the system load, the wind generation is considered a
negative load and summed to create the net load. The net load maxima are shifted down and the
local minima are shifted down and stretched. Depending on the amount of wind power being

6

produced, the range increases significantly. This mainly occurs when the load is at its local
minima and the wind power is at its local maxima.
2.1.3 Solutions
The two predominant solutions to the reliability concerns of wind power fall in the
spectrum of geospatial and temporal. A geospatial solution involves an analysis of the entire
country for wind power output, transmission capabilities, and load demands. For areas that
produce high wind power and a low load demand, would have excess power. This power can be
transferred to other areas with a higher load demand. However, this solution is not plausible
without building more transmission in the current system infrastructure. The cost to update the
infrastructure with new transmission lines to meet the target of 20% of wind by 2030 is about 20
billion dollars [26]. This is a long term and expensive solution to the increasing energy demands.
The other option is a temporal solution for mitigating the variability of wind energy with ESSs.
The solution is an immediate compared to upgrading the system infrastructure which take
significant time and resources. However, ESSs are also expensive due to the material, size, and
sometimes installation. In addition, ESSs are developing as technology changes and some
systems are in development or prototype stages. Proper planning and decision making based on
the system needs would determine what are the best ESSs needed to be use. The specifications
these ESSs need to meet to be beneficial are size capacity, charging/discharging, power density,
life cycle, and a few others. The next chapter will discuss different ESSs and how they compare
with these specifications.

7

Figure 2.1: Primary Energy Consumption with Renewable Energy Breakdown [27]

8

Chapter 3
Different Techniques of Estimating Capacity Value

3.1 Capacity Value
Capacity value is of great importance to wind power since it assesses the risk of
generation capacity deficiency [1]. The metrics that are used for evaluation are loss of load
probability (LOLP) and loss of load expectation (LOLE). LOLP is the probability of when the
load exceeds the available generation but does not explain for how long. LOLE is the expected
number of hours or days the load will not be met over a period of time. The effective load
carrying capacity (ELCC) is the metric used to find the capacity value [4]. There are numerous
methodologies to finding the capacity value of wind generation such as: Garver Approximation,
Multi-State Unit Representation, Annual Peak Calculations, Peak Period Capacity Factors, and
Z-Statistic Method.

3.1.1 Garver Approximation
Garver Approximation has been important in the calculation of capacity values but is
outdated by the advancements in computing power [1]. Garver proposed using the graphical
approach when calculating the capacity value of an additional generator. This same methodology
can be used to estimation the ELCC of a wind generator added to a power system. The main two
assumptions of Garver’s methodology are:


The probability distribution for wind availability is the same at all times.



LOLE is calculated as

, where

is the peak demand, and

parameters.
9

and

are fitting

The ELCC of wind generation is calculated by the following:

(3.1)

where

is the probability that the available wind capacity is

[1].

3.1.2 Multi-state Unit Representation
Multi-state unit representation uses probabilistic representation of wind plants similar to
conventional units. The wind plant is model with partial capacity outage states with
corresponding probabilities. When calculating the LOLP, the wind generation is included in the
capacity outage probabilities table (COPT) calculation. The wind power is typically model by a
histogram of wind power output during a period of time. A major concern with this approach is
the loss of wind/load correlation [1]. There are factors that affect the wind energy availability
such as, seasonal and diurnal variation, which cannot be modeled in one probability density
function. However, this concern has been address by using different probability distributions at
categorized hours.

3.1.3 Annual Peak Calculations
Annual peak calculations use the LOLP at the time of annual peak demand to represent
the system’s risk. The calculations are the same as year-around risk calculations besides using
the LOLP at the time of annual peak. Peak demand only occurs once a year which is problematic
when choosing a probability distribution for wind capacity. The two approaches to evaluate wind
resources at annual peak are:

10



Use a histogram to represent the hourly load for the entire peak season. However, this
approach is limited, due to many days not being close to the annual peak demand.



Use a histogram to represent the load factors from hours where demand is a percentage of
the year’s peak. This has greater relevance to the peak demand and reduces the amount of
data.

The major disadvantages of the annual peak calculations are that it does not consider the loss of
load at other times of the year, and it is difficult to designate appropriate distributions for wind
resources, as well as peak load [1].

3.1.4 Peak-Period Capacity Factor
Peak-period capacity factor calculates the capacity factor over peak periods to estimate
the capacity value of wind. Using this approximation has been fairly accurate with the exception
of hydro and transaction schedules [5]. This is due to hydro and transaction schedules being
positively correlated to the load. However, capacity factors do not provide the short term or
annual variability of wind power or the correlation to demand. Capacity factors gives insight to
the capacity value of that system such as a high capacity factor yields a high capacity value for
the same system.

3.1.5 Z-Statistic Method
The z-statistic method is taking the difference between available resources and load
during peak hours as a random variable with a probability distribution. The z-statistic is kept
constant by adding the load carrying capacity for an added power plant in increments. When

11

keeping the z-statistic constant, it is equivalent to keeping the LOLP constant. There are two
following assumptions involved in this approach:


When adding wind, the shape of the probability distribution for the margin of available
capacity does not change significantly. However, the mean and standard deviation may
change.



The standard deviation for available wind capacity is small versus the available capacity
from existing generation. Therefore, the z-statistical method is only valid for small
penetrations of wind power.

From these assumptions, the expression for ELCC is derived as:

(3.2)

where
mean wind load factor over peak hours
z-statistic representing LOLP.
This method provides a greater insight on what influences the level of the ELCC than iterative
calculations [1]. In addition, this method is relatively simple for determining how wind
variability and correlations among wind power affects the load caring capacity.

3.2

Preferred Method
The preferred method is modeled after calculations for conventional thermal generation

in power systems. They are modeled by their perspective forced outage rates (FOR) and
capacities. This information is convolved through an iterative method to produce the COPT. The
COPT is comprised of the capacity levels and their respective probabilities [6]. However, wind
12

cannot be modeled by capacity and FOR because of its variability. For calculating the ELCC
with wind, the following three steps can be used:
1. Compute the LOLPs of a system without wind using the hourly load time series and
the COPT of the system. The annual LOLE should be calculated and meet the
reliability target for that period. However, if the target is not met, reduce the load to
achieve the reliability target.
2. The time series of wind power is treated like a negative load and combined with load
time series. The LOLE is recalculated and should be lower (better) than in step 1.
3. The load data is increased by a ΔL across all hours using an iterative process until the
LOLE reaches the LOLE in step two. The increase in load that achieves the target
reliability is the ELCC/capacity value of the wind plant.

This chapter has provided different studies of capacity value which explain the different
levels of capacity value. However, there are factors that influence the calculation of capacity
value that these methods have in common. The key factor in calculating the capacity value is to
capture the relationship between wind and the load. When modeling the relationship between
wind and the load, it is critical to use information from the same year [1]. In addition, for
variable generation, it is best practice to use a year or more data to help with prediction of wind’s
capacity value.

13

Chapter 4
Energy Storage Systems

Energy storage can be divided into two categories: direct and indirect energy storage.
Direct energy storage is a storage system that stores electric energy in an electric or magnetic
field. Indirect storage is a storage system that converts electric energy to kinetic, chemical or
mechanical energy. Depending on the power system’s parameters, these would determine what
ESS is suitable. There are many ESSs such as, compressed air, flywheel, battery,
supercapacitors, superconducting magnetic, pumped hydro, and many more. The next sections
will provide a review of these types of ESSs.

4.1

Compressed Air Energy Storage
Compressed air energy storage (CAES) is storing energy by compressing and cooling air

at high pressures in a sealed geological vessel. When needed, the compressed air is combined
with natural gas and used to rotate a turbine for generating electricity. The efficiency of a CAES
system is comprised of two components: energy ratio and heat rate [9]. The energy ratio is
defined as:
(4.1)

where
Egen

energy out from generation; and

Ecomp

energy consumed by the compressor.

The efficiency of the system is defined as:
14

(4.2)
where
HR

amount of gas consumed per unit of energy generated.

There are different types of CAES systems according to generation: first, second, and
third generation. The first generation is the only CAES system in conventional use. This system
is composed of compression and generation components that mix natural gas with compressed
air and burned as the same in conventional turbine plants. The second generation systems are
similar to the first but the design is coupled with technological advancements. The second
generation system has a turnaround efficiency of about 54% compared with 48-50% with first
generation [9]. The third generation systems do not use natural gas but instead stores heat during
compression to be used for heating the compressed air later. This system is still in research and
development stage due to technology issues with heat storage. The efficiency of this particular
system is directly connected to the energy ratio since gas turbines are not used. A great benefit is
zero carbon emissions and no fuel consumption. However, the downside is the increased capital
cost and technology immaturity. Since the second generation system is more efficient than the
first and the third generation is in research and development, investments towards CAES systems
are for a second generation hybrid system [9].

4.2

Flywheel Energy Storage
Flywheel energy storage stores kinetic energy of a rotating mass in an electromechanical

storage system. There are major two types of flywheel system: steel- and composite-rotor
flywheels. In both systems, the energy is stored in the momentum of the rotating rotor. In a
15

vacuum, the rotor spins on the bearings to reduce friction and increase efficiency. This energy
can be converted to electrical or mechanical forms due to the rotor containing a motor/generator
[10]. Steel-rotor systems depend on the mass of the rotor for energy storage while composite
systems depend on speed. When charging the system, the current flows through the motor which
increase the speed of the flywheel. When the system is discharging, current flow is produced by
the generator which slows down the wheel. The sizes range from 40kW to 1.6MW for times of
5-120 seconds [10]. The advantages of a flywheel energy system are low-life cycle costs due to
minimal maintenance, compact and can be arrange in small area, and contain no harmful
chemicals or produce any gases [10]. A disadvantage is the containment of all the rotating energy
equipment in the system. To minimize the issue, a particular site, design, and material have to be
chosen or the system must be tested and rated.

4.3

Superconducting Magnetic Energy Storage
Superconducting magnetic energy storage (SMES) stores electrical energy in a magnetic

field without converting to chemical or mechanical form. In SMES, there is a coil made of
superconducting wires that allows direct current to flow through without any significant loss.
This produces a magnetic field which stores the energy. SMES have to operate in cryogenic
temperatures to maintain superconductivity. The efficiency of an SMES system is dependent on
the duration of which it is stored. Thus, the efficiency of the system is defined as:
(4.3)
where
0.025, charge/discharge inefficiency

16

0.04, daily loss
number of days between charge and discharge
Typically, a SMES system will not charge and discharge completely. Therefore,

will not be an

integer. The expected range of efficiency for SMES for long turn use is 80 and 90 percent and
about 98% for short-term use [11]. Most SMES units fielded to date provide 1MW for 1 second
and can be paralleled for more power [10]. The SMES coil has the ability to release large
quantities of power and fully recharge in minutes. There are numerous advantages of a SMES
system such as controllability and reliability, no degradation in performance overtime, and very
economical. The system is compacted, self-contained, highly mobile, and can be kept in remote
locations. There are no harmful chemical or produces any flammable gases and the estimated life
expectancy of a typical system is at least 20 years [10].

4.4

Battery Energy Storage
Battery energy storage (BES) stores electrical energy through the reverse chemical

reactions of a typical battery. There are many types of BES systems such as lead-acid, sodiumsulfur, zinc/bromine, lithium-ion, and many more.

4.4.1 Lead-Acid Battery
Lead-acid batteries have been in industry use for over 20 years. They are traditionally
used as form of backup power and in automobiles [12]. The lifetime of lead –acid batteries
depend on the usage, discharge rate, and number of deep discharge cycles. When charging, at the
negative electrode, an electrochemical reaction is converting lead sulfate into a soft lead. Also, at
the positive electrode, the lead oxide is formed due to the current flow. One main advantage of
17

lead-acid batteries is the low cost. The cost of lead-acid batteries is mainly influenced by the
price of lead at the time. Lead-acid batteries have a relatively short calendar/cycle life of
approximately 4-5 years/ 750 cycles, meaning for heavily usage, they would need to be replaced
frequently [12]. Therefore, they are mainly used for power quality management and emergency
power.

4.4.2 Sodium-sulfur (NaS) Battery
Sodium-sulfur (NaS) batteries contain sodium at the negative electrode and sulfur at the
positive electrode. NaS batteries are used for a wide variety of applications such as, peak
shaving, emergency power, power quality management, and peak shaving [12]. Some advantages
of NaS batteries are long cycle life, low material cost, temperature stability, and high power and
energy density. The energy density is three times as much as tradition lead-acid batteries [12].
The efficiency of DC conversion for NaS batteries is approximately 85% which make them ideal
for future DC distribution systems [12].

4.4.3 Zinc/Bromine (Za/Br) Battery
Zinc/Bromine (Za/Br) batteries are a type of flow battery that uses a pump system to
circulate reactants throughout the battery. The battery modules are rated to discharge at 150A at
an average voltage of 96V for 4 hours [12]. Since the battery is stacked in parallel, individual
stacks can be replaced without replacing the whole system. The life cycle of Za/Br batteries are
rated at 2500 cycles and uses less toxic electrolytes than traditional lead-acid batteries [12]. Due
to the lack of technological maturity, there are very few real world installations.

18

4.4.5 Lithium-ion Battery
Lithium-ion batteries allow lithium ions to flow between the cathode and anode
producing a current flow. There are many chemical configurations and this battery is one of the
newer technologies. The advantages of lithium-ion batteries are long calendar life, high energy
density, low self-discharge, and no memory effect. However, lithium is very scarce which makes
them more expensive than others [13].

4.5

Supercapacitor Energy Storage
Supercapacitor energy storage stores electric energy through the double layer effect. A

supercapacitor consists of an ion permeable membrane that separates two electrodes and the
electrodes are electrically connected by an electrolyte. When voltage is applied, an electric
double layer is formed at both electrodes. The charges are separated by the double layer which
has no conventional solid dielectric. The capacitance values are determined by two different
storage principles: double-layer and pseudocapacitance. In double-layer capacitance, the
separation of charge by a Helmholtz double layer, stores electric energy electrostatically [14]. In
pseudocapacitance, there is a faradaic redox reaction that electrochemically stores the electric
energy. The efficiency of a supercapacitor is comprised of the following: charging current, total
energy consumed during charging, and energy stored. The charging current is defined as:
(4.4)

where
C

ideal capacitor.

The total energy stored consumed during the charging period is defined as:

19

.

(4.5)

The energy stored of the supercapacitor is defined as:
(4.6)
The charging efficiency is defined as:
.

(4.7)

In many ways a supercapacitor is incomparable to a battery, in that, it’s light weight, great
capacity, higher power density, faster charging rate, longer shelf and cycle life, and wider
temperature characteristic [14].

4.6

Pumped Hydroelectric Storage
Pumped hydroelectric storage (PHS) stores energy using gravitational potential energy of

water, where water is pumped from a lower elevated reservoir to a higher one. During low
demand, energy is used to pumped water to the high reservoir. When demand is high, the water
is release back to the lower reservoir through a turbine to generate electricity. The roundtrip
efficiency of PES for old units are 60% and for newer unit as high as 78% [15]. PHS units are
very flexible and will start up from 200 to 1600 times each year [15]. The time to charge to full
generation ranges from 5 to 22 minutes and a cold start takes approximately one to four minutes.
The disadvantages are the environmental impact, site location due to the size, and investment
cost versus traditional combustion turbines. Currently, PES is the most favorable and proven to
be combined with wind resources to optimize operating benefits and meet required supply level
[15].

20

4.7

Application
The purpose of using the ESS is to mitigate the reliability issues discussed in chapter 2.

Based on the information given, the ESS needs to discharge quickly to help reduce the ramping
range and rate. Discharging the stored energy will help level the minima and maxima of wind
energy to a more consistent output. If maturity was not a factor, flywheels, superconductors,
batteries, and SMES would be suitable options. In addition, appropriate is needed based on the
installed capacity. ESSs such as CAES or PHS are not suitable options because of their
geological restrictions and the discharging time frame. Even though PHS is the only grid-scale
storage technology, the sizing of PHS is not appropriate for wind farms. The difference in scale
would not provide the advantages of integrating with wind.
4.8

Summary
Table 1 presents the advantages and disadvantages of each system. Table 2 discusses

some of the important parameters to determining a suitable system for a power system based on
cost, efficiency, nominal power and others.

21

Table 4.1: Advantages and Disadvantages of Energy Storage Systems [17]
Systems

Battery Technologies

CAES

Advantages
-High power density
-High energy density
-Mature Technology

Flywheel

-High power density
-Long cycle life
-Very fast recharge

SMES

-High power density

Lead-Acid

-Mature technology
-Inexpensive

Soduim-Sulphur
(NaS)

-High power density
-High energy density
-Mature technology
-High efficiency

Zinc-Bromine
(Zn-Br)

Lithium-Ion

Supercapacitor

Pumped Hydroelectric

-High power density
-High energy density

-High power density
-High energy density
-High efficiency
-High power density
-Long cycle life
-Very fast recharge
-High power density
-High energy density
-Mature technology

22

Disadvantages
-Geographical limitations
-Expensive to site and build
-Requires fuel input
-Only for huge applications
-Low energy density
-Large standby losses
-Low energy density
-Expensive technologies
-High auxiliary needs
-Low power density
-Low energy density
-Short life cycle
-Environmental hazard
-High maintenance
requirements
-Expensive technology

-Early stage technology
-Toxic components
-Fast corrodible components
-High maintenance cost
-Early stage technology
-Expensive technology
-Low energy density
-Requires advanced power
electronics
Expensive technology
-Geographical limitations
-Expensive to site and build
-Requires fuel input
-Only for huge applications

Table 4.2: Comparison of Technical Aspects Energy Storage Systems [16, 17, 18]
Systems

Maturity

Unit Costs

Efficiency
(charging/
discharging)

CAES

Mature

700 $/kW

60%-70%

1300 $/kW

90%

0.01-2
MW

secs-mins

700 $/kW4000$/kW

90-95%

10-100
MW

secs-mins

1100 $/kW3100 $/kW

70-90%

0.001-10
MW

mins-days

Flywheel
SMES
Battery

Mature and
commercial
products
Prototypes
units
Mature and
commercial
products

Typical
Nominal
power
0.1-400
MW

Discharging
at Nominal
Power
hrs-days

Supercapacitor

Prototypes

$30 (per cell)

80-95%

0.01-1
MW

secs-mins

Pumped
Hydroelectric

Mature

2300 $/kW

80%

0.1-1 GW

hrs-days

23

Chapter 5
Reliability

5.1

System Description
The test system examined in this thesis is the IEEE RTS-79 [20] with a peak demand of

2850 MW and a total installed generation of 3405 MW. The wind speed data has been collected
by the government of Canada records for historical climate data. The data was recorded for
providence Alberta, Canada and the assumption is that this data would generate the amount of
power as Halkirk Wind Farm. Halkirk wind farm is a 150 MW farm in Halkirk, Alberta, Canada.
The parameters of Halkirk Wind Farm are used to determine the wind power. The minimum
wind speed for Halkirk’s wind turbines to produce power is 11 km/hr. The optimum wind power
speed for power production is 60 km/hr. If winds are 90 km/hr or greater, sensors signal the
main computer to power down for safety reasons. These parameters are converted to m/sec to be
use in equation 1.1.
The charging and discharging of the energy storage system (ESS) is based on the peak
and off-peak hours. Since wind produces the most power during the morning and the late
evening, the system will be charging after 5pm until 8am. During peak hours, more power is
consumed and the ESS can optimize the wind power by discharging during peak hours. The
storage volume is an integer multiple of the installed capacity which determines when the ESS
will be exhausted. ESSs of various installed capacities are combined with the wind power to
determine the impact on the capacity value. The wind power and energy storage are treated as a
negative load and combined with the system load. The ESS can only be fully charged if the
amount of wind power produced is equal to or greater than the ESS’s volume.

24

Reliability can be defined as the probability that a component or system will perform its
designated functions for a given period of time under the conditions in which it was designed to
operate. In this thesis, two approaches will be used to find the reliability of a system without
wind power, with wind power, and wind power with storage. The reliability index (LOLE) will
yield the capacity value to be compared amongst the systems.

5.2

Discrete Convolution
Discrete convolution consists of building a generation and load model containing

probability and frequency distributions to represent the total generation and load. The generation
model is constructed with the capacity outage levels better known as the capacity outage
probability and frequency table (COPAFT). The load model is constructed with the probability
and frequency distribution of the total system load. These two models together are used to build
the generation reserve model. The generation reserve model represents the probability and
frequency distributions for excess system generation. This information is obtained by performing
discrete convolution of the generation and load model.

5.2.1 Generation Model
The generation model is constructed using the following equation (5.1) and (5.2)
(5.1)

(5.2)
where
u

repair rate of the new unit

25

p

availability probability

q

failure probability

The variable X represents the outage state and subscript i is the index of the existing capacity
outage state Ci such that Cj=X. If a state does not exist in the “old” COPAFT, then Ci is the
smallest of the existing states that are larger than X. The subscript j is the index of the existing
capacity out state Cj such that Cj=X-CN. If a state does not exist, then Cj is the smallest of the
existing states that are larger than X-CN.

Table 5.1: Initialization of COPAFT
State
1
2
3
4
5
6
7
8
9
10
11
12
13

Co=Ci (MW)
0
1
2
3
4
5
6
7
8
9
10
11
12

Pi=P{Co>=Ci}
1
0.02
0.02
0.02
0.02
0.02
0.02
0.02
0.02
0.02
0.02
0.02
0.02

Fi=F{Co>=Ci}
0
0.000333333
0.000333333
0.000333333
0.000333333
0.000333333
0.000333333
0.000333333
0.000333333
0.000333333
0.000333333
0.000333333
0.000333333

Table 3.1 represents the blank COPAFT to start with where, q=0.02, f=qu=0.000333333, for the
first 12 MW generator. There are 32 generators in the system and therefore, the COPAFT needs
to be updated 32 times.

5.2.2 Load Model

26

The load model is constructed by the following two equations (5.3) and (5.4).
(5.3)

(5.4)
For the generation reserve, most of the work was done doing the generation system and load
model. To compute the probability and frequency the following equations (5.5) through (5.11)
were used. M can be expressed as:
(5.5)
Cc is the total generation capacity of the system; Co is the forced outage capacity; L is the load
capacity. The probability of M can be expressed as:
(5.6)
where
(5.7)
(5.8)
The frequency of M can be expressed as:
(5.9)
where
(5.10)
(5.11)

The maximum load is 2850 from the data sheet and since it does not get any greater the
probability and frequency for the load will be zero until it reaches 2850.
27

5.2.3 Generation Reserve Model
From generation reserve model, the system reliability index is computed and this will give
the hourly LOLE to be compared with the other systems.

HLOLE=LOLP x D

(5.12)

where

5.3

LOLP

loss of load probability

D

length of time interval in hours

Monte Carlo
Monte Carlo consists of imitating a physical system with a stochastic behavior and is an

alternative to analytical methods [19]. When systems are too complex to solve with analytical
methods, simulation such as Monte Carlo are used. The components of a system are divided into
elements whose behavior can be predicted either deterministically or probability distributions.
Therefore, mathematical models are used and preceded by performing sampling experiments.
Monte Carlo simulation methods can be classified as sequential or non-sequential. Sequential
simulation methods replicate the system behavior over time [19]. The mathematical model
generates the artificial history to be analyzed over time and a statistical conclusion is drawn.
Non-sequential simulation methods perform random sampling over all states that can be assumed
during the period of interest. After an appropriate number of samples have been drawn, a
statistical conclusion can be drawn. For this thesis, sequential methods will be used for the
system described in section 3.1.

28

5.3.1 Mathematical Model
The IEEE-79 RTS has 32 generators, which yields 232 states. Using a vector of 32
dimensions to represent the generators’ states, 1 represents the up state and 0 represents the down
state in the following equality.

Geni =

Drawing 32 random numbers for each generator

, for every hour, time

is calculated using

equation 5.13

(5.13)

When generator i is in the up state, p is the failure rate of that generator. If generator i is in the
down state p is the repair of that generator. The generator vector is updated by choosing the
smallest t amongst the 32 times and changes the state of the generator that corresponds to that
time. The interval of the system load is every hour and the simulation must track the amount of
time passed. When the generator vector is updated, a new generator capacity is calculated. At
that particular hour, the load capacity can be determined as well. A comparison between the
generator capacity and load capacity is used to determine if the entire system has failed. The
reliability indices are updated accordingly.
5.3.2 Sequential Methods
Sequential methods can be divided into three categories: synchronous timing, asynchronous
timing, and mixed timing. Synchronous timing method is also known as fixed time interval
29

method [19]. It is a two-step method and a time interval Δt is determined by the system’s
characteristics. In the initial state, time is increased by Δt and the simulation determines what
event has occurred. Asynchronous timing method, also known as next event method, is advanced
by a variable amount of time versus fixed. The simulation keeps a record of the next couple of
simulated events to occur. Once the imminent event has occurred, the time is changed to the
occurrence of that event. Mixed timing incorporates both synchronous and asynchronous timing.
This is the most commonly used in the industry for the purposes of both reliability evaluation
and production [19]. Typically, the form consists of hourly load curve over a period of time and
advancing states of system components asynchronously. For this thesis, the mixed timing method
will be used. There is an algorithm that describes a typical embodiment of this mixed timing
method [19]. The algorithm has been modified to fit the needs of this thesis as follows:
1. Input system data, including generation capacity, generation probability data, repair rate,
failure rate, initial state, load data and initial indices.
2. Initial sample path index at i=1
3. Draw random numbers for the system
4. From those random numbers, compute the time and find tmin.
5. Find the generator that corresponds to tmin, which is the imminent event, and change its
state. Find the total generation capacity.
6. Find the load during tmin and determine in each hour, if the system failed.
7. Update the failure time and compute hourly LOLE.
5.4

Results and Discussion

30

Using Discrete Convolution and Monte Carlo described above, the results of the IEEERTS 79 system is displayed in table 5.2. This data is used as a marker to compare the results of
other systems. Also, the results of the system with wind power is displayed which determines the
capacity value.

Table 5.2: LOLE of Original and Wind Integrated Systems
Method
Discrete Convolution
Monte Carlo

Original System (LOLE)
9.303
9.423

Wind Integrated System (LOLE)
7.165
7.270

In chapter 2.7, the preferred method was introduced to find the capacity value of wind
power. When finding the capacity value of wind, the LOLE should be lower (better) than the
original system. In the preferred method, the first step is to compute the LOLE of a system
without wind using the hourly load time series and the COPT of the system, which was
accomplished in the previous section. The next step is to calculate the capacity value of wind
power. The reliability methods will be repeated to complete the second step, where the wind
power is a negative load and combined with the system load. Once the new LOLE is acquired,
the wind power is removed. In the third step, the system load is increased by a small factor to
obtain the same LOLE as the first step. This is accomplished by performing iterations of the
system and comparing the LOLEs until the difference is a tolerance of 10

-4

for discrete
-2

convolution. Since Monte Carlo is based on randomness, the tolerance was increase to 50 . In
addition, a number of one million samples of random numbers were generated to produce the
results for the Monte Carlo method.

31

Table 5.3: Factor Increases and Capacity Values

Method
Discrete
Convolution
Monte Carlo

Load Factor Percentage
Increase (LFPI)

LOLE after LFPI

Capacity Value

1.472%

9.302

41.952 MW

2.110%

9.356

60.135 MW

The charging and discharging of the ESS is based on the peak and off-peak hours. Since
wind produces the most power during the morning and the late evening, the system will be
charging after 5pm until 8am. During peak hours, more power is consumed and the ESS can
optimize the wind power by discharging during peak hours. The storage volume is an integer
multiple of the installed capacity which determines when the ESS will be exhausted. ESSs of
various installed capacities are combined with the wind power to determine the impact on the
capacity value. The wind power and energy storage are treated as a negative load and combined
with the system load. The ESS can only be fully charged if the amount of wind power produced
is equal to or greater than the ESS’s volume. Table 5.4 and 5.5 displays storage volumes to
discharge power for all peak hours and table 5.6 and 5.7 displays storage volumes for a portion
of the peak hours. Figures 5.1 and 5.2 show the comparison of HLOLEs for all peak hours and a
portion of peak hours. Figures 5.3 and 5.4 show the comparison of capacity values for the system
without wind, with wind, and wind with ESS.

32

Table 5.4: Discrete Convolution Results for ESSs Discharging during All Peak Hours
Installed
Capacity
100kW
150kW
200kW
400kW
1MW
2MW
3MW
4MW
5MW
10MW

Storage
Volume
800kW
1.2MW
1.6MW
3.2MW
8MW
16MW
24MW
32MW
40MW
80MW

Factor
Increase
4.816%
4.816%
4.816%
4.784%
4.550%
3.473%
2.730%
2.367%
2.233%
1.733%

HLOLE
4.389
4.389
4.389
4.404
4.496
5.186
5.811
6.164
6.312
6.864

Capacity
Value
137.256
137.256
137.256
136.344
129.675
98.981
77.805
67.460
63.641
49.391

Percentage
Comparison
327%
327%
327%
325%
309%
254%
185%
161%
152%
118%

Table 5.5: Monte Carlo Results for ESSs during All Peak Hours
Installed
Capacity
100kW
150kW
200kW
400kW
1MW
2MW
3MW
4MW
5MW
10MW

Storage
Volume
800kW
1.2MW
1.6MW
3.2MW
8MW
16MW
24MW
32MW
40MW
80MW

Factor
Increase
7.0%
6.3%
6.4%
6.7%
5.6%
4.1%
2.6%
1.6%
0.9%
0.1%

HLOLE
4.655
4.872
4.840
4.827
4.786
4.967
5.257
5.608
5.865
6.240

33

Capacity
Value
199.50
179.55
182.40
190.95
159.60
116.85
74.10
45.6
25.65
2.85

Percentage
Comparison
332%
299%
303%
318%
265%
194%
123%
75.8%
42.7%
4.74%

Figure 5.1: HLOLE for ESS Discharging During All Peak Hours
(For interpretation of the references to color in this and all other figures, the reader is
referred to the electronic version of this thesis.)
Table 5.6: Discrete Convolution Results for ESSs Discharging During a Portion of Peak Hours
Installed
Capacity
100kW
150kW
200kW
400kW
1MW
2MW
3MW
4MW
5MW
10MW

Storage
Volume
500kW
750kW
1MW
2MW
5MW
10MW
15MW
20MW
25MW
50MW

Factor
Increase
3.139%
3.139%
3.139%
3.139%
3.051%
2.606%
1.933%
1.623%
1.489%
1.467%

HLOLE
5.542
5.542
5.542
5.542
5.610
5.986
6.713
7.020
7.147
7.162

34

Capacity
Value
89.462
89.462
89.462
89.462
86.953
74.271
55.091
46.256
42.437
41.809

Percentage
Comparison
113%
113%
113%
113%
107%
77%
31%
10%
1%
-0.4%

Table 5.7: Monte Carlo Results for ESSs Discharging During a Portion of Peak Hours
Installed
Capacity
100kW
150kW
200kW
400kW
1MW
2MW
3MW
4MW
5MW
10MW

Storage
Volume
500kW
750kW
1MW
2MW
5MW
10MW
15MW
20MW
25MW
50MW

Factor
Increase
4.0%
3.4%
3.1%
3.4%
3.5%
3.2%
2.1%
1.2%
0.8%
0.0%

HLOLE
6.325
6.853
6.837
6.716
6.642
6.586
6.660
6.758
6.947
7.122

Capacity
Value
114.00
96.90
88.35
96.90
99.75
91.20
59.85
34.20
22.80
00.00

Percentage
Comparison
190.0%
161.0%
147.0%
161.0%
166.0%
152.0%
99.5%
56.9%
37.9%
00.0%

Figure 5.2: HLOLE for ESS Discharging During a Portion of Peak Hours

35

Figure 5.3: Capacity Value Comparison for Discrete Convolution

Figure 5.4: Capacity Value Comparison for Monte Carlo
36

From the figures and tables above, the capacity values started at their highest and
significantly decreased. Most of the data in the beginning is similar due to the small increment
increase in installed capacity compared to the system load. There is not a significant change until
the installed capacity reaches 1 MW. As the volume size increases, the capacity value decreases
almost exponentially until it reaches its original capacity value of about 42 MW for discrete
convolution. However, Monte Carlo results show the capacity values decreasing beyond the
original capacity value of 60 MW. This data shows it is important to choose a volume size most
optimum for the wind farm’s generation to benefit from the ESSs. In addition, the distribution of
power during a select number of hours can affect the capacity value. The capacity values of the
ESS being discharged during the full peak hour period is higher than those acquired for a portion
of the peak hour period. The data results using Monte Carlo is different from discrete
convolution because of the randomness of the method and the higher tolerance difference.
However, the Monte Carlo method still shows the reliability and capacity value decreasing as the
ESS volume increases.
These two methods are well-known in the reliability world. Discrete convolution is the
more favorable method because it takes in account all the parameters of the system. However,
this method can get very time consuming and unmanageable the large the system becomes. In
this case, Monte Carlo can be used as an alternative method. This method uses randomness to
emulate the system. Between the two, discrete convolution is the more accurate method.
However, the more trials performed with Monte Carlo, the more accurate the results will be or
fairly close to discrete convolution. In this test system, the Monte Carlo simulation only
performed one million trials. To get closer results, the simulation may have needed to run a

37

hundred million trials or more. In spite of this, the proper programming language that has the
computing power could yield better results for the Monte Carlo method.

38

Chapter 6
Conclusion and Recommendations

The objective of this research was to evaluate the integration of ESS with wind power
and how that would improve the capacity value. This research shows that ESS is a solution to the
variability of wind power and the mismatch between load and generation. The results of the ESS
combined with wind improved the capacity value 2 to 3 times the capacity value of wind alone.
With ESSs discharging during peak hours, it responds to the load increase and provides a
consistent level of power. Peak hours are the most critical periods because of the spike in load
and it’s the period when wind farms are not generating as much power. Based on the test system,
ESSs that are mature, 10MW or less in size, discharges and charges fast, and cost effective
would be suitable storage options. Possible options using tables 3.1 and 3.2 could be flywheels
and batteries. However, more research in the area of Halkirk Wind Farm would have to be taken
into account to consider CAES or PHS because of their geological restrictions.
For recommendations, the data used in this study should be data from an actual wind
farm. In addition, multiple years of this data should be used as well. For this study, there was
only one year of load data and it is necessary to have equal time frames of load and wind speed
data. This would make the analysis have real world application and a more extensive comparison
can be done for the capacity value estimation. For analyzing purposes, another method should be
used instead of Monte Carlo or program with high computing power to yield closer results to
discrete convolution. Since the method is based on randomness, it is more difficult to perform
iterations to determine the percentage load increase to reach the previous LOLE value. Further
research can be performed on different configurations of ESSs. For example, possibly using

39

flywheels to energize pumps to pump water from the lake to the reservoir of a PHS with wind
power could be a configuration. Some PHSs use its own power to energize its pump. With
configurations, there is a need for network communication amongst these systems to work
together.

40

Appendices

41

Appendix A

Matlab Code: Discrete Convolution of Original System
clc
clear all
format long
PL=zeros(24,364); Pl=zeros(60, 365); k=1; Fl=zeros(60, 365); COPAFT_os=
zeros(110, 3);COPAFT_old= zeros(110, 3);COPAFT= zeros(110, 3); NUM=3000;A=
zeros(1, 32);Load_model=zeros(3406,1);
MTTF=[2940 2940 2940 2940 2940 450 450 450 450 1980 1980 1980 1980 1980 1980
1960 1960 1960 1960 1200 1200 1200 960 960 960 960 950 950 950 1150 1100
1100]';
MTTR=[60 60 60 60 60 50 50 50 50 20 20 20 20 20 20 40 40 40 40 50 50 50 40 40
40 40 50 50 50 100 150 150]';
generator=[12 12 12 12 12 20 20 20 20 50 50 50 50 50 50 76 76 76 76 100 100
100 155 155 155 155 197 197 197 350 400 400];
%Reading in LDTS data
data=fopen('LDTS.txt', 'r');
data1= fscanf(data, '%f', [24 365]);
PL=(data1)*2850;
PL=[PL,PL(:,1)]; %Adding the 365th day
%Load Model
for l=1:3000;
for j=1:365;
x=PL(:,j)>= l;
Pl(k,j)=sum(x);
for n=1:23;
y(n,j)= (PL(n,j)<l && PL(n+1,j)>=l);
end
Fl(k,j)=sum(y(:,j));
end
k=k+1;
Load_model(l,1)=l;
Load_model(l,2)= sum(Pl(l,:))/8760;
Load_model(l,3)= sum(Fl(l,:))/8760;
Pl_f(l)=(sum(Pl(l,:))/8760)';
Fl_f(l)=(sum(Fl(l,:))/8760)';
end
%Generation System Model
for i=1:32;
q(i)=(1/MTTF(i))/((1/MTTF(i))+(1/MTTR(i)));
p(i)=(1/MTTR(i))/((1/MTTF(i))+(1/MTTR(i)));
end
%-Intialize blank COPAFT
COPAFT(1,2)=1;

42

COPAFT_old= COPAFT_os + COPAFT;
for i=1:3500
COPAFT_old(i+1,1)=COPAFT_old(i)+1;
if i<=12
COPAFT_old(i+1,2)=q(1,1);
COPAFT_old(i+1,3)=q(1,1)*(1/MTTR(1,1));
end
end
Cn=13;
%Updating COPAFT
COPAFT_os=COPAFT_old(2:length(COPAFT_old),:);
for i=2:32
Cn=Cn+generator(i);
C=COPAFT_os;
for k=1:Cn
Pi(k)=COPAFT_os(k,2);
Fi(k)=COPAFT_os(k,3);
if (k-generator(i))<=0
Pj(k)=1;
Fj(k)=0 ;
else
Pj(k)=COPAFT_os((k-generator(i)),2);
Fj(k)=COPAFT_os((k-generator(i)),3);
end
P(k)=Pi(k)*p(i)+Pj(k)*q(i);
F(k)=(Fi(k)*p(i))+(Fj(k)*q(i))+((Pj(k)-Pi(k))*q(i)*(1/MTTR(i)));
C(k,2)=P(k);
C(k,3)=F(k);
end
COPAFT_os=C;
End
%Ignore data less than 10e-7
for i=1:length(COPAFT_os);
if COPAFT_os(i,2)>(1*10^-7)
COPAFT(i+1,:)=COPAFT_os(i,:);
end
end
%Generation Reserve Model
P_marg(1)=0;
F_marg(1)=0;
j=2;
for i=1:length(COPAFT)
P_r(i)=P(i)-P(i+1);
F_r(i)=F(i)-F(i+1);
R=3405-i+j;
if R>=2851
P_load(i)=0;
F_load(i)=0;
else
P_load(i)=Pl_f(3405-i+j);
F_load(i)=Fl_f(3405-i+j);
end
P_marg(i+1)=P_marg(i)+(P_r(i)*P_load(i));

43

F_marg(i+1)=F_marg(i)+(F_r(i)*P_load(i))+(P_r(i)*F_load(i));
end
P_marg_o= sort(P_marg,'descend');
F_marg_o= sort(F_marg,'descend');
for i=1:length(P_marg)
T(i)=(P_marg_o(i)>(1*10^-7) && F_marg_o(i)>(1*10^-7));
if T==1
Gen_reserve(i,1)=P_marg_o(i);
Gen_reserve(i,2)=F_marg_o(i);
else
end
end
%Results
HLOLE_1=P_marg_o(1)*8760;
LOLF=F_marg_o(1)*8760;

44

Appendix B

Matlab Code: Discrete Convolution of System with Wind Power
clc
clear all
format long
PL=zeros(24,364); Pl=zeros(60, 365); k=1; Fl=zeros(60, 365); COPAFT_os=
zeros(110, 3);COPAFT_old= zeros(110, 3);COPAFT= zeros(110, 3); NUM=3000; A=
zeros(1, 32);Load_model=zeros(3406,1); P=zeros(24,365);
MTTF=[2940 2940 2940 2940 2940 450 450 450 450 1980 1980 1980 1980 1980 1980
1960 1960 1960 1960 1200 1200 1200 960 960 960 960 950 950 950 1150 1100
1100]';
MTTR=[60 60 60 60 60 50 50 50 50 20 20 20 20 20 20 40 40 40 40 50 50 50 40 40
40 40 50 50 50 100 150 150]';
generator=[12 12 12 12 12 20 20 20 20 50 50 50 50 50 50 76 76 76 76 100 100
100 155 155 155 155 197 197 197 350 400 400];
r=0.986; Ar=6362; C=0.5;
%Reading in LDTS data
data=fopen('LDTS.txt', 'r');
data1= fscanf(data, '%f', [24 365]);
data_=fopen('alberta.txt', 'r');
data_2=(fscanf(data_, '%f', [24 365]))*0.27777;
PL=(data1)*2850;
PL=[PL,PL(:,1)]; %Adding the 365th day
%Calculating Wind Power
for i=1:24
for j= 1:365
if data_2(i,j)<11 || data_2(i,j)>25
Pw(i,j)=0;
end
if data_2(i,j)>=11 && data_2(i,j)<=16
Pw(i,j)=.5*r*Ar*C*data_2(i,j)^3;
end
if data_2(i,j)>16 && data_2(i,j)<=25
Pw(i,j)=.5*r*Ar*C*16^3;
end
end
end
%Calculating New Load
PL_= (PL-Pw);
%Load Model
for l=1:3000;
for j=1:365;
x=PL_(:,j)>= l;
Pl(k,j)=sum(x);
for n=1:23;
y(n,j)= (PL_(n,j)<l && PL_(n+1,j)>=l);
end
Fl(k,j)=sum(y(:,j));

45

end
k=k+1;
Load_model(l,1)=l;
Load_model(l,2)= sum(Pl(l,:))/8760;
Load_model(l,3)= sum(Fl(l,:))/8760;
Pl_f(l)=(sum(Pl(l,:))/8760)';
Fl_f(l)=(sum(Fl(l,:))/8760)';
end
%Generation System Model
for i=1:32;
q(i)=(1/MTTF(i))/((1/MTTF(i))+(1/MTTR(i)));
p(i)=(1/MTTR(i))/((1/MTTF(i))+(1/MTTR(i)));
end
%-Intialize blank COPAFT
COPAFT(1,2)=1;
COPAFT_old= COPAFT_os + COPAFT;
for i=1:3500
COPAFT_old(i+1,1)=COPAFT_old(i)+1;
if i<=12
COPAFT_old(i+1,2)=q(1,1);
COPAFT_old(i+1,3)=q(1,1)*(1/MTTR(1,1));
end
end
Cn=13;
%Updating COPAFT
COPAFT_os=COPAFT_old(2:length(COPAFT_old),:);
for i=2:32
Cn=Cn+generator(i);
C=COPAFT_os;
for k=1:Cn
Pi(k)=COPAFT_os(k,2);
Fi(k)=COPAFT_os(k,3);
if (k-generator(i))<=0
Pj(k)=1;
Fj(k)=0 ;
else
Pj(k)=COPAFT_os((k-generator(i)),2);
Fj(k)=COPAFT_os((k-generator(i)),3);
end
P(k)=Pi(k)*p(i)+Pj(k)*q(i);
F(k)=(Fi(k)*p(i))+(Fj(k)*q(i))+((Pj(k)-Pi(k))*q(i)*(1/MTTR(i)));
C(k,2)=P(k);
C(k,3)=F(k);
end
COPAFT_os=C;
end
%Ignore data less than 10e-7
for i=1:length(COPAFT_os);
if COPAFT_os(i,2)>(1*10^-7)
COPAFT(i+1,:)=COPAFT_os(i,:);
end
end

46

%Generation Reserve Model
P_marg(1)=0;
F_marg(1)=0;
j=2;
for i=1:length(COPAFT)
P_r(i)=P(i)-P(i+1);
F_r(i)=F(i)-F(i+1);
R=3405-i+j;
if R>=2851
P_load(i)=0;
F_load(i)=0;
else
P_load(i)=Pl_f(3405-i+j);
F_load(i)=Fl_f(3405-i+j);
end
P_marg(i+1)=P_marg(i)+(P_r(i)*P_load(i));
F_marg(i+1)=F_marg(i)+(F_r(i)*P_load(i))+(P_r(i)*F_load(i));
end
P_marg_o= sort(P_marg,'descend');
F_marg_o= sort(F_marg,'descend');
for i=1:length(P_marg)
T(i)=(P_marg_o(i)>(1*10^-7) && F_marg_o(i)>(1*10^-7));
if T==1
Gen_reserve(i,1)=P_marg_o(i);
Gen_reserve(i,2)=F_marg_o(i);
else
end
end
%Results
HLOLE_1=P_marg_o(1)*8760;
LOLF=F_marg_o(1)*8760;

47

Appendix C

Matlab Code: Iterations of Discrete Convolution to Find LPI
clc
clear all
format long
MTTF=[2940 2940 2940 2940 2940 450 450 450 450 1980 1980 1980 1980 1980 1980
1960 1960 1960 1960 1200 1200 1200 960 960 960 960 950 950 950 1150 1100
1100]';
MTTR=[60 60 60 60 60 50 50 50 50 20 20 20 20 20 20 40 40 40 40 50 50 50 40 40
40 40 50 50 50 100 150 150]';
generator=[12 12 12 12 12 20 20 20 20 50 50 50 50 50 50 76 76 76 76 100 100
100 155 155 155 155 197 197 197 350 400 400];
r=0.986; Ar=6362; C=0.5; tlr=0.001; Fa=1; Q=50; tol=zeros(1,Q);
%Reading in LDTS data
data=fopen('LDTS.txt', 'r');
data1= fscanf(data, '%f', [24 365]);
data_=fopen('alberta.txt', 'r');
data_2=(fscanf(data_, '%f', [24 365]))*0.27777;
PL=(data1)*2850;
PL=[PL,PL(:,1)]; %Adding the 365th day
%Calculating Wind Power
for i=1:24
for j= 1:365
if data_2(i,j)<11 || data_2(i,j)>25
Pw(i,j)=0;
end
if data_2(i,j)>=11 && data_2(i,j)<=16
Pw(i,j)=.5*r*Ar*C*data_2(i,j)^3;
end
if data_2(i,j)>16 && data_2(i,j)<=25
Pw(i,j)=.5*r*Ar*C*60^3;
end
end
end
%Iterations
for U=1:Q
COPAFT_os= zeros(110, 3);COPAFT_old= zeros(110, 3);COPAFT= zeros(110, 3);
Pl=zeros(60, 365);
k=1; Fl=zeros(60, 365); NUM=3000; A= zeros(1,
32);Load_model=zeros(3406,1); P=zeros(24,365);
PL_=(PL-Pw)*Fa;
%Load Model
for l=1:3000;
for j=1:365;
x=PL_(:,j)>= l;
Pl(k,j)=sum(x);
for n=1:23;

48

y(n,j)= (PL_(n,j)<l && PL_(n+1,j)>=l);
end
Fl(k,j)=sum(y(:,j));
end
k=k+1;
Load_model(l,1)=l;
Load_model(l,2)= sum(Pl(l,:))/8760;
Load_model(l,3)= sum(Fl(l,:))/8760;
Pl_f(l)=(sum(Pl(l,:))/8760)';
Fl_f(l)=(sum(Fl(l,:))/8760)';
end
%Generation System Model
for i=1:32;
q(i)=(1/MTTF(i))/((1/MTTF(i))+(1/MTTR(i)));
p(i)=(1/MTTR(i))/((1/MTTF(i))+(1/MTTR(i)));
end
%-Intialize blank COPAFT
COPAFT(1,2)=1;
COPAFT_old= COPAFT_os + COPAFT;
for i=1:3500
COPAFT_old(i+1,1)=COPAFT_old(i)+1;
if i<=12
COPAFT_old(i+1,2)=q(1,1);
COPAFT_old(i+1,3)=q(1,1)*(1/MTTR(1,1));
end
end
Cn=13;
%Updating COPAFT
COPAFT_os=COPAFT_old(2:length(COPAFT_old),:);
for i=2:32
Cn=Cn+generator(i);
C=COPAFT_os;
for k=1:Cn
Pi(k)=COPAFT_os(k,2);
Fi(k)=COPAFT_os(k,3);
if (k-generator(i))<=0
Pj(k)=1;
Fj(k)=0 ;
else
Pj(k)=COPAFT_os((k-generator(i)),2);
Fj(k)=COPAFT_os((k-generator(i)),3);
end
P(k)=Pi(k)*p(i)+Pj(k)*q(i);
F_(k)=(Fi(k)*p(i))+(Fj(k)*q(i))+((Pj(k)-Pi(k))*q(i)*(1/MTTR(i)));
C(k,2)=P(k);
C(k,3)=F_(k);
end
COPAFT_os=C;
end
%Ignore data less than 10e-7
for i=1:length(COPAFT_os);
if COPAFT_os(i,2)>(1*10^-7)

49

COPAFT(i+1,:)=COPAFT_os(i,:);
end
end
%Generation Reserve Model
P_marg(1)=0;
F_marg(1)=0;
j=2;
for i=1:length(COPAFT)
P_r(i)=P(i)-P(i+1);
F_r(i)=F_(i)-F_(i+1);
R=3405-i+j;
if R>=2851
P_load(i)=0;
F_load(i)=0;
else
P_load(i)=Pl_f(3405-i+j);
F_load(i)=Fl_f(3405-i+j);
end
P_marg(i+1)=P_marg(i)+(P_r(i)*P_load(i));
F_marg(i+1)=F_marg(i)+(F_r(i)*P_load(i))+(P_r(i)*F_load(i));
end
P_marg_o= sort(P_marg,'descend');
F_marg_o= sort(F_marg,'descend');
for i=1:length(P_marg)
T(i)=(P_marg_o(i)>(1*10^-7) && F_marg_o(i)>(1*10^-7));
if T==1
Gen_reserve(i,1)=P_marg_o(i);
Gen_reserve(i,2)=F_marg_o(i);
else
end
end
%Results
HLOLE_1=P_marg_o(1)*8760;
LOLF=F_marg_o(1)*8760;
tol(U)=9.302784188074227-HLOLE_1
if tol(U)<0
Fa=Fa-0.0001
elseif tol(U)>1
Fa=Fa+0.01
elseif tol(U)>0.1
Fa=Fa+0.001
elseif tol(U)>0.01
Fa=Fa+0.0001
elseif tol(U)>0.001
Fa=Fa+0.00001
else tol(U)<=tlr
break
end
end

50

Appendix D

Matlab Code: Discrete Convolution of Wind Power and ESS
clc
clear all
format long
MTTF=[2940 2940 2940 2940 2940 450 450 450 450 1980 1980 1980 1980 1980 1980
1960 1960 1960 1960 1200 1200 1200 960 960 960 960 950 950 950 1150 1100
1100]';
MTTR=[60 60 60 60 60 50 50 50 50 20 20 20 20 20 20 40 40 40 40 50 50 50 40 40
40 40 50 50 50 100 150 150]';
generator=[12 12 12 12 12 20 20 20 20 50 50 50 50 50 50 76 76 76 76 100 100
100 155 155 155 155 197 197 197 350 400 400];
r=0.986; Ar=6362; C=0.5; HLOLE_1=zeros(1,10);LOLF=zeros(1,10)
Pw=zeros(24,365); Es=zeros(1,365); CAES=zeros(240,365); n=0; z=0;
CAES_1=zeros(24,365);
CAES_2=zeros(24,365);CAES_3=zeros(24,365);CAES_4=zeros(24,365);CAES_5=zeros(2
4,365);CAES_6=zeros(24,365);
CAES_7=zeros(24,365);
CAES_8=zeros(24,365);CAES_9=zeros(24,365);CAES_10=zeros(24,365);
r=0.986; Ar=6362; C=0.5; y=[24, 48, 72, 96, 120, 144];
size=[100*10^3 150*10^3 200*10^3 400*10^3 1*10^6 2*10^6 3*10^6 4*10^6 5*10^6
10*10^6];
vol=[800*10^3 1.2*10^6 1.6*10^6 3.2*10^6 8*10^6 16*10^6 24*10^6 32*10^6
40*10^6 80*10^6];
EF=[.54 .85 .85 .85 .95 .78];
%Reading in LDTS data
data=fopen('LDTS.txt', 'r');
data1= fscanf(data, '%f', [24 365]);
data_=fopen('alberta.txt', 'r');
data_2=(fscanf(data_, '%f', [24 365]))*0.27777;
PL=(data1)*(2850);
PL=[PL,PL(:,1)]; %Adding the 365th day
%Calculating Wind Power
for i=1:24
for j= 1:365
if data_2(i,j)<11 || data_2(i,j)>25
Pw(i,j)=0;
end
if data_2(i,j)>=11 && data_2(i,j)<=16
Pw(i,j)=.5*r*Ar*C*data_2(i,j)^3;
end
if data_2(i,j)>16 && data_2(i,j)<=25
Pw(i,j)=.5*r*Ar*C*16^3;
end
end
end
%Summing wind power during off peak hours

51

for k= 1:24
for h=1:365
if k<9 || k>16
Es(h)= Pw(k,h)+Es(h);
end
end
end
for w=1:10
for q=1:365
hour(w,q)= floor(Es(q)/size(w));
end
end
%Determining power to discharge during peak hours
for w=1:10
for p=9:16
for q=1:365
if (hour (w,q)==1)
CAES(9+n:16+n-7,q)=size(w);
elseif (hour (w,q)==2)
CAES(9+n:16+n-6,q)=size(w);
elseif (hour (w,q)==3)
CAES(9+n:16+n-5,q)=size(w);
elseif (hour (w,q)==4)
CAES(9+n:16+n-4,q)=size(w);
elseif (hour (w,q)==5)
CAES(9+n:16+n-3,q)=size(w);
elseif (hour (w,q)==6)
CAES(9+n:16+n-2,q)=size(w);
elseif (hour (w,q)==7)
CAES(9+n:16+n-1,q)=size(w);
elseif (hour (w,q)>=8)
CAES(9+n:16+n-0,q)=size(w);
else
end
end
end
n=n+24;
end
CAES_1=CAES(1:24,1:365);
CAES_2=CAES(25:48,1:365);
CAES_3=CAES(49:72,1:365);
CAES_4=CAES(73:96,1:365);
CAES_5=CAES(97:120,1:365);
CAES_6=CAES(121:144,1:365);
CAES_7=CAES(145:168,1:365);
CAES_8=CAES(169:192,1:365);
CAES_9=CAES(193:216,1:365);
CAES_10=CAES(217:240,1:365);
for b= 1:10
%Calculating New Load
k=1; Fl=zeros(60, 365); COPAFT_os= zeros(110, 3);COPAFT_old= zeros(110,
3);COPAFT= zeros(110, 3); NUM=3000; A= zeros(1, 32);Load_model=zeros(3406,1);
P=zeros(24,365);
PL_= (PL-Pw-CAES(1+z:24+z,1:365));

52

%Load Model
for l=1:3000;
for j=1:365;
x=PL_(:,j)>= l;
Pl(k,j)=sum(x);
for e=1:23;
y(e,j)= (PL_(e,j)<l && PL_(e+1,j)>=l);
end
Fl(k,j)=sum(y(:,j));
end
k=k+1;
Load_model(l,1)=l;
Load_model(l,2)= sum(Pl(l,:))/8760;
Load_model(l,3)= sum(Fl(l,:))/8760;
Pl_f(l)=(sum(Pl(l,:))/8760)';
Fl_f(l)=(sum(Fl(l,:))/8760)';
end
%Generation System Model
for i=1:32;
q(i)=(1/MTTF(i))/((1/MTTF(i))+(1/MTTR(i)));
p(i)=(1/MTTR(i))/((1/MTTF(i))+(1/MTTR(i)));
end
%-Intialize blank COPAFT
COPAFT(1,2)=1;
COPAFT_old= COPAFT_os + COPAFT;
for i=1:3500
COPAFT_old(i+1,1)=COPAFT_old(i)+1;
if i<=12
COPAFT_old(i+1,2)=q(1,1);
COPAFT_old(i+1,3)=q(1,1)*(1/MTTR(1,1));
end
end
Cn=13;
%Updating COPAFT
COPAFT_os=COPAFT_old(2:length(COPAFT_old),:);
for i=2:32
Cn=Cn+generator(i);
C=COPAFT_os;
for k=1:Cn
Pi(k)=COPAFT_os(k,2);
Fi(k)=COPAFT_os(k,3);
if (k-generator(i))<=0
Pj(k)=1;
Fj(k)=0 ;
else
Pj(k)=COPAFT_os((k-generator(i)),2);
Fj(k)=COPAFT_os((k-generator(i)),3);
end
P(k)=Pi(k)*p(i)+Pj(k)*q(i);
F(k)=(Fi(k)*p(i))+(Fj(k)*q(i))+((Pj(k)-Pi(k))*q(i)*(1/MTTR(i)));
C(k,2)=P(k);
C(k,3)=F(k);
end
COPAFT_os=C;
end

53

%Ignore data less than 10e-7
for i=1:length(COPAFT_os);
if COPAFT_os(i,2)>(1*10^-7)
COPAFT(i+1,:)=COPAFT_os(i,:);
end
end
%Generation Reserve Model
P_marg(1)=0;
F_marg(1)=0;
j=2;
for i=1:length(COPAFT)
P_r(i)=P(i)-P(i+1);
F_r(i)=F(i)-F(i+1);
R=3405-i+j;
if R>=2851
P_load(i)=0;
F_load(i)=0;
else
P_load(i)=Pl_f(3405-i+j);
F_load(i)=Fl_f(3405-i+j);
end
P_marg(i+1)=P_marg(i)+(P_r(i)*P_load(i));
F_marg(i+1)=F_marg(i)+(F_r(i)*P_load(i))+(P_r(i)*F_load(i));
end
P_marg_o= sort(P_marg,'descend');
F_marg_o= sort(F_marg,'descend');
for i=1:length(P_marg)
T(i)=(P_marg_o(i)>(1*10^-7) && F_marg_o(i)>(1*10^-7));
if T==1
Gen_reserve(i,1)=P_marg_o(i);
Gen_reserve(i,2)=F_marg_o(i);
else
end
end
%Results
HLOLE_1(b)=P_marg_o(1)*8760;
LOLF(b)=F_marg_o(1)*8760;
z=z+24
end

54

Appendix E

Matlab Code: Iterations of Discrete Convolution of Wind Power and ESS to LPI increase
clc
clear all
format long
MTTF=[2940 2940 2940 2940 2940 450 450 450 450 1980 1980 1980 1980 1980 1980
1960 1960 1960 1960 1200 1200 1200 960 960 960 960 950 950 950 1150 1100
1100]';
MTTR=[60 60 60 60 60 50 50 50 50 20 20 20 20 20 20 40 40 40 40 50 50 50 40 40
40 40 50 50 50 100 150 150]';
generator=[12 12 12 12 12 20 20 20 20 50 50 50 50 50 50 76 76 76 76 100 100
100 155 155 155 155 197 197 197 350 400 400];
r=0.986; Ar=6362; C=0.5; HLOLE_1=zeros(1,10);LOLF=zeros(1,10);Fa=ones(1,10);
Q=50; tol=zeros(1,Q);tlr=0.001;
Pw=zeros(24,365); Es=zeros(1,365); CAES=zeros(240,365); n=0;z=0;
CAES_1=zeros(24,365);
CAES_2=zeros(24,365);CAES_3=zeros(24,365);CAES_4=zeros(24,365);CAES_5=zeros(2
4,365);CAES_6=zeros(24,365);
CAES_7=zeros(24,365);CAES_8=zeros(24,365);CAES_9=zeros(24,365);CAES_10=zeros(
24,365);
r=0.986; Ar=6362; C=0.5; y=[24, 48, 72, 96, 120, 144];
size=[100*10^3 150*10^3 200*10^3 400*10^3 1*10^6 2*10^6 3*10^6 4*10^6 5*10^6
10*10^6];
vol=[800*10^3 1.2*10^6 1.6*10^6 3.2*10^6 8*10^6 16*10^6 24*10^6 32*10^6
40*10^6 80*10^6];
EF=[.54 .85 .85 .85 .95 .78];
%Reading in LDTS data
data=fopen('LDTS.txt', 'r');
data1= fscanf(data, '%f', [24 365]);
data_=fopen('alberta.txt', 'r');
data_2=(fscanf(data_, '%f', [24 365]))*0.27777;
PL=(data1)*2850;
PL=[PL,PL(:,1)]; %Adding the 365th day
%Calculating Wind Power
for i=1:24
for j= 1:365
if data_2(i,j)<11 || data_2(i,j)>25
Pw(i,j)=0;
end
if data_2(i,j)>=11 && data_2(i,j)<=16
Pw(i,j)=.5*r*Ar*C*data_2(i,j)^3;
end
if data_2(i,j)>16 && data_2(i,j)<=25
Pw(i,j)=.5*r*Ar*C*16^3;
end
end
end

55

%Summing wind power during off peak hours
for k= 1:24
for h=1:365
if k<9 || k>17
Es(h)= Pw(k,h)+Es(h);
end
end
end
for w=1:10
for q=1:365
hour(w,q)= floor(Es(q)/size(w));
end
end
%Determining power to discharge during peak hours
for w=1:10
for p=9:16
for q=1:365
if (hour (w,q)==1)
CAES(9+n:16+n-7,q)=size(w);
elseif (hour (w,q)==2)
CAES(9+n:16+n-6,q)=size(w);
elseif (hour (w,q)==3)
CAES(9+n:16+n-5,q)=size(w);
elseif (hour (w,q)==4)
CAES(9+n:16+n-4,q)=size(w);
elseif (hour (w,q)==5)
CAES(9+n:16+n-3,q)=size(w);
elseif (hour (w,q)==6)
CAES(9+n:16+n-2,q)=size(w);
elseif (hour (w,q)==7)
CAES(9+n:16+n-1,q)=size(w);
elseif (hour (w,q)>=8)
CAES(9+n:16+n-0,q)=size(w);
else
end
end
end
n=n+24;
end
CAES_1=CAES(1:24,1:365);
CAES_2=CAES(25:48,1:365);
CAES_3=CAES(49:72,1:365);
CAES_4=CAES(73:96,1:365);
CAES_5=CAES(97:120,1:365);
CAES_6=CAES(121:144,1:365);
CAES_7=CAES(145:168,1:365);
CAES_8=CAES(169:192,1:365);
CAES_9=CAES(193:216,1:365);
CAES_10=CAES(217:240,1:365);
for b= 1:10
for U=1:Q
%Calculating New Load

56

k=1; Fl=zeros(60, 365); COPAFT_os= zeros(110, 3);COPAFT_old= zeros(110,
3);COPAFT= zeros(110, 3); NUM=3000; A= zeros(1, 32);Load_model=zeros(3406,1);
P=zeros(24,365);
PL_= (PL-Pw-CAES(1+z:24+z,1:365))*Fa(b);
%Load Model
for l=1:3000;
for j=1:365;
x=PL_(:,j)>= l;
Pl(k,j)=sum(x);
for e=1:23;
y(e,j)= (PL_(e,j)<l && PL_(e+1,j)>=l);
end
Fl(k,j)=sum(y(:,j));
end
k=k+1;
Load_model(l,1)=l;
Load_model(l,2)= sum(Pl(l,:))/8760;
Load_model(l,3)= sum(Fl(l,:))/8760;
Pl_f(l)=(sum(Pl(l,:))/8760)';
Fl_f(l)=(sum(Fl(l,:))/8760)';
end
%Generation System Model
for i=1:32;
q(i)=(1/MTTF(i))/((1/MTTF(i))+(1/MTTR(i)));
p(i)=(1/MTTR(i))/((1/MTTF(i))+(1/MTTR(i)));
end
%-Intialize blank COPAFT
COPAFT(1,2)=1;
COPAFT_old= COPAFT_os + COPAFT;
for i=1:3500
COPAFT_old(i+1,1)=COPAFT_old(i)+1;
if i<=12
COPAFT_old(i+1,2)=q(1,1);
COPAFT_old(i+1,3)=q(1,1)*(1/MTTR(1,1));
end
end
Cn=13;
%Updating COPAFT
COPAFT_os=COPAFT_old(2:length(COPAFT_old),:);
for i=2:32
Cn=Cn+generator(i);
C=COPAFT_os;
for k=1:Cn
Pi(k)=COPAFT_os(k,2);
Fi(k)=COPAFT_os(k,3);
if (k-generator(i))<=0
Pj(k)=1;
Fj(k)=0 ;
else
Pj(k)=COPAFT_os((k-generator(i)),2);
Fj(k)=COPAFT_os((k-generator(i)),3);
end
P(k)=Pi(k)*p(i)+Pj(k)*q(i);
F(k)=(Fi(k)*p(i))+(Fj(k)*q(i))+((Pj(k)-Pi(k))*q(i)*(1/MTTR(i)));
C(k,2)=P(k);
C(k,3)=F(k);

57

end
COPAFT_os=C;
end
%Ignore data less than 10e-7
for i=1:length(COPAFT_os);
if COPAFT_os(i,2)>(1*10^-7)
COPAFT(i+1,:)=COPAFT_os(i,:);
end
end
%Generation Reserve Model
P_marg(1)=0;
F_marg(1)=0;
j=2;
for i=1:length(COPAFT)
P_r(i)=P(i)-P(i+1);
F_r(i)=F(i)-F(i+1);
R=3405-i+j;
if R>=2851
P_load(i)=0;
F_load(i)=0;
else
P_load(i)=Pl_f(3405-i+j);
F_load(i)=Fl_f(3405-i+j);
end
P_marg(i+1)=P_marg(i)+(P_r(i)*P_load(i));
F_marg(i+1)=F_marg(i)+(F_r(i)*P_load(i))+(P_r(i)*F_load(i));
end
P_marg_o= sort(P_marg,'descend');
F_marg_o= sort(F_marg,'descend');
for i=1:length(P_marg)
T(i)=(P_marg_o(i)>(1*10^-7) && F_marg_o(i)>(1*10^-7));
if T==1
Gen_reserve(i,1)=P_marg_o(i);
Gen_reserve(i,2)=F_marg_o(i);
else
end
end
%Results
HLOLE_1(b)=P_marg_o(1)*8760;
LOLF(b)=F_marg_o(1)*8760;
EPNS=sum(P_marg_o)-(0.5*(P_marg_o(1)+P_marg_o(1311)));
EUE=EPNS*8760;

tol(Q)=9.302784188074227-HLOLE_1(b)
if tol(Q)<0
Fa(b)=Fa(b)-0.0001
elseif tol(Q)>1
Fa(b)=Fa(b)+0.01
elseif tol(Q)>0.1
Fa(b)=Fa(b)+0.001
elseif tol(Q)>0.01

58

Fa(b)=Fa(b)+0.0001
elseif tol(Q)>0.001
Fa(b)=Fa(b)+0.00001
elseif tol(Q)<=tlr
break
end
end
z=z+24
end

59

Appendix F

Matlab Code: Monte Carlo of Original System
clc
clear all
format long
PL=zeros(24,364); k=1; S= ones(1,32);
N=10000000;s=1;s_1=1;T0=0.000000001;T_s=0;T_f=0;C_s=0;C_f=0;
MTTF=[2940 2940 2940 2940 2940 450 450 450 450 1980 1980 1980 1980 1980 1980
1960 1960 1960 1960 1200 1200 1200 960 960 960 960 950 950 950 1150 1100
1100]';
MTTR=[60 60 60 60 60 50 50 50 50 20 20 20 20 20 20 40 40 40 40 50 50 50 40 40
40 40 50 50 50 100 150 150]';
generator=[12 12 12 12 12 20 20 20 20 50 50 50 50 50 50 76 76 76 76 100 100
100 155 155 155 155 197 197 197 350 400 400];
%Reading in LDTS data
data=fopen('LDTS.txt', 'r');
data1= fscanf(data, '%f', [24 365]);
PL=(data1)*2850;
PL=[PL,PL(:,1)]; %Adding the 365th day

%Calculate Repair/Failure Rate
R_r = 1./MTTF;
F_r= 1./MTTR;
M=R_r;
%Simulate Mixed time sequential Monte Carlo
for i=1:N
R=rand(1,32);
M1=diag(M);
t=-log(R)/M1;
%Determine the minimum time for the generator to
change state
[tmin index]=min(t); %Locate the specific generator
G=generator*(diag(S));
Gen=sum(G);
%Total generation
X1=ceil(mod(T0,8760));
T1=T0+tmin;
X2=ceil(mod(T1,8760));
for j=X1:X2
Load=PL(j);
if Gen>Load
s=1;
if j==X1
T_s=T_s+ceil(T0)-T0;
elseif j==X2
T_s=T_s+T1-floor(T1);
else
T_s=T_s+1;
end
if s_1==0
C_s=C_s+1;

60

end
s_1=s;
else
s=0;
if j==X1
T_f=T_f+ceil(T0)-T0;
elseif j==X2
T_f=T_f+T1-floor(T1);
else
T_f=T_f+1;
end
if s_1==1
C_f=C_f+1;
end
s_1=s;
end
end
T0=T1;
if S(index)==1
M(index)=F_r(index);
S(index)=0;
else
M(index)=R_r(index);
S(index)=1;
end
end
%Results
LOLP=T_f/T0;
HLOLE=LOLP*8760
HLOLF=(C_s/T0)*8760

61

Appendix G

Matlab Code: Monte Carlo of System with Wind Power
clc
clear all
format long
PL=zeros(24,364); k=1; S= ones(1,32);
N=1000000;s=1;s_1=1;T0=0.000000001;T_s=0;T_f=0;C_s=0;C_f=0;
MTTF=[2940 2940 2940 2940 2940 450 450 450 450 1980 1980 1980 1980 1980 1980
1960 1960 1960 1960 1200 1200 1200 960 960 960 960 950 950 950 1150 1100
1100]';
MTTR=[60 60 60 60 60 50 50 50 50 20 20 20 20 20 20 40 40 40 40 50 50 50 40 40
40 40 50 50 50 100 150 150]';
generator=[12 12 12 12 12 20 20 20 20 50 50 50 50 50 50 76 76 76 76 100 100
100 155 155 155 155 197 197 197 350 400 400];
r=0.986; Ar=6362; C=0.5;
%Reading in LDTS data
data=fopen('LDTS.txt', 'r');
data1= fscanf(data, '%f', [24 365]);
data_=fopen('alberta.txt', 'r');
data_2=(fscanf(data_, '%f', [24 365]))*0.277778;
PL=(data1)*2850;
PL=[PL,PL(:,1)]; %Adding the 365th day
%Calculating Wind Power
for i=1:24
for j= 1:365
if data_2(i,j)<11 || data_2(i,j)>25
P(i,j)=0;
end
if data_2(i,j)>=11 && data_2(i,j)<=16
P(i,j)=.5*r*Ar*C*data_2(i,j)^3;
end
if data_2(i,j)>16 && data_2(i,j)<=25
P(i,j)=.5*r*Ar*C*16^3;
end
end
end
%Combining the system load and wind
PL_=PL-P;
%Calculate Repair/Failure Rate
R_r = 1./MTTF;
F_r= 1./MTTR;
M=R_r;
%Simulate Mixed time sequential Monte Carlo
for i=1:N
R=rand(1,32);
M1=diag(M);

62

t=-log(R)/M1;
%Determine the minimum time for the generator to
change state
[tmin index]=min(t); %Locate the specific generator
G=generator*(diag(S));
Gen=sum(G);
%Total generation
X1=ceil(mod(T0,8760));
T1=T0+tmin;
X2=ceil(mod(T1,8760));
for j=X1:X2
if Gen>PL_(j)
s=1;
if j==X1
T_s=T_s+ceil(T0)-T0;
elseif j==X2
T_s=T_s+T1-floor(T1);
else
T_s=T_s+1;
end
if s_1==0
C_s=C_s+1;
end
s_1=s;
else
s=0;
if j==X1
T_f=T_f+ceil(T0)-T0;
elseif j==X2
T_f=T_f+T1-floor(T1);
else
T_f=T_f+1;
end
if s_1==1
C_f=C_f+1;
end
s_1=s;
end
end
T0=T1;
if S(index)==1
M(index)=F_r(index);
S(index)=0;
else
M(index)=R_r(index);
S(index)=1;
end
end
%Results
LOLP=T_f/T0;
HLOLE=LOLP*8760
LOLF=(C_s/T0)*8760

63

Appendix H

Matlab Code: Iterations of Monte Carlo to Find LPI
clc
clear all
format long
PL=zeros(24,364); k=1; S= ones(1,32);
N=1000000;s=1;s_1=1;T0=0.000000001;T_s=0;T_f=0;C_s=0;C_f=0;
MTTF=[2940 2940 2940 2940 2940 450 450 450 450 1980 1980 1980 1980 1980 1980
1960 1960 1960 1960 1200 1200 1200 960 960 960 960 950 950 950 1150 1100
1100]';
MTTR=[60 60 60 60 60 50 50 50 50 20 20 20 20 20 20 40 40 40 40 50 50 50 40 40
40 40 50 50 50 100 150 150]';
generator=[12 12 12 12 12 20 20 20 20 50 50 50 50 50 50 76 76 76 76 100 100
100 155 155 155 155 197 197 197 350 400 400];
r=0.986; Ar=6362; C=0.5; Q=50; tlr=0.1 ;Fa=1;tol=zeros(1,Q);
%Reading in LDTS data
data=fopen('LDTS.txt', 'r');
data1= fscanf(data, '%f', [24 365]);
data_=fopen('alberta.txt', 'r');
data_2=(fscanf(data_, '%f', [24 365]))*0.277778;
PL=(data1)*2850;
PL=[PL,PL(:,1)]; %Adding the 365th day
%Calculating Wind Power
for i=1:24
for j= 1:365
if data_2(i,j)<11 || data_2(i,j)>25
Pw(i,j)=0;
end
if data_2(i,j)>=11 && data_2(i,j)<=16
Pw(i,j)=.5*r*Ar*C*data_2(i,j)^3;
end
if data_2(i,j)>16 && data_2(i,j)<=25
Pw(i,j)=.5*r*Ar*C*16^3;
end
end
end
%Calculate Repair/Failure Rate
R_r = 1./MTTF;
F_r= 1./MTTR;
M=R_r;
for

U=1:Q
COPAFT_os= zeros(110, 3);COPAFT_old= zeros(110, 3);COPAFT= zeros(110, 3);
Pl=zeros(60, 365);
k=1; Fl=zeros(60, 365); NUM=3000; A= zeros(1,
32);Load_model=zeros(3406,1); P=zeros(24,365);
PL_=(PL-Pw)*Fa;
%Simulate Mixed time sequential Monte Carlo
for i=1:N

64

R=rand(1,32);
M1=diag(M);
t=-log(R)/M1;
%Determine the minimum time for the generator to
change state
[tmin index]=min(t); %Locate the specific generator
G=generator*(diag(S));
Gen=sum(G);
%Total generation
X1=ceil(mod(T0,8760));
T1=T0+tmin;
X2=ceil(mod(T1,8760));
for j=X1:X2
if Gen>PL_(j)
s=1;
if j==X1
T_s=T_s+ceil(T0)-T0;
elseif j==X2
T_s=T_s+T1-floor(T1);
else
T_s=T_s+1;
end
if s_1==0
C_s=C_s+1;
end
s_1=s;
else
s=0;
if j==X1
T_f=T_f+ceil(T0)-T0;
elseif j==X2
T_f=T_f+T1-floor(T1);
else
T_f=T_f+1;
end
if s_1==1
C_f=C_f+1;
end
s_1=s;
end
end
T0=T1;
if S(index)==1
M(index)=F_r(index);
S(index)=0;
else
M(index)=R_r(index);
S(index)=1;
end
end
%Results
LOLP=T_f/T0;
HLOLE=LOLP*8760
LOLF=(C_s/T0)*8760
tol(U)=9.422655345350872-HLOLE

65

if tol(U)<0
Fa=Fa-0.0001
elseif tol(U)>1
Fa=Fa+0.01
elseif tol(U)>0.1
Fa=Fa+0.001
elseif tol(U)>0.01
Fa=Fa+0.0001
elseif tol(U)>0.001
Fa=Fa+0.00001
else tol(U)<=0.1
break
end
end

66

Appendix I

Matlab Code-Monte Carlo of Wind Power and ESS
clc
clear all
format long
PL=zeros(24,364); k=1; S= ones(1,32);
N=1000000;s=1;s_1=1;T0=0.000000001;T_s=0;T_f=0;C_s=0;C_f=0;
MTTF=[2940 2940 2940 2940 2940 450 450 450 450 1980 1980 1980 1980 1980 1980
1960 1960 1960 1960 1200 1200 1200 960 960 960 960 950 950 950 1150 1100
1100]';
MTTR=[60 60 60 60 60 50 50 50 50 20 20 20 20 20 20 40 40 40 40 50 50 50 40 40
40 40 50 50 50 100 150 150]';
generator=[12 12 12 12 12 20 20 20 20 50 50 50 50 50 50 76 76 76 76 100 100
100 155 155 155 155 197 197 197 350 400 400];
r=0.986; Ar=6362; C=0.5; HLOLE=zeros(1,10);LOLF=zeros(1,10);Fa=ones(1,10);
Q=50; tol=zeros(1,Q);tlr=0.5;
Pw=zeros(24,365); Es=zeros(1,365); CAES=zeros(240,365); n=0;z=0;
CAES_1=zeros(24,365);
CAES_2=zeros(24,365);CAES_3=zeros(24,365);CAES_4=zeros(24,365);CAES_5=zeros(2
4,365);CAES_6=zeros(24,365);
CAES_7=zeros(24,365);CAES_8=zeros(24,365);CAES_9=zeros(24,365);CAES_10=zeros(
24,365);
r=0.986; Ar=6362; C=0.5; y=[24, 48, 72, 96, 120, 144];
size=[100*10^3 150*10^3 200*10^3 400*10^3 1*10^6 2*10^6 3*10^6 4*10^6 5*10^6
10*10^6];
vol=[800*10^3 1.2*10^6 1.6*10^6 3.2*10^6 8*10^6 16*10^6 24*10^6 32*10^6
40*10^6 80*10^6];
EF=[.54 .85 .85 .85 .95 .78];
%Reading in LDTS data
data=fopen('LDTS.txt', 'r');
data1= fscanf(data, '%f', [24 365]);
data_=fopen('alberta.txt', 'r');
data_2=(fscanf(data_, '%f', [24 365]))*0.27777;
PL=(data1)*2850;
PL=[PL,PL(:,1)]; %Adding the 365th day
%Calculating Wind Power
for i=1:24
for j= 1:365
if data_2(i,j)<11 || data_2(i,j)>25
P(i,j)=0;
end
if data_2(i,j)>=11 && data_2(i,j)<=16
P(i,j)=.5*r*Ar*C*data_2(i,j)^3;
end
if data_2(i,j)>16 && data_2(i,j)<=25
P(i,j)=.5*r*Ar*C*16^3;
end
end
end

67

%Summing wind power during off peak hours
for k= 1:24
for h=1:365
if k<9 || k>17
Es(h)= P(k,h)+Es(h);
end
end
end
for w=1:10
for q=1:365
hour(w,q)= floor(Es(q)/size(w));
end
end
%Determining power to discharge during peak hours
for w=1:10
for p=9:16
for q=1:365
if (hour (w,q)==1)
CAES(9+n:16+n-7,q)=size(w);
elseif (hour (w,q)==2)
CAES(9+n:16+n-6,q)=size(w);
elseif (hour (w,q)==3)
CAES(9+n:16+n-5,q)=size(w);
elseif (hour (w,q)==4)
CAES(9+n:16+n-4,q)=size(w);
elseif (hour (w,q)==5)
CAES(9+n:16+n-3,q)=size(w);
elseif (hour (w,q)==6)
CAES(9+n:16+n-2,q)=size(w);
elseif (hour (w,q)==7)
CAES(9+n:16+n-1,q)=size(w);
elseif (hour (w,q)>=8)
CAES(9+n:16+n-0,q)=size(w);
else
end
end
end
n=n+24;
end
CAES_1=CAES(1:24,1:365);
CAES_2=CAES(25:48,1:365);
CAES_3=CAES(49:72,1:365);
CAES_4=CAES(73:96,1:365);
CAES_5=CAES(97:120,1:365);
CAES_6=CAES(121:144,1:365);
CAES_7=CAES(145:168,1:365);
CAES_8=CAES(169:192,1:365);
CAES_9=CAES(193:216,1:365);
CAES_10=CAES(217:240,1:365);
for b=1:10
k=1; NUM=3000; A= zeros(1, 32);Load_model=zeros(3406,1); P=zeros(24,365);

68

PL_= (PL-P-CAES(1+z:24+z,1:365))*Fa(b);
%Calculate Repair/Failure Rate
R_r = 1./MTTF;
F_r= 1./MTTR;
M=R_r;
%Simulate Mixed time sequential Monte Carlo
for i=1:N
R=rand(1,32);
M1=diag(M);
t=-log(R)/M1;
%Determine the minimum time for the generator to
change state
[tmin index]=min(t); %Locate the specific generator
G=generator*(diag(S));
Gen=sum(G);
%Total generation
X1=ceil(mod(T0,8760));
T1=T0+tmin;
X2=ceil(mod(T1,8760));
for j=X1:X2
Load=PL_(j);
if Gen>Load
s=1;
if j==X1
T_s=T_s+ceil(T0)-T0;
elseif j==X2
T_s=T_s+T1-floor(T1);
else
T_s=T_s+1;
end
if s_1==0
C_s=C_s+1;
end
s_1=s;
else
s=0;
if j==X1
T_f=T_f+ceil(T0)-T0;
elseif j==X2
T_f=T_f+T1-floor(T1);
else
T_f=T_f+1;
end
if s_1==1
C_f=C_f+1;
end
s_1=s;
end
end
T0=T1;
if S(index)==1
M(index)=F_r(index);
S(index)=0;
else
M(index)=R_r(index);

69

S(index)=1;
end
end
%Results
LOLP=T_f/T0;
HLOLE(b)=LOLP*8760
LOLF(b)=(C_s/T0)*8760
z=z+24
end

70

Appendix J

Matlab Code- Monte Carlo Iterations of Wind Power with ESS to Find LPI
lc
clear all
format long
PL=zeros(24,364); k=1; S= ones(1,32);
N=1000000;s=1;s_1=1;T0=0.000000001;T_s=0;T_f=0;C_s=0;C_f=0;
MTTF=[2940 2940 2940 2940 2940 450 450 450 450 1980 1980 1980 1980 1980 1980
1960 1960 1960 1960 1200 1200 1200 960 960 960 960 950 950 950 1150 1100
1100]';
MTTR=[60 60 60 60 60 50 50 50 50 20 20 20 20 20 20 40 40 40 40 50 50 50 40 40
40 40 50 50 50 100 150 150]';
generator=[12 12 12 12 12 20 20 20 20 50 50 50 50 50 50 76 76 76 76 100 100
100 155 155 155 155 197 197 197 350 400 400];
r=0.986; Ar=6362; C=0.5; HLOLE=zeros(1,10);LOLF=zeros(1,10);Fa=ones(1,10);
Q=50; tol=zeros(1,Q);tlr=0.5;
Pw=zeros(24,365); Es=zeros(1,365); CAES=zeros(240,365); n=0;z=216;
CAES_1=zeros(24,365);
CAES_2=zeros(24,365);CAES_3=zeros(24,365);CAES_4=zeros(24,365);CAES_5=zeros(2
4,365);CAES_6=zeros(24,365);
CAES_7=zeros(24,365);CAES_8=zeros(24,365);CAES_9=zeros(24,365);CAES_10=zeros(
24,365);
r=0.986; Ar=6362; C=0.5; y=[24, 48, 72, 96, 120, 144];
size=[100*10^3 150*10^3 200*10^3 400*10^3 1*10^6 2*10^6 3*10^6 4*10^6 5*10^6
10*10^6];
vol=[800*10^3 1.2*10^6 1.6*10^6 3.2*10^6 8*10^6 16*10^6 24*10^6 32*10^6
40*10^6 80*10^6];
EF=[.54 .85 .85 .85 .95 .78];
%Reading in LDTS data
data=fopen('LDTS.txt', 'r');
data1= fscanf(data, '%f', [24 365]);
data_=fopen('alberta.txt', 'r');
data_2=(fscanf(data_, '%f', [24 365]))*0.27777;
PL=(data1)*2850;
PL=[PL,PL(:,1)]; %Adding the 365th day
%Calculating Wind Power
for i=1:24
for j= 1:365
if data_2(i,j)<11 || data_2(i,j)>25
P(i,j)=0;
end
if data_2(i,j)>=11 && data_2(i,j)<=16
P(i,j)=.5*r*Ar*C*data_2(i,j)^3;
end
if data_2(i,j)>16 && data_2(i,j)<=25
P(i,j)=.5*r*Ar*C*16^3;
end
end
end

71

%Summing wind power during off peak hours
for k= 1:24
for h=1:365
if k<9 || k>17
Es(h)= P(k,h)+Es(h);
end
end
end
for w=1:10
for q=1:365
hour(w,q)= floor(Es(q)/size(w));
end
end
%Determining power to discharge during peak hours
for w=1:10
for p=9:16
for q=1:365
if (hour (w,q)==1)
CAES(9+n:16+n-7,q)=size(w);
elseif (hour (w,q)==2)
CAES(9+n:16+n-6,q)=size(w);
elseif (hour (w,q)==3)
CAES(9+n:16+n-5,q)=size(w);
elseif (hour (w,q)==4)
CAES(9+n:16+n-4,q)=size(w);
elseif (hour (w,q)==5)
CAES(9+n:16+n-3,q)=size(w);
elseif (hour (w,q)==6)
CAES(9+n:16+n-2,q)=size(w);
elseif (hour (w,q)==7)
CAES(9+n:16+n-1,q)=size(w);
elseif (hour (w,q)>=8)
CAES(9+n:16+n-0,q)=size(w);
else
end
end
end
n=n+24;
end
CAES_1=CAES(1:24,1:365);
CAES_2=CAES(25:48,1:365);
CAES_3=CAES(49:72,1:365);
CAES_4=CAES(73:96,1:365);
CAES_5=CAES(97:120,1:365);
CAES_6=CAES(121:144,1:365);
CAES_7=CAES(145:168,1:365);
CAES_8=CAES(169:192,1:365);
CAES_9=CAES(193:216,1:365);
CAES_10=CAES(217:240,1:365);
for b=10:10
for U=1:Q
k=1; NUM=3000; A= zeros(1, 32);Load_model=zeros(3406,1); P=zeros(24,365);

72

PL_= (PL-P-CAES(1+z:24+z,1:365))*Fa(b);
%Calculate Repair/Failure Rate
R_r = 1./MTTF;
F_r= 1./MTTR;
M=R_r;
%Simulate Mixed time sequential Monte Carlo
for i=1:N
R=rand(1,32);
M1=diag(M);
t=-log(R)/M1;
%Determine the minimum time for the generator to
change state
[tmin index]=min(t); %Locate the specific generator
G=generator*(diag(S));
Gen=sum(G);
%Total generation
X1=ceil(mod(T0,8760));
T1=T0+tmin;
X2=ceil(mod(T1,8760));
for j=X1:X2
Load=PL_(j);
if Gen>Load
s=1;
if j==X1
T_s=T_s+ceil(T0)-T0;
elseif j==X2
T_s=T_s+T1-floor(T1);
else
T_s=T_s+1;
end
if s_1==0
C_s=C_s+1;
end
s_1=s;
else
s=0;
if j==X1
T_f=T_f+ceil(T0)-T0;
elseif j==X2
T_f=T_f+T1-floor(T1);
else
T_f=T_f+1;
end
if s_1==1
C_f=C_f+1;
end
s_1=s;
end
end
T0=T1;
if S(index)==1
M(index)=F_r(index);
S(index)=0;
else
M(index)=R_r(index);

73

S(index)=1;
end
end
%Results
LOLP=T_f/T0;
HLOLE(b)=LOLP*8760
LOLF(b)=(C_s/T0)*8760
tol(Q)=9.422655345350872-HLOLE(b)
if tol(Q)<0
Fa(b)=Fa(b)-0.0001
elseif tol(Q)>1
Fa(b)=Fa(b)+0.01
elseif tol(Q)>0.5
Fa(b)=Fa(b)+0.001
elseif tol(Q)<=tlr
break
end
end
z=z+24
end

74

Appendix K

Test System Data
Table A.1: Generating Unit Reliability Data
Unit Size MW
12
20
50
76
100
155
197
350
400

Number of
Unit
5
4
6
4
3
4
3
1
2

Forced outage
Rate
0.02
0.10
0.01
0.02
0.04
0.04
0.05
0.08
0.12

MTTF

MTTR

2940
450
1980
1960
1200
960
950
1150
1100

60
50
20
40
50
40
50
100
150

Figure A.1: IEEE RTS Load Data

1

0.9

0.8

0.7

0.6

0.5

0.4

0

1000

2000

3000

4000

75

5000

6000

7000

8000

9000

REFERENCES

76

REFERENCES
1. Keane, A., Milligan, M., Dent, C. J., Hasche, B., D’Annunzio, C., Dragoon, K..,
O’Malley, M. (2011). Capacity Value of Wind Power. IEEE Transactions on Power
Systems, 26(2), 564–572. doi:10.1109/TPWRS.2010.2062543
2. Lu, M.-S., Chang, C.-L., Lee, W.-J., & Wang, L. (2008). Combining the Wind Power
Generation System with Energy Storage Equipments. In IEEE Industry Applications
Society Annual Meeting, 2008. IAS ’08 (pp. 1–6). doi:10.1109/08IAS.2008.139
3. Tuohy, A., & O’Malley, M. (2009). Impact of pumped storage on power systems with
increasing wind penetration. In IEEE Power Energy Society General Meeting, 2009. PES
’09 (pp. 1–8). doi:10.1109/PES.2009.5275839
4. Garver, L. L. (1966). Effective Load Carrying Capability of Generating Units. IEEE
Transactions on Power Apparatus and Systems, PAS-85(8), 910–919.
doi:10.1109/TPAS.1966.291652
5. Kirby, B., Milligan, M., Lovekin, J., Jackson, K. J., Makarov, Y., & Hawkins, D. (2004,
July). California Renewables Portfolio Standard Renewable Generation Integration Cost
Analysis. California Energy Commission.
6. Mitra, J. (n.d.). Power System Reliability Using Discrete Convolution (pp. 151–169).
7. Abdullah, M. ., Muttaqi, K. M., Agalgaonkar, A., & Sutanto, D. (2013). Estimating the
capacity value of energy storage integrated with wind power generation. In 2013 IEEE
Power and Energy Society General Meeting (PES) (pp. 1–5).
doi:10.1109/PESMG.2013.6672996
8. Milligan, M., & Porter, K. (2005). Determining the Capacity Value of Wind: A Survey of
Methods and Implementation. National Renewable Energy Laboratory.
9. Taylor, J., & Halnes, A. (2010). Analysis Of compressed air energy storage. In PCIC
Europe 2010 Conference Record (pp. 1–5).

77

10. Boyes, J. D., & Clark, N. H. (2000). Technologies for energy storage. Flywheels and
super conducting magnetic energy storage. In IEEE Power Engineering Society Summer
Meeting, 2000 (Vol. 3, pp. 1548–1550 vol. 3). doi:10.1109/PESS.2000.868760
11. Hassenzahl, W. V. (1983). Superconducting magnetic energy storage. Proceedings of the
IEEE, 71(9), 1089–1098. doi:10.1109/PROC.1983.12727
12. Sparacino, A., Reed, G. F., Kerestes, R. J., Grainger, B. M., & Smith, Z. T. (2012).
Survey of battery energy storage systems and modeling techniques. In 2012 IEEE Power
and Energy Society General Meeting (pp. 1–8). doi:10.1109/PESGM.2012.6345071
13. Vazquez, S., Lukic, S. M., Galvan, E., Franquelo, L. G., & Carrasco, J. M. (2010).
Energy Storage Systems for Transport and Grid Applications. IEEE Transactions on
Industrial Electronics, 57(12), 3881–3895. doi:10.1109/TIE.2010.2076414
14. Jiang, X., Zhang, J., & Jian, W. (2013). The Analysis of Ultracapacitor Charging
Efficiency. In 2013 Fifth International Conference on Computational and Information
Sciences (ICCIS) (pp. 1198–1201). doi:10.1109/ICCIS.2013.317
15. Bradshaw, D. T. (2000). Pumped hydroelectric storage (PHS) and compressed air energy
storage (CAES). In IEEE Power Engineering Society Summer Meeting, 2000 (Vol. 3, pp.
1551–1573). doi:10.1109/PESS.2000.868761
16. Nomura, S., Shintomi, T., Akita, S., Nitta, T., Shimada, R., & Meguro, S. (2010).
Technical and Cost Evaluation on SMES for Electric Power Compensation. IEEE
Transactions on Applied Superconductivity, 20(3), 1373–1378.
doi:10.1109/TASC.2009.2039745
17. Kaldellis, J. K., Zafirakis, D., & Kavadias, K. (2007). Techno-economic comparison of
energy storage systems for island autonomous electric networks. Renewables &
Sustainable Energy Reviews.
18. Simpson, A., & Walker, G. R. (2002). Lifecycle costs of ultracapacitors in electric
vehicle applications. In Power Electronics Specialists Conference, 2002. pesc 02. 2002
IEEE 33rd Annual (Vol. 2, pp. 1015–1020 vol.2). doi:10.1109/PSEC.2002.1022588
19. Mitra, J., & Singh, C. (2005). Monte Carlo Simulation and Intelligent Search Methods
(pp. 23–38). IEEE.

78

20. Subcommittee, P. M. (1979). IEEE Reliability Test System. IEEE Transactions on Power
Apparatus and Systems, PAS-98(6), 2047–2054. doi:10.1109/TPAS.1979.319398
21. Xian, W., Yuan, W., Yan, Y., & T. A., C. (2009). Minimize Frequency Fluctuations of
Isolated Power System with Wind Farm by Using Superconducting Magnetic Energy
Storage.
22. Zhou, F., Abbey, C., Jiao, L., & Ooi Teck, B. (n.d.). Use of Large Capacity SMES to
Improve the Power Quality and Stability of Wind Farms.
23. Tripathy, S. C., Kalantar, M., & Balasubramanian, R. (1991). Dynamics and Stability of
Wind and Diesel Turbine Generators with Superconducting Magnetic Energy Storage
Unit on an Isolated System. Presented at the IEEE Transactions on Energy Conversion.
24. Takahashi, R., & Tamura, J. (2008). Frequency Control of Isolated Power System with
Wind Farm by Using Flywheel Energy Storage System. Presented at the 2008
International Conference on Electrical Machines.
25. Renewable Energy Policy Network for the 21st Century. (2013). Renewables 2013
Global Status Report.
26. U.S. Department of Energy. (2008). Energy Efficiency and Renewable Energy.
27. "Primary Energy Consumption by Source." U.S. Energy Information Administration
(EIA). 1 Sept. 2014. Web. 13 Oct. 2014. <http://www.eia.gov>.

79