UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

The value of stochastic optimization in reducing the cost of day-ahead wind forecast uncertainty Tadesse, Mehretab Getachew 2018

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata

Download

Media
24-ubc_2018_september_tadesse_mehretab.pdf [ 6MB ]
Metadata
JSON: 24-1.0369051.json
JSON-LD: 24-1.0369051-ld.json
RDF/XML (Pretty): 24-1.0369051-rdf.xml
RDF/JSON: 24-1.0369051-rdf.json
Turtle: 24-1.0369051-turtle.txt
N-Triples: 24-1.0369051-rdf-ntriples.txt
Original Record: 24-1.0369051-source.json
Full Text
24-1.0369051-fulltext.txt
Citation
24-1.0369051.ris

Full Text

THE VALUE OF STOCHASTIC OPTIMIZATION IN REDUCING THE COST OF DAY-AHEAD WIND FORECAST UNCERTAINTY  by Mehretab Getachew Tadesse B.Sc., Addis Ababa University, 2013  A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF  MASTER OF APPLIED SCIENCE in THE FACULTY OF GRADUATE AND POSTDOCTORAL STUDIES (Civil Engineering) THE UNIVERSITY OF BRITISH COLUMBIA (Vancouver) July 2018     © Mehretab Getachew Tadesse, 2018  ii   The following individuals certify that they have read, and recommend to the Faculty of Graduate and Postdoctoral Studies for acceptance, a thesis/dissertation entitled: The Value of Stochastic Optimization in Reducing the Cost of Day-Ahead Wind Forecast Uncertainty  submitted by Mehretab Tadesse in partial fulfillment of the requirements for the degree of Master of Applied Science in Civil Engineering  Examining Committee: Dr. Ziad Shawwash Supervisor  Dr. Steven Weijs Supervisory Committee Member    iii  Abstract Wind power is gaining popularity around the globe mainly because of its advantages such as its renewability and its low environmental impact. However, wind power has some operational disadvantages. Wind power is uncertain, variable, and non-dispatchable. These three properties add cost on wind power producers and electric utilities. In this thesis, stochastic optimization (SO) is applied to minimize the cost of day-ahead (DA) wind forecast uncertainty. SO requires the statistical properties of an uncertain parameter to be accurately represented. In this thesis, scenarios were used to represent the wind forecast error (WFE). To generate WFE scenarios that properly characterize the wind forecast uncertainty, first the statistical properties of the DA WFE was investigated. Then two competing scenario generation algorithms, one based on Markov chain (MC) and the other based on autoregressive moving average (ARMA), were proposed and tested. Both models were successful in generating representative WFE scenarios, but the ARMA-based model was found to be better at recreating the statistical properties of the WFE. These two algorithms required the generation of numerous scenario to confidently capture the WFE probability distribution. Since using all the generated scenarios in a SO problem was infeasible, a scenario reduction algorithm based on probability distance was applied to reduce the number of scenarios while preserving the information they contain. The scenarios obtained after the reduction process were used as inputs into a two-state linear SO model. The results of the optimization showed that, under normal operating conditions, SO can reduce the cost of DA wind forecast uncertainty by up to 70%. The SO models were also found to be better at avoiding load shedding and wind power curtailment. When comparing the two scenario generation algorithms with respect to cost savings, the MC-based model resulted in higher cost savings because the scenarios it generated were better at capturing extreme WFEs, which led to less load shedding events. Sensitivity analysis conducted by varying input assumptions showed that the savings from SO vary considerably based on the modeling assumptions and that care should be taken when designing a cost savings evaluation strategy.   iv  Lay Summary Wind power is gaining popularity around the world because of its advantageous traits like its renewability and low pollutant emission. However, wind energy has some disadvantages. One of the main disadvantages is that it is uncertain, which means it is difficult to forecast. In this thesis, stochastic optimization (SO), a scheduling approach that explicitly considers the uncertainty of wind forecast, is investigated as an approach to minimize the effect of wind power uncertainty in the BC Hydro system. To achieve this, the property of the wind forecast uncertainty in BC was first studied. Then, methods of producing representative possible outcome scenarios of the wind generation for a given forecast were devised. SO was then carried out considering these possible outcomes of wind generation instead of the single wind generation forecast. The results showed that SO can significantly decrease the effect of wind uncertainty in the BC Hydro system.    v  Preface The thesis is an original work of the author, Mehretab Tadesse. The author developed and applied both MC-based and ARMA-based scenario generation algorithms presented in Chapter 3. The stochastic linear programming model discussed in Chapters 3 is based on a propriety BC Hydro model, but the author made the necessary modifications to convert it into a two-stage stochastic linear optimization model. Dr. Ziad Shawwash, the supervisor of the thesis, has provided several valuable suggestions to improve the stochastic model.  The author also conducted the writing of all chapters of the thesis. Mr. Jonathan van Groll has proofread the thesis and has given helpful comments. Dr. Ziad Shawwash has also provided constructive edits and comments that have been incorporated in the final version of the thesis.  vi  Table of Contents  Abstract ................................................................................................................................... iii Lay Summary .......................................................................................................................... iv Preface ...................................................................................................................................... v Table of Contents ................................................................................................................... vi List of Tables ........................................................................................................................... xi List of Figures ......................................................................................................................... xii List of Abbreviations ............................................................................................................. xvi Acknowledgements ..............................................................................................................xix  Introduction .......................................................................................................... 1 1.1 Background and Motivation ................................................................................................................ 3 1.2 Goals and Objectives............................................................................................................................ 7 1.3 Organization of the Thesis .................................................................................................................. 7  Literature Review ................................................................................................. 9 2.1 Wind in BC ............................................................................................................................................. 9 2.2 Wind Integration Cost ....................................................................................................................... 10 2.2.1 The BC Hydro Wind Integration Study .................................................................................... 11 2.2.2 Review of Other Wind Integration Studies.............................................................................. 14 vii  2.3 Day-Ahead Wind Forecast Error .................................................................................................... 17 2.4 Scenario Generation ........................................................................................................................... 21 2.4.1 Review of the Stochastic Wind Data and Scenario Generation Techniques ................... 23 2.4.2 Scenario Reduction ........................................................................................................................ 30 2.5 Stochastic Optimization of Wind Power ....................................................................................... 34  Modeling Methodology .......................................................................................40 3.1 Statistical Analysis of the Day Ahead Wind Forecast Error ..................................................... 40 3.1.1 Statistical Moments ........................................................................................................................ 40 3.1.2 Stationarity and Seasonality .......................................................................................................... 41 3.1.3 Probability Distributions ............................................................................................................... 42 3.1.4 Autocorrelation and Partial Autocorrelation Functions ....................................................... 43 3.2 A Markov Chain Based Scenario Generation Model .................................................................. 44 3.3 An ARMA Based Scenario Generation Model.............................................................................. 47 3.4 Evaluation of Scenario Quality ......................................................................................................... 50 3.5 Scenario Reduction Approach ......................................................................................................... 51 3.6 Optimization Models .......................................................................................................................... 55 3.6.1 Hydraulic Simulation Model (HYSIM) ........................................................................................ 55 3.6.2 The Generalized Optimization Model (GOM) ........................................................................ 56 3.7 Modifications Made to the Generalized Optimization Model (GOM) .................................... 60 viii  3.7.1 Modification to Optimization Algorithm ................................................................................... 61 3.7.2 Added Variables and Parameters ................................................................................................ 63 3.7.3 Modification to the Objective Function .................................................................................... 65 3.7.4 Modifications to the Constraints ................................................................................................ 67 3.7.5 Modification to Operating Reserves .......................................................................................... 69 3.8 Estimation of Stochastic Savings ...................................................................................................... 70  Case Study and Inputs ........................................................................................73 4.1 The BC Hydro System ....................................................................................................................... 73 4.2 Model Inputs and Assumptions ........................................................................................................ 74 4.2.1 Load ................................................................................................................................................... 74 4.2.2 Reservoir Inflows ............................................................................................................................ 75 4.2.3 Price ................................................................................................................................................... 77 4.2.4 Wind Data ........................................................................................................................................ 78 4.2.5 Forebay and Storage Targets ....................................................................................................... 80 4.2.6 Optimized Plants ............................................................................................................................. 80 4.2.7 Reserves............................................................................................................................................ 81 4.3 Other Inputs and Assumptions ........................................................................................................ 81 4.4 Cases Investigated ............................................................................................................................... 82  Results and Discussion .......................................................................................84 ix  5.1 Day-Ahead Wind Forecast Error Properties ............................................................................... 84 5.1.1 Statistical Moments ........................................................................................................................ 85 5.1.2 Stationarity and Seasonality .......................................................................................................... 85 5.1.3 Probability Distribution ................................................................................................................. 87 5.1.4 Autocorrelation and Partial Autocorrelation Functions ....................................................... 90 5.2 Scenario Generation Using Markov Chains .................................................................................. 91 5.3 Scenario Generation Using ARMA.................................................................................................. 98 5.4 Scenario Quality Evaluation ............................................................................................................ 103 5.4.1 Statistical Moments ...................................................................................................................... 105 5.4.2 Marginal and Conditional Distributions................................................................................... 105 5.4.3 ACF and PACF .............................................................................................................................. 108 5.4.4 Stability ............................................................................................................................................ 108 5.5 Scenario Reduction ........................................................................................................................... 110 5.6 Optimization Results ........................................................................................................................ 114 5.6.1 Assessment Using the Expected Optimal Values .................................................................. 114 5.6.2 Out-of-Sample Assessment ........................................................................................................ 116 5.6.3 Load Shedding ............................................................................................................................... 121 5.6.4 Wind Power Curtailment ........................................................................................................... 122 5.7 Sensitivity Analysis ............................................................................................................................ 123 x  5.7.1 Sensitivity to Hydrologic Conditions ....................................................................................... 123 5.7.2 Sensitivity to Load Shedding Penalty Cost .............................................................................. 125 5.7.3 Sensitivity to AS Market Depth Assumption ......................................................................... 126  Conclusion and Recommendations for Future Work ................................. 128 6.1 Statistical Properties of the WFE .................................................................................................. 128 6.2 Scenario Generation ......................................................................................................................... 129 6.3 The Value of Stochastic Optimization in Reducing the DA Wind Integration Cost ......... 130 6.4 The Effect of the Scenarios Generation Method in the Stochastic Savings ......................... 131 6.5 Future Work ...................................................................................................................................... 131 Bibliography ......................................................................................................................... 133  xi  List of Tables Table 1 -  Wind Farms Operating in BC ........................................................................................................... 10 Table 2 - Summary of wind integration cost for the BC Hydro system ................................................... 14 Table 3 - Some of BC Hydro's generation resources modeled .................................................................. 80 Table 4 – Other inputs and assumptions that went into GOM .................................................................. 82 Table 5 - Water, wind and load year combination ......................................................................................... 83 Table 6 – WFE statistics ....................................................................................................................................... 85 Table 7 - Pj,i transition probability matrix – ending WFE state given WFE state at t-1 ......................... 95 Table 8 - Pk,i transition probability matrix– ending WFE state given WFE state at t-2 ......................... 95 Table 9 - Pf,i transition probability matrix– ending WFE state given Forecast state at t....................... 95 Table 10 - Cumulative transition probability matrix Cj,i ............................................................................... 96 Table 11 - AIC and -log likelihood comparisons ............................................................................................. 99 Table 12 - Statistical moments comparison ................................................................................................... 105 Table 13 - Summary of sensitivity analysis results ........................................................................................ 127  xii  List of Figures Figure 1 - Global installed wind capacity ............................................................................................................ 1 Figure 2 - Average LCOE in good to excellent wind sites ............................................................................. 3 Figure 3 - Effects of wind power in the power system ................................................................................. 12 Figure 4 – A typical wind power curve with commonly used terminology .............................................. 19 Figure 5 - Scenario tree (on the left) and scenario fan (on the right) ....................................................... 22 Figure 6 - Probability transformation................................................................................................................. 49 Figure 7 - A simplified flowchart showing the process of scenario generation ....................................... 54 Figure 8 - Renewable generation resources in British Columbia ................................................................ 74 Figure 9 - BC Hydro load forecast for 2035 .................................................................................................... 75 Figure 10 - Cumulative inflows for BC Hydro's large reservoirs (1964-1973) ........................................ 76 Figure 11 - Average inflows into the Williston Reservoir (GMS) (1930-2009) ....................................... 76 Figure 12 - Price data for the Mid-C trading hub (2035/36) ........................................................................ 77 Figure 13 - AS prices (2016/2017) ...................................................................................................................... 78 Figure 14 - Wind reserves ................................................................................................................................... 81 Figure 15 - Cases Investigated ............................................................................................................................. 83 Figure 16 – The DA wind forecast error time series .................................................................................... 84 Figure 17 – Hourly MAE ....................................................................................................................................... 86 Figure 18 – Hourly std of the WFE.................................................................................................................... 86 xiii  Figure 19 - Monthly average wind generation ................................................................................................. 86 Figure 20 - Monthly MAE ..................................................................................................................................... 86 Figure 21 - WFE histogram .................................................................................................................................. 88 Figure 22 – Scatter plot of the wind forecast and WFE ............................................................................... 88 Figure 23 - MAE for each forecast bin .............................................................................................................. 88 Figure 24 - Std for each forecast bin ................................................................................................................. 88 Figure 25 – WFE marginal histogram and best fit distributions .................................................................. 89 Figure 26 – WFE conditional histogram and best fit distributions for the second forecast bin [0.05,0.1) .................................................................................................................................................................................... 89 Figure 27 - WFE histogram and best fit distributions for the 10th forecast bin [0.45,0.5) ................... 90 Figure 28 - WFE histogram and best fit distributions for the 19th forecast bin [0.9,0.95) ................... 90 Figure 29 - Sample ACF of the WFE time series ............................................................................................ 91 Figure 30 - Sample PACF of the WFE time series .......................................................................................... 91 Figure 31 - ACF comparison ............................................................................................................................... 92 Figure 32 -PACF comparison .............................................................................................................................. 92 Figure 33 - Mean (bottom half) and std (upper half) computed after considering n number of scenarios for a typical day (MC-based model). The various lines represent the 24 hours in a day. ..................... 98 Figure 34 – The marginal ECDF used for transforming the WFE to standard normal and the various conditional ECDFs used for reverse transforming normally distributed ARMA output back to WFE .................................................................................................................................................................................. 101 xiv  Figure 35 - Number of scenarios vs statistical parameters: standard deviation (upper half) and mean (lower half) for the 24 hours of a typical medium wind forecast day (ARMA model) ........................ 103 Figure 36 - Scenarios generated using the MC-based model. .................................................................... 104 Figure 37 – 50 Scenarios generated using the ARMA-based model. ........................................................ 104 Figure 38 – Scatter plot of WFE obtained from the MC-based model ................................................... 107 Figure 39 - Scatter plot of WFE obtained from the ARMA-based model .............................................. 107 Figure 40 - Comparison between the marginal ECDFs of the observed WFE time series, MC-based model simulated time series, and ARMA-based model simulated time series. ...................................... 107 Figure 41 - Autocorrelation function comparison ........................................................................................ 109 Figure 42 - Partial autocorrelation function comparison ............................................................................ 109 Figure 43 - Scenario stability analysis – means (lower half) and standard deviations (upper half) of 25 simulations from the MC-based model for a typical day ............................................................................ 110 Figure 44 - Scenario stability analysis – means (lower half) and standard deviations (upper half) of 25 simulations from the ARMA-based model for a typical day ....................................................................... 110 Figure 45- Relationship between number of scenarios and percentage of information captured .... 111 Figure 46 - 40 reduced scenarios from 1000 (MC). The line thickness indicates relative probability .................................................................................................................................................................................. 113 Figure 47 - 40 reduced scenarios from 1000 (ARMA). The line thickness indicates relative probability .................................................................................................................................................................................. 113 Figure 48 - Statistical parameter comparison between the full and reduced scenario sets for a typical day ............................................................................................................................................................................ 114 xv  Figure 49 – Difference in expected optimal values (the optimal value of model (b) is the reference) .................................................................................................................................................................................. 115 Figure 50 – Relative cost comparison based on the expected optimal values....................................... 116 Figure 51 - Difference in the daily actual optimal values ............................................................................. 117 Figure 52 – Relative cost comparison based on the actual optimal values ............................................. 118 Figure 53 - Comparison of the objective function components of models (c) and (d) with model (b). .................................................................................................................................................................................. 118 Figure 54 - Forecasted, observed and the 40 reduced wind generation scenarios for the 110th day (the MC-based model) ................................................................................................................................................. 120 Figure 55 - Forecasted, observed and the 40 reduced wind generation scenarios for the 110th day (ARMA-based model) .......................................................................................................................................... 120 Figure 56 - Cumulative load shedding ............................................................................................................. 121 Figure 57 - Cumulative wind power curtailment .......................................................................................... 122 Figure 58 – Cost comparison for different water years ............................................................................. 124 Figure 59 – Cost comparison for the low load shedding penalty cost assumption .............................. 125 Figure 60 – Cost comparison for low AS market depth assumption ...................................................... 127  xvi  List of Abbreviations ACC  Autocorrelation Coefficient ACF  Autocorrelation Function AGC  Automatic Generation Control AIC  Akaike Information Criterion ANN  Artificial Neural Network AR  Auto-Regressive ARD  Arrow Lake generating plant ARIMA  Auto-Regressive Integrated Moving Average  ARMA  Auto-Regressive Moving Average AS  Ancillary Service BC Hydro British Columbia Hydro and Power Authority BC  British Columbia BIC  Bayesian Information Criterion CAISO  California Independent System Operator CAPEX Capacity Expansion CDF  Cumulative Density Function CSTPM Cumulative State Transition Probability Matrix DA  Day-Ahead DOE  Department of Energy ECDF  Empirical Cumulative Density Function ERCOT Electric Reliability Council of Texas GMS  GM Shrum generating plant GOM  Generalized Optimization Model GW  Gigawatt GWh   Gigawatt Hour HYSIM  Hydraulic Simulation IPP  Independent Power Producer IRP  Integrated Resource Plan ISO  Independent System Operator xvii  K  Kurtosis LARIMA Limited Autoregressive Integrated Moving Average LCOE  Levelized Cost of Energy LFE  Load Forecast Error LLN  Law of Large Numbers LP  Linear Programming MA  Moving Average MAE  Mean Absolute Wind Forecast Error MARIMA Multivariate Autoregressive Integrated Moving Average MC  Markov Chain MCA  Mica plant Mid-C  Middle Columbia MVW  Marginal Value of Water MW   Megawatt MWh  Megawatt Hour NWP  Numerical Weather Prediction ORO  Operating Reserve Obligation PACC  Partial Autocorrelation Coefficient PACF  Partial Autocorrelation Function PCN  Peace Canyon generating plant PDF  Probability Density Function PSCo   Public Service Company of Colorado REV  Revelstoke generating plant ROR  Run-of-River RPS  Renewable Portfolio Standards RT  Real-Time SARIMA Seasonal Autoregressive Integrated Moving Average SCORE Statistical Correction to Output from a Record Extension SO  Stochastic Optimization SPP  Southern Power Pool STC  Site-C generating plant xviii  Std  Standard Deviation STOM  Short Term Optimization Model UC  Unit Commitment UCC  Unit Capacity Cost UEC  Unit Energy Cost WFE  Wind Forecast Error WIS  Wind Integration Study WPFE  Wind Power Forecast Error WPTS  Wind Power Time Series WRF  Weather Research and Forecasting WSFE  Wind Speed Forecast Error WSTS  Wind Speed Time Series  xix  Acknowledgements First of all, I would like to thank my supervisor Dr. Ziad Shawwash for providing continued support, advice, and encouragement throughout my thesis. I would also like to thank NSERC and BC Hydro for providing the funding required for me to pursue my post-graduate education. I offer special gratitude to the other UBC students at BC Hydro, Abhishek Agrawal, and Jonathan van Groll, for helping to come up with ideas, proofreading the thesis and offering support. I would also like to thank Postdoctoral Fellow Dr. Quentin Desreumaux, whose suggestions have been very valuable for the thesis. I also do not want to pass the opportunity to thank several BC Hydro employees who have helped me to understand BC Hydro’s in-house models including Wun Kin, Doug Robinson, Herbert Louie, Jiyi Zhou, Jian Li, Ryan Rasmussen and Leonie Tharratt. I would also like to thank my second reader Dr. Steven Weijs for taking the time from his busy schedule to review the thesis.  Finally, I would like to thank family and friends, both here in Canada and back home in Ethiopia, for their continued support and motivation through my master’s program. 1   Introduction The adoption of wind power as a source of electric energy is increasing at an unprecedented rate around the world. Figure 1 depicts the rise in the total global installed wind capacity from less than 25 GW in 2001 to more than 480 GW in 2016. In Canada, an additional 0.71 GW of wind was installed in 2016, bringing the total capacity of the 285 wind farms and their 6288 individual wind turbines to 12 GW. In the United States, the total installed wind capacity has increased from 4.15 GW in 2001 to 84 GW in early 2017. This rise in installed capacity is observed in all parts of the world, but Asia has experienced the most significant developments in recent years. (American Wind Energy Association, 2017; Canadian Wind Energy Association, 2016; GWEC, 2017).   Figure 1 - Global installed wind capacity Source: GWEC (2017) One of the main reasons for the rise in popularity of wind energy is the attention being given to climate change and greenhouse gas emissions. Wind energy is considered to be more environmentally friendly than conventional power generation sources that require the combustion of fossil fuels.  Except for the emissions resulting from the construction of the wind farms and the production of the turbines, wind energy has no energy production emissions. Wind energy’s impact on the land it is built on, its interference with wildlife habitat, and its production by-products are also minimal when compared to other power generation resources that do not use fossil fuels such as hydropower and nuclear power. In many countries, this has led to several state and federal policies that promote the development of wind energy. In addition, wind energy’s lack of dependency on a commoditized fuel 0501001502002503003504004505002001 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012 2013 2014 2015 2016GW2  source like coal or natural gas and its reliance on a renewable fuel source (i.e. wind) makes it well suited for achieving energy security (IEA 2015; Bird et al. 2005). Wind energy incentives, which are one of the policy actions taken by some governments, have been an essential driving force for increased research and investment in wind energy. For example, in the US, the incentives come in the form of federal tax incentives like the Production Tax Credits, Business Energy Investment Tax Credit, and Advanced Energy Manufacturing Tax Credit. In addition, the Department of Energy offers loan guarantees to help green energy companies secure loans. Several U.S. states also provide additional tax incentives. Renewable Portfolio Standards (RPS), which creates an obligation for utilities to generate a specific percentage of their energy using renewable energy sources, are also one of the policy drivers for the growth of wind. In Canada, the incentives come in the form of investment subsidies, RPS, accelerated capital cost allowance and indirectly through carbon credits and carbon tax (U.S. Department of Energy 2013; Surespanwind 2017). In addition to environmental benefits, another factor contributing to the rising popularity of wind energy is its declining cost. The levelized cost of energy (LCOE) for wind energy, which expresses the total cost of construction and operation of an electric plant divided by the total energy it produces in its lifetime, has decreased significantly since the 1980s. Figure 2 shows the LCOE in the US for the last 35 years. It shows the steady fall in cost from 1980 to the early 2000s, where the cost starts to stabilize. A significant drop in cost is also seen after 2008. In some locations with excellent wind resources and transmission access, the price of wind has become competitive with fossil fuel generation resources. This is especially true for states with wind incentives. The cost reduction is mainly driven by advancements and scale-up of wind turbine technology that has led to improvements in performance and reliability. These advancements have made wind farms an attractive investment even in locations which were previously thought to be unfavorable for wind power. The slight rise in cost seen from the early 2000s to 2008 is explained by rising material prices, increased labor costs, higher profits from manufacturers, and turbine upscaling (GWEC 2017; Schwabe et al. 2011). According to Schwabe et al. (2011), the price is expected to fall going forward and is expected to stabilize around 2030.  3   Figure 2 - Average LCOE in good to excellent wind sites    Source: U.S. Department of Energy (DOE) (2015) 1.1 Background and Motivation Although wind power has become an attractive resource option because of its advantages and declining cost, it has some unique disadvantages as a power generation resource. Wind is an intermittent (i.e., varying minute to minute and hour to hour), uncertain (i.e., difficult to accurately forecast), and non-dispatchable (i.e., only generates electricity when the wind is blowing) power source. The intermittency necessitates the use of capacity reserves from fast-moving generators that are capable of responding and firming up the minute-to-minute fluctuations of wind power generation. The uncertainty in the forecast at different time scales also obligates utilities to withhold generating capacity to manage the forecast error. The non-dispatchability property can also force system operators to curtail wind power generation or spill water from hydropower generating facilities. These three properties can have significant financial impacts on power systems (Dowds et al. 2015; Charles River Associates 2009; BC Hydro 2013). Utilities usually conduct wind integration studies to assess the effects of adding more wind energy to their system and these studies often have sections dedicated to estimating the wind integration cost. This thesis focuses on the financial cost arising from the uncertainty of wind power in the day-ahead (DA) scheduling time frame. The DA wind uncertainty is more concerning to system operators than uncertainties in other timeframes (e.g., the hour ahead timeframe) for two main reasons. First, most 4  unit commitment (UC) and scheduling processes of the power generation system is carried out in the DA timeframe. UC is the process used by operators to schedule the generation units to meet the forecasted load with minimum cost while meeting reliability criteria and system constraints (Greenhall 2013). Most of the UC plans are developed in the DA timeframe because it is easier to acquire transmission access and because some generators have a long start-up time. Market prices are also less volatile in the DA time frame, which reduces market risk (Guzman 2010; Conejo et al. 2010). The forecast uncertainty brought in by wind power alters the optimal UC schedules since more flexible generators (which usually are more expensive) have to be committed and operated to respond to the forecast error. The second reason that the DA uncertainty of wind is more concerning is that it is more significant than the real-time (RT) (i.e., hour ahead or 15-minutes ahead) forecast uncertainty. The average forecast error in the day-ahead timeframe can be orders of magnitudes higher than the average error in the RT timeframe. This is because the forecast confidence and accuracy decrease as the forecast lead time increases (Hodge & Milligan 2011; B.-M. Hodge et al. 2012).  The average DA mean absolute wind power forecast error can be can vary significantly depending on several factors such as wind penetration (i.e., the percentage level of wind capacity in the system), the terrain surrounding wind farms, the geographic diversity of wind farms, and the quality of the forecasting model. High level of uncertainty has a negative impact on the UC. It can lead to the over-commitment of fast reacting units that have high running cost. A significant amount of capacity might also have to be reserved, which reduces energy trading opportunities (Miettinen & Holttinen 2017; Giebel & Sørensen, Poul Ejnar Holttinen 2007; Greenhall 2013). Because of the factors that influence the WFE, the cost of DA wind forecast uncertainty (also called the DA wind integration cost or DA wind uncertainty cost in this thesis) is site specific. For example, the DA wind uncertainty cost has been estimated to be as low as $1/MWh (Pacific Corp 2015), as high as $19/MWh (Idaho Power 2013), or somewhere in between at $3.5/MWh for the Eastern Interconnection (EnerNex Corporation 2011).   The cost of DA wind forecast uncertainty has also been investigated in British Columbia (BC). The BC Hydro and Power Utility (BC Hydro) examined the costs in a wind integration study and reported that the DA wind integration costs ranged from $2.78/MWh to $11.94/MWh depending on the number and geographic diversity of wind farms that could be integrated into the BC Hydro system. 5  The average cost was $6.45/MWh, which was about 9% of the price of electric energy used for the study (BC Hydro 2013). This shows that the DA wind integration cost can be an issue when integrating wind into the BC Hydro system.  Researchers and industry experts have proposed several approaches to minimize the effects of wind power uncertainty. Commonly suggested solutions include improved forecasting, stochastic optimization, more flexible generating fleet, more coordination between balancing authorities, sub-hourly scheduling, rolling planning (i.e. scheduling units more frequently with updated forecast), developing wind farms in geographically diverse areas, energy storage, market restructuring, hedging losses with financial tools, transmission improvements, and  better demand response (Charles River Associates 2009; Greenhall 2013; Bruninx 2016; Idaho Power 2013; Conejo et al. 2010; GE Energy 2010). The focus of this thesis will be on the stochastic optimization approach. Stochastic optimization (SO), sometimes called stochastic programming, is a type of optimization algorithm where the uncertainty of the stochastic parameter(s) is explicitly modeled. It can be used to generate optimal bidding strategies that consider the uncertainty in the wind, load, and price forecasts. Several researchers have investigated SO’s benefit in minimizing the wind integration cost. Rahimiyan et al. (2011) have shown that wind power bidding strategies from stochastic optimization models generally outperform bids based on deterministic forecasts. Matevosyan and Söder (2006), Morales et al. (2010), Angarita et al. (2009), Dai and Qiao, (2015), Vagropoulos et al. (2016) and several other researchers have also arrived at the similar conclusion that SO models provide robust solutions that minimize the cost of wind uncertainty.  Some researchers have estimated the percentage of the uncertainty cost that can be saved by applying stochastic optimization. The term stochastic saving is used in this thesis to describe the percentage of the cost that can be saved by using SO instead of deterministic optimization. Stochastic saving is usually estimated by comparing three different optimization models: a deterministic model with a perfect forecast assumption, a deterministic model without a perfect foresight assumption, and a stochastic model. The perfect forecast assumption provides the lower bound of the cost while the deterministic model without the perfect forecast assumption provides the upper bound. The stochastic model’s cost is usually between the two. Tuohy et al. (2008) reported that 35% of the cost of uncertainty could be recovered with stochastic optimization, Greenhall (2013) reported savings 6  around 25% and Ruiz et al. (2009) estimated that only 12% savings could be expected. The variation in stochastic savings indicates that the savings are dependent on many factors, such as the electric system modeled and the accuracy of the statistical characterization of the uncertain parameter(s). Uncertain parameters in SO problems are usually characterized by scenarios that are obtained from a scenario generation process.  Scenario generation is the process of creating a limited set of representative outcomes of an uncertain parameter. This is one of the main challenges modelers face when they decide to use SO. Creating a limited set of possible outcomes (i.e., scenarios) that can accurately capture the properties of stochastic (i.e. uncertain) parameters is not a trivial task and it is one of the obstacles that limit the development and use of stochastic UC scheduling algorithms (Greenhall 2013). Therefore, scenario generation is the focus of many research topics in the field of stochastic optimization. Creating a representative set of scenarios is crucial because the cost saving expected from SO (i.e. stochastic saving) is dependent on the quality of the information captured in the scenarios. For example, the two scenario generation methods used in Greenhall (2013) resulted in significantly different stochastic savings. Other researchers including Conejo et al. (2010), Morales et al. (2010) and  Rahimiyan et al. (2011) have also demonstrated that the benefit of SO is dependent on the quality of information that is captured in the scenarios.  Two of the commonly adopted modeling approaches in the SO literature to generate wind power related scenarios are auto-regressive moving average (ARMA) based models and Markov chain (MC) based models. These models have been found to be flexible and capable of capturing the unique statistical properties of time series data such as wind power and wind power forecast error time series. These two approaches require the generation of a substantial number of scenarios to guarantee a good representation of the underlying stochastic process. The caveat of generating numerous scenarios is that the size of the optimization problem scales with the number of scenarios used. Some large-scale optimization problems might require a high computational power, a long solution time, or become intractable (i.e. unsolvable) when formulated as a SO problem with many scenarios. One of the commonly used approaches to deal with this problem is scenario reduction (Conejo et al. 2010; Bruninx 2016). 7  Scenario reduction is the process of reducing the cardinality of scenarios generated by a scenario generation algorithm. It works by deleting scenarios that have similar properties based on some predefined metric and transferring the probabilities of the deleted scenarios to the preserved scenario that are most similar to them. Scenario reduction is an attempt to find a compromise between solution quality and computational burden. The algorithms proposed in Dupačová et al. (2003) and Heitsch and Römisch (2003) have been widely used by many researchers in the field of power system optimization. This thesis will investigate appropriate scenario generation and reduction methods that can be used in the SO model developed in this research.  1.2 Goals and Objectives The main goals of this thesis are to investigate two scenario generation algorithms that can be used to accurately capture the statistical properties of the wind forecast error and to determine the savings in the DA wind integration cost that can be expected from using SO. More specifically, the main goals of this thesis are: 1. To characterize the statistical properties of the wind forecast error. 2. To propose and compare two scenario generation approaches, one based on MC and the other based on ARMA, for creating scenario sets that can accurately capture the properties of the wind forecast error.  3. To quantify the value of stochastic optimization in the BC Hydro system given the system’s numerous operational constraints, the level of wind power uncertainty that is unique to the region and the generation resource mix that is dominated by hydropower.   4. To assess the effect of the scenario generation technique on the optimization result. The ARMA and Markov chain based models will be used to generate scenarios for the SO and the results will be compared. 1.3 Organization of the Thesis The thesis is presented in six chapters. In this chapter, the motivation and goals of this thesis are discussed. In the next chapter, literature that is found to be relevant to this thesis is thoroughly reviewed. In Chapter 3, the modeling methodology applied in this thesis is described. The proposed scenario generation approaches are covered in detail in this chapter.  In Chapter 4, a case study that 8  is used to test the proposed models on the BC Hydro system is presented. The results of this case study are presented in Chapter 5. In Chapter 6, the important observations and conclusions of the thesis are presented along with recommendations for future work. 9   Literature Review This section reviews the relevant literature on the following fields: development of wind power in BC, wind forecast error properties, wind integration studies, scenario generation methods, scenario reduction methods, and stochastic optimization. The focus of the review will be on the state-of-the-art methods applied to management wind integration in large-scale power systems with particular emphasis on hydropower systems. 2.1 Wind in BC The province of British Columbia (BC) is rich with renewable energy resources. Almost 95% of BC’s electricity is generated from renewables. The capacity of the electric power system in BC exceeds 16 GW, 90% of which is renewable. Hydropower is the most important source of power generation, followed by biomass, natural gas, wind, and diesel. Two provincial Crown corporations, British Columbia Hydro and Power Authority (BC Hydro) and the Columbia Power Corporation own and operate approximately 80% of the installed capacity in the Province. The remaining capacity is provided by FortisBC and independent power producers including municipalities and First Nations. BC Hydro is responsible for most of the installed capacity and is the primary electric power producer in BC (National Energy Board 2016). The BC Energy Plan which was issued by the Ministry of Energy, Mines, and Petroleum Resources in 2007 sets out conservation, energy efficiency, and clean energy goals. The Plan also provides strategies and policy actions to be followed to achieve these goals. Some of the policy actions set that all new electricity generation projects should have zero net greenhouse gas emissions, that at least 90% of the electricity produced in BC must continue to come from clean or renewable sources and that the procurement of electricity appropriately recognizes the value of intermittent resources. The plan’s definition of clean or renewable resources includes sources of energy that are renewed by natural processes, such as water, solar, wind, tidal, geothermal, wood residue, and energy from organic municipal waste (MEMPR 2007). In accordance with the Energy Plan, BC Hydro issued a Clean Power Call and awarded several contracts to wind projects (BC Hydro 2010). As of June 2017, five utility-scale wind farms are operating in BC with a total nameplate capacity of 672 MW. The first farm to start operation, the Bear Mountain Wind Park, began operation in 2009, and the last wind farm to 10  come online was the Meikle Windfarm, started operating in 2017. The five operating wind farms in BC are listed in the table below along with their location, capacity, and the start date of their commercial service.  Table 1 -  Wind Farms Operating in BC Wind Farm Location Start date of Commercial Operation Nameplate Capacity (MW) Bear Mountain Wind Park Near Dawson Creek (Peace Region) October 2009 102 Dokie Wind Farm Near Chetwynd (Peace Region) February 2011 144 Quality Wind Farm Near Tumbler Ridge (Peace Region) November 2012 142.2 Cape Scott Wind Farm Near Port Hardy (Vancouver Island) November 2013 99 Meikle Wind Farm Near Tumbler Ridge(Peace Region) January 2017 184.6  The installed capacity of wind power is expected to increase in BC since BC Hydro considers wind power to be one of the resource options available to meet BC’s energy deficiency. The BC Energy Plan and BC Hydro’s vision state the goal of achieving energy independence by being capable of producing adequate energy in the province. BC Hydro’s latest Integrated Resource Plan (IRP)  (BC Hydro 2013) was released in 2013 and it forecasts that there will likely be an energy deficit in the near future. BC Hydro is continuously investigating several resource options to meet this shortfall while meeting BC Energy Plan’s requirements. Wind is one of the resource options that is investigated since it offers competitive price and satisfies the requirements set for new power plants.  2.2 Wind Integration Cost It is important to properly evaluate the cost of wind power when considering it as a resource option. This evaluation involves the estimation of the wind integration cost. Most of the cost of integrating wind into an electric system arises from two fundamental properties of wind generation: variability and uncertainty. The variability of wind induces some operational added cost even under a hypothetical situation where perfect knowledge of the future is assumed. This added cost comes from increased regulation/following reserve requirements that are needed to manage the variability of wind and increased ramping/cycling of conventional generation units to balance its variability. The uncertainty of wind power causes added cost by increasing the operating/imbalance reserves required 11  to manage the WFE. This uncertainty can cause a significant additional cost in the hour-ahead and DA time frames (Milligan et al., 2013).  As pointed out in Milligan et al. (2013), isolating the wind integration cost is challenging. For example, one of the preferred methods of assessing wind integration cost is to run a production costing model with and without wind power and compare the costs. If a conventional generation provides the power in the 'without wind' case, the added cost might not be representative as it could understate the benefits of wind power such as the water, fuel, and emission savings. New wind farms may also cause a decrease in capacity factors for some of the conventional generators since wind is usually given priority for power generation because of its lower running cost and the availability of incentives like Renewable Energy Credits. This reduction of capacity factor might increase the total cost and LCOE for some of the other generation units. The effects of integrating wind into the system are not limited to financial costs. Wind power can affect the system reliability, grid stability, system efficiency, transmission congestion, etc. The total effect of wind on the system depends on the wind penetration level, the geographic diversity of the wind farms, type of other generating units in the system, and the overall size of the system. The negative effects tend to be higher on the smaller power system and tend to increase as the wind penetration level increases. Diversified wind farm locations and flexible generation resources like large hydropower plants usually reduce the impact of wind power (Holttinen et al. 2009). Some of the impacts of wind power on a power system for time and area scales are shown in Figure 3. 2.2.1 The BC Hydro Wind Integration Study BC Hydro conducted a wind integration study in 2010 to isolate the financial impacts of wind power integration at different penetration levels. The results of this study were used to compare the financial attributes of wind power with other resource options. Some of the other resource options that were compared with wind power include biomass, run-of-river hydro, large hydro, geothermal, natural gas, and solar power. These resource options were compared and ranked on four attributes: technical, environmental, economic, and financial. Technical attributes refer to resource characteristics like dependable generating capacity and average annual energy generation. The environmental attributes evaluate the effects the resource options have on land, atmosphere, and marine life. The economic  12  development attributes of a resource option consider the contribution of the project to the local economy. Financial attributes mainly deal with the levelized cost of resource options as unit energy cost (UEC) in dollars per megawatt hour ($/MWh) and unit cost of capacity (UCC) in dollars per megawatt ($/MW) (BC Hydro 2013). When estimating the financial attributes (i.e. UEC and UCC) of a resource option, the non-dispatchability cost and the integration cost are not included. UEC and UCC mostly focus on the investment and financing cost. A separate study is conducted to investigate the integration costs. The final cost that is used to rank the financial performance of resource options is the sum of all these costs and the main reason BC Hydro conducts the wind integration study is to determine the total cost of wind power.   Figure 3 - Effects of wind power in the power system           Source: Holttinen et al. (2009) The BC Hydro wind integration study was conducted for two future study years (2011 and 2021) with forecasted loads and six wind scenarios (see Table 2). The scenarios were developed based on the combination of two wind farm selection criteria labeled ‘CAPEX’ and ‘high diversity’ for three wind penetration levels: 15%, 25%, and 35%. Penetration levels are measured by dividing the wind installed capacity by the peak demand in the corresponding load year. In the CAPEX scenarios, the wind farms were selected based solely on their UEC, which is mostly based on the capacity factor 13  (i.e., average generation/total capacity). Since the highest capacity factor in BC is found in the Peace region, most of the wind farms selected in the CAPEX scenarios are in the Peace region. High diversity scenarios were also considered under the assumption that high diversity lowers the variability and uncertainty of wind power. For these scenarios, equal installed capacity was distributed across four different regions in BC that experience different weather patterns.  BC Hydro (2013) identified two types of wind integration costs: operating reserves cost and day-ahead opportunity costs. The operating reserve cost resulted from an added requirement of the three types of operating reserves (regulation, following, and imbalance). BC Hydro forgoes the flexibility of trading the capacity withheld for these reserves in electricity markets. The DA opportunity cost arises from the loss of flexibility when trading in the DA electricity market with DA wind power forecasts that are uncertain.  This uncertainty can cause reliability issues and therefore is managed by keeping reserves that could have been traded in the DA electricity market.   The results of the BC Hydro wind integration study is shown in Table 2.  The table shows that the reserve cost ranged from $3.24/MWh for the high geographic diversity scenario at 15% penetration to $7.68/MWh for the CAPEX scenario at 25% penetration. The table also shows that the DA opportunity cost ranged from $2.78/MWh for the high geographic diversity scenario at 15% penetration to $11.94/MWh for the CAPEX scenario at 25% penetration. These costs are significant given that the reference electricity price for the study was $74.45/MWh. On average, about 55% of the wind integration cost was attributed to the DA uncertainty of the wind forecast. The table also shows that the geographic diversity decreases the wind integration cost. It can also be seen that the cost of DA wind uncertainty was more affected by the geographic diversity than the cost of wind variability. The CAPEX 35% scenarios have lower integration cost than the 25% because wind farms from different regions were added to achieve a penetration level of 35%.     14  Table 2 - Summary of wind integration cost for the BC Hydro system Source: BC Hydro, (2013)  2.2.2 Review of Other Wind Integration Studies Numerous wind integration studies (WISs) have been conducted by regional utilities, government agencies, and interested bodies. The study methodologies used in these WISs are not usually consistent since each study is designed to assess the integration impacts for a particular region, electricity market, utility, energy source mix, etc. Dowds et al. (2015) reviewed several large-scale wind integration studies conducted in the US and found some patterns in the methodologies used. Some of the common steps in conducting WISs include: - Collection of wind data from numerical weather prediction (NWP) models, anemometer readings and historical wind plant outputs - Collection of historical load data - The assumption of Gaussian distribution to represent the wind power variability and uncertainty - Using power system models to investigate the effect of wind power on the operation and reliability of the transmission system - Reporting the findings of the study. Wind integration cost is usually reported in dollars per megawatt-hour ($/MWh).  15  Some large-scale WISs are discussed below, with more emphasis given to those which contain DA wind integration costs. a. Eastern Wind Integration and Transmission Study EnerNex Corporation (2011) carried out a detailed wind integration study for the Eastern Interconnection system. The Eastern Interconnection covers a large land area in the US and serves a significant portion of the county’s population. It covers most of the mid-west, eastern and southern states of the US. The study considered four wind scenarios, three of which have 20% penetration and one scenario that has 30% penetration. The scenarios with the same penetration differed from one another in the geographic siting of the wind farms. The reference scenario which represented the wind capacity at the time of the study had a 6% wind penetration. The study concluded that a substantial portion of the wind integration cost is caused by the large-scale variability and uncertainty of wind. The integration cost for the 20% penetration scenario was estimated to range from US$8/MWh for the low diversity scenario to US$5.77/MWh for the high diversity scenario. The 30% penetration scenario, which is the combination of the above two scenarios (low and high diversity), had an integration cost of $7.07/MWh. The DA uncertainty was found to be the cause of 28% of the integration cost for the low diversity scenario while it increased to 50% for the high diversity scenario.  b. Western Solar and Wind Integration Study  GE Energy (2010) conducted a similarly extensive wind and solar integration study for the WestConnect group of utilities. WestConnect is a collaborative group of transmission providers that operate in the Southwestern United States. The study was conducted for different wind penetration scenarios ranging from 3% (the penetration at that time, base case) to 30%. The results showed that at 30% penetration, fuel and emission costs were reduced by 25-45% when compared to the base case. For the studied system, the cost of wind variability was discovered to be minimal since the reduced net load (i.e. load less wind, the amount of load that has to be served by conventional generation) meant that more capacity resources were made available to provide operating reserves. The cost of uncertainty was found to be much higher. The average system mean absolute forecast error (MAE) for the different wind scenarios was about 8% of the installed capacity, which increased the operating cost by $13-22/MWh depending on the scenario investigated. Using the state-of-the-16  art forecast for the renewable energy sources in the DA unit commitment process was shown to reduce operating cost by $12-20/MWh. The perfect forecast was shown to reduce the cost further by $1-2/MWh. Reserve shortfalls and wind curtailment were observed, especially for the 30% penetration scenario.  c. Pan-Canadian Wind Integration Study GE Energy (2016) performed a Pan-Canadian wind integration study for the Canadian electric system. The purpose of the study was to identify the operation impacts, costs, and benefits of integrating substantial amounts of wind power into the Canadian system. The study was conducted for 4 penetration scenarios ranging from 5% (the base case) to 35%. The study concluded that 35% penetration is possible in Canadian electric system with some transmission improvements and more regulation reserves additions. In the study, it was assumed that hydro units are scheduled based on the DA load and wind forecasts and that these schedules were fixed for real-time operation. Under these assumptions, the WFE was shown to have a significant cost. The study showed that re-dispatching available flexible hydro units in the real-time operation could counteract the effects of the WFE, reducing the system cost by up to C$228M in the Eastern Interconnection (for the 20% scenario). The cost estimates and savings reported in this study do not consider the forgone market opportunity of selling the generators’ flexible capacity in the ancillary markets.  d. SPP WITF Wind Integration Study Charles River Associates (2010) studied the operational and reliability impacts of wind integration on the Southern Power Pool (SPP). The SPP is a regional transmission organization that operates in eight Southwestern states. The study considered four wind penetration scenarios at 4% (base case), 10%, 20%, and 30% penetrations. The conclusion regarding the DA market was that integration of wind does not necessarily result in lower market prices and that the WFE increases the wind integration cost. Some of the suggestions given to minimize the negative market impacts include comprehensive forecasts and the use of stochastic unit commitment. e. Idaho Power Wind Integration Study  17  Idaho Power (2013) conducted a wind integration study to investigate the expected cost that arises from integrating 800 MW to 1200 MW of wind power into the Idaho power system. At the time of the study, Idaho Power had a total wind capacity of 678 MW. Hydro and coal generation dominate the Idaho Power system. Since the coal-fired plants have low flexibility, Hydro plants were assigned the job of managing the variability and uncertainty of wind power. The wind reserves required to manage the DA WFE were calculated as the capacity required to manage 90% of the WFE. The study arrived at an average wind integration cost of $8.06/MWh for 800MW, $13.06/MWh for 1000 MW and $19.01/MWh for 1200 MW. This cost does not include the cost associated with wind curtailment.  The cost was estimated by running a production model with and without incremental DA wind reserves. f. Pacific Corp Wind Integration Study Pacific Corp (2015)  conducted a WIS for the thermal power dominated Pacific Corp system using their Planning and Risk (PaR) model. Pacific Corp is an electric power company that operates in the Northwestern United States. The result of the study estimated the total wind integration cost of $3.06/MWh. The report stated that 77% of the increase in operational cost in the Pacific Corp’s balancing areas was caused by the increase in incremental regulating reserves that are needed to manage short-term deviations. The rest of the integration cost was the result of the DA WFE.  In summary, the wind integration studies reviewed in this section show that the estimated wind integration costs vary significantly by region and by electric system. The integration cost was shown to vary for different utilities and interconnections depending on many factors like generation fleet, geography, balancing area size, electricity market structure, wind integration cost estimation methodology, etc. For example, the wind integration cost reported by Pacific Corp (2015) is much lower than the one reported by Idaho Power (2013). 2.3 Day-Ahead Wind Forecast Error As mentioned in Chapter 1, one of the commonly suggested approaches to minimize the negative impacts of wind uncertainty is improved forecasting. Several forecasting techniques and tools are continuously being investigated by researchers and utility companies who wish to improve the forecast accuracy. Physical, statistical and combined approaches have been studied, with different 18  approaches performing better than others at different jurisdictions and under different settings, such as the forecast time horizon or extreme event forecasting (Jung & Broadwater 2014).  Although some progress has been made in reducing the wind forecast error, completely eliminating it and having a perfect forecast is unachievable. Therefore, the magnitude and property of WFE is an important topic that has been the focus of many studies. Some of these studies are reviewed below.  Bludszuweit et al. (2008) argued that, contrary to what is usually assumed in many wind uncertainty studies, the WFE is not normally distributed. In this study, persistence forecasting was used to create a one-year long sub-hourly WFE time series from recorded wind power generation time series. Persistence is a forecasting method that assumes the wind for the next timestep is the same as the wind at the current timestep. The authors claim that persistence forecasting would create a WFE distribution with comparable properties to other more sophisticated forecasting methods. The inaccuracy of assuming a normal distribution in approximating the WFE was shown to be greater at the tails of the distribution. The authors instead fitted a newly proposed Beta probability density function (PDF) to the WFE time series. The Beta PDF was chosen over the normal distribution because it has variable kurtosis that can be specified. The authors proposed new approximation methods to estimate the Beta distribution parameters and showed that the Beta pdf is better at capturing the fat-tailed WFE distribution. Their work focused on the within-the-day forecast (i.e., forecast horizon less than 24 hours) but their approach is also applicable for the DA WFE.  Tewari et al. (2011) showed that the WFE distribution resembles the normal distribution when the WFE data is collected from wind farms that have geographic diversity. However, the assumption of normality does not hold for individual wind farms. The authors stress that the WFE distribution is not continuous since small errors in the wind speed forecast at speeds lower than the cut-in speed or higher than the steady wind speed region in the wind power curve (refer to Figure 4) have no impact in the power generation forecast. This is because the wind power curve is horizontal at these wind speeds. To capture this non-continuity, the authors suggested a mixed distribution function that models the continuous parts of the distribution using the Laplace distribution and assigns the discrete points (i.e. low and high speeds) a probability of zero.  19   Figure 4 – A typical wind power curve with commonly used terminology    Source: http://www.wind-power-program.com/popups/powercurve.htm Hodge and Milligan (2011) investigated the WFE distribution at different time scales ranging from a few minutes to three hours. The Cauchy distribution was selected to model the WFE since it can capture the high kurtosis of the WFE. The WFE has high kurtosis because most of the variance is caused by few large errors instead of many small errors. The properties of the WFE distribution, especially the kurtosis, were found to be heavily dependent on the forecast timescale.  The WFE distributions at smaller time steps (less than 15 minutes) exhibit much larger kurtosis values than timescales of more than an hour. The Cauchy distribution was successfully applied to create confidence bands around the wind forecast for all the time scales considered and it was shown to provide more reliable outputs when compared to the normal distribution.  Hodge et al. (2011) investigated the WFE distributions of different countries at DA and hour-ahead timeframes. The real WFE data that was gathered from several sources was fitted to a hyperbolic distribution, a distribution that can accommodate the heavy tails and skewness of the observed WFE data. The study showed that the DA WFE has a much lower kurtosis than the hour-ahead WFE, confirming that the hour-ahead forecasts are much more accurate than the DA forecast. With regards to installed capacity and geographic diversity of the wind farms, the overall trend was that higher 20  installed wind capacity reduces the variance of the WFE distribution. This is mostly because higher wind capacity indicates more geographic diversity of wind farms, which means that various wind farms are subjected to different and uncorrelated weather conditions. Hodge et al. (2012) further compared the load forecast error (LFE) and WFE distributions in the DA timescale for the New York and California independent system operators (ISOs). In the comparison, it was shown that both LFE and WFE distributions are leptokurtic (i.e., have higher kurtosis than the normal distribution) and that the hyperbolic distribution is suitable to model them. The most significant difference between the two is that the LFE has a much smaller range of values and much higher kurtosis values, which indicates that it has lower forecast error magnitudes when compared to the WFE.  Miettinen and Holttinen (2017) investigated the characteristics of DA WFE in Nordic countries. The results showed the value of geographic diversity in reducing the WFE. The study found that the correlation of the DA WFE from different wind farm decreased as the distance between the farms increased. The mean absolute error (MAE) of the DA WFE was shown to drop from 8% when a single wind farm was analyzed to 2.5% when the aggregated wind generation of the Nordic region was analyzed. The maximum observed WFE also decreased from more than 80% to 13.5%.  The effects of the WFE statistics on the operational cost of an electric system was discussed in Lowery and O’Malley (2012). In order to isolate the most critical statistical parameters of the WFE, the authors simulated several sets of WFE scenarios with different statistical characteristics by using a flexible, moment matching based scenario generation tool (see Section 2.4.1 for details). They controlled the property of the WFE by varying the first four moments (i.e., mean, variance, skewness, and kurtosis) of the distribution in the scenario generation tool. Then, a stochastic optimization algorithm that minimizes the system cost was run with these different sets of scenarios.  Some of the observations made in their study suggest that the system cost is the lowest when all the four moments are considered in the optimization and that considering the first two moments only gave better results than considering the first three moments. This indicates that the use of skewness and kurtosis should be carried out together, as including skewness without kurtosis distorts the representation of the stochastic process. They also investigated the effect of the accuracy of the stochastic information on the solution quality by intentionally overestimating and underestimating the moments. It was shown that the system cost is generally dependent on all moments and that there is an interaction between the accuracy of the moments (i.e., the effect of underestimating one moment can be countered or 21  exacerbated by overestimating another moment). For example, an underestimation of the kurtosis can offset an overestimation of variance.  The authors also studied the solution dependence on the variance only by fixing the skewness and the kurtosis with their true values but by varying the accuracy of the variance. The system cost was seen to increase when the variance was overestimated and underestimated, but the cost rose at a higher rate when the variance was overestimated by more than 25%.  In summary, the studies discussed above show some of the important properties that must be considered when analyzing the WFE. The studies have established that the WFE property is a function of many factors such as geographic location, number of wind turbines, and forecast time scale. The facts that the WFE increases with forecast lead time, decreases with geographic diversity, and that it has a leptokurtic distribution have been determined from the literature.  2.4 Scenario Generation As discussed earlier in Chapter 1, uncertainty in stochastic optimization (SO) models can be represented by a set of probabilistic scenarios. Scenario generation is the process of creating these scenarios. When wind is the uncertain parameter in a SO model, it is treated as a stochastic process. A stochastic process, sometimes called a random process, is a collection of random variables that are realized (or indexed) over time. Examples of other stochastic processes in the context of energy systems include commodity prices, reservoir inflows, and electricity demand. A scenario is one possible realization of this stochastic process. During scenario generation, these possible realizations are created, and probabilities are assigned to them. If the stochastic parameter is discrete (e.g., the number of generators available), the possible outcomes of that random variable are finite and therefore less effort is usually required to characterize the process. However, for continuous stochastic processes (like wind speed or market prices), the possible outcomes are infinite, which makes their accurate representation a challenge. Since computational requirements limit the number of scenarios that can be used in a SO algorithm, a finite set of scenarios has to be produced. Care must then be taken when converting a continuous stochastic process to a non-continuous one to preserve its statistical characteristics. Many researchers have shown that the accuracy, reliability, and savings from SO models increase when the quality of the information contained in the generated 22  scenarios increases (Keutchayan et al. 2016; Sharma et al. 2013; Conejo et al. 2010; Bruninx 2016; Rahimiyan et al. 2011).   Sets of scenarios for multi-period problems are usually presented in the form of a ‘scenario tree’ or ‘scenario fan’. A scenario tree (refer to Figure 5) consists of a finite set of nodes and it starts from the root node in the first period, where it is assumed that there is no uncertainty. From there, the tree branches into nodes at the next period. Each branch has its own probability of occurrence and realization value. The branches then form other nodes and continue to branch out to the next time step as shown in Figure 5. The sum of the probabilities of all the branches originating from a node is one. The branching continues up to the final timestep (Latorre et al. 2007). The number of scenarios in a scenario tree increases exponentially with the number of timesteps (stages) considered. A scenario fan is different as several scenarios branch out from the root node and follow a single path up to the final timestep (stage) without subsequent branching at the intermediate time steps (stages) or nodes. This means the total number of scenarios in a scenario fan will remain constant throughout the timesteps. Scenario fans are more suited for two-stage SO problems while scenario trees are more suited for multi-stage SO problems. The difference between the can be seen in Figure 5.   Figure 5 - Scenario tree (on the left) and scenario fan (on the right) Source: Hibiki (2006) Some wind power or wind speed scenario generation techniques are adopted from algorithms whose primary purpose is to generate stochastic wind time series data for other purposes. Stochastic wind time series that accurately capture the unique properties of wind power are required for informed 23  planning and simulation of an electric system that has a significant amount of wind power. Therefore, the creation of stochastic wind speed or wind power time series is also an active area of research. The algorithms that generate these time series can produce as many unique wind realizations as desired, which make them ideal tools for creating scenarios. Some of the popular methods for the generation of stochastic wind data include approaches based on autoregressive moving average and Markov chain models. Therefore, relevant literature discussing the creation of stochastic wind time series is reviewed here even if scenarios are not explicitly mentioned since the methodology can be applied to generate wind scenarios.   The general procedure of using a stochastic time series generation method is to generate wind scenarios by repeatedly sampling from a stochastic model for the same time horizon (e.g., 24 hours). Each sample produces a unique time series of wind speed or wind power. The sampling is repeated until the underlying distribution of the wind speed or power generation is confidently captured. All the realizations can be considered scenarios with equal probability of occurrence. Since the size of the stochastic optimization problem scales with the number of scenarios, it is often preferred to reduce the number of scenarios without losing much of the information they contain.  This might necessitate the use of scenario reduction techniques that aim to reduce the cardinality of scenarios while preserving the characteristics of the stochastic parameters. Some of the scenario reduction techniques that are used for this purpose include moment (property) matching, clustering, and probability distance based methods. It is worth noting that some of the algorithms like moment matching can be used for the generation and reduction of scenarios. Explanations of various scenario generation techniques along with their comparison can be found in Löhndorf (2016).  In the following sections, the literature on scenario generation and reduction is reviewed. More focus is given to the methods that are found to be relevant to the study conducted in this thesis.  The algorithms that can be used for both scenario generation and reduction are discussed in the section that represents the purpose they are most commonly used for. 2.4.1 Review of the Stochastic Wind Data and Scenario Generation Techniques Several scenario generation approaches have been used in the literature to generate wind power or wind speed scenarios. The statistical properties of a stochastic process (i.e. wind speed, power, or 24  forecast error) that researchers attempt to replicate when creating scenarios include its statistical moments, probability distribution and autocorrelation function (ACF). The ACF measures the similarity of a given time series with a time-lagged version of itself. Stochastic time series such as wind speed that have high autocorrelations tend to have realizations that persist for multiple time steps. This means that if the wind speed is high in the current timestep, it is likely to stay high for the next timestep.  The approaches that have been used to create scenarios while preserving these properties include approaches based on Markov chain, ARMA, moment matching, and machine learning.  a. Markov chain models Models based on Markov chain (MC) are one of the widely adopted approaches for the generation of stochastic wind data. MC is a stochastic process where the possible realizations at a given stage can be inferred from information about the previous stage(s). In the literature discussing the generation of stochastic wind scenarios, these stages are usually discretized time steps. MCs are classified based on their order, which is the number of previous stages considered when constructing the models. First-order MCs have the property that the realization at any timestep depends only on the previous timestep, with no dependence on the realizations prior to that timestep. Similarly, second-order MCs need information from the previous two timesteps only. Generally, in an nth order MC, the realization of the process at timestep t depends on the realizations of the n timesteps (i.e., t-1, t-2...t-n). MC modeling of a continuous stochastic process usually requires grouping its continuous outcomes into a limited number of discretized states. The MC model is then parameterized by creating a transition probability matrix that quantifies the probability transitions between these states (Shamshad et al. 2005; Tagliaferri et al. 2016). Several researchers have used Markov chains to generate synthetic wind speed time series (WSTS). Sahin and Sen (2001) used a first-order Markov chain to produce a WSTS for locations in Turkey with high wind energy potential. Their work showed that first-order MCs can sufficiently preserve most of the statistical parameters and that second-order MCs provide improved results. Nfaoui et al. (2004) used an MC model to generate stochastic WSTS for a location in Morocco. The model was successful in reproducing most of the statistical characteristics, such as the marginal distribution of the observed wind time series, but failed to replicate the ACF.  Shamshad et al. (2005) used first and second order MC models to generate WSTS and compared the results. They concluded that both 25  models were able to preserve most of the statistical properties of the observed WSTS, with the second-order model yielding a slightly better result. The authors of this study were also unable to replicate the autocorrelation and spectral density properties of the observed data. They suggest that removing the seasonal component of the time series could help to solve the problem. It is apparent from the studies discussed above that one of the major drawbacks of using MCs is their underestimation of the autocorrelation factor. Brokish and Kirtley (2009) further investigated the pitfalls of using MCs to generate wind power time series (WPTS) synthetically. For this study, the authors used data from a single wind turbine to train first, second and third-order MCs, which were then used to generate WPTSs data at different temporal resolutions. Then the properties of the generated and original WPTSs were compared. The results showed that although PDF of the wind power was accurately captured, even the third order MC failed to display the same persistence as the original data. The authors concluded that MC models are not appropriate for modeling a process whose time evolution is important. They stated that the shortcomings are especially noticeable with wind data when the time steps simulated are less than 15-minutes in length.  Several researchers have provided recommendations to improve MC models. Hocaoglu et al. (2008) suggest the use of higher number of states (i.e., finer discretization) to increase the quality of the generated wind data. In their study, they showed that an MC model with 26 states yielded a better statistical characterization than a model with 13 states. Tagliaferri et al. (2016) proposed a ‘nested’ MC to better preserve the ACF. The proposed model generates WSTS at 1-second timestep using the hourly average and the 1-sec wind speed records from the past hour. The authors claim that this approach helps to solve the inaccurate autocorrelation replication of MC models.  Papaefthymiou and Klockl (2008) suggest that working in the wind power domain yields better results than working in the wind speed domain.  They provided several reasons why this is the case. When working in the wind speed domain, the wind speeds are first generated using MCs and then converted to wind power using power curves. The authors state that this two-step process unnecessarily complicates the model and increase inaccuracies. The number of wind power states required for the MC is also lower since a many wind speed states that contain high and low wind speeds (i.e., less than the cut-in speed and after the steady wind speed, see Figure 4) can be grouped into the same wind power state. The lower persistence observed in the WPTS also enables the use of lower order MCs. 26  Some of the tests conducted in the study indicate that higher order Markov chains are better at reproducing the ACF of the observed data. The authors note that MCs are better at capturing the PDF of the stochastic process than ARMA based models, which makes them preferable in situations where the PDF is more important than the ACF. b. Auto-Regressive Moving Average (ARMA) models Auto-regressive moving averages (ARMA) based models have also been widely used for the creation of stochastic wind data. ARMA models combine the features of moving average and autoregressive models to capture the relevant properties of a stochastic process. 𝐴𝑅𝑀𝐴(𝑝, 𝑞) models are represented by p, the number of autoregressive parameters, and q, the number of moving average parameters. One of the main advantages of ARMA models is their ability to explicitly capture the autocorrelation property of a stochastic time series. Therefore, they are extensively used in disciplines like finance and climate science for forecasting, simulation, and creation of autocorrelated time series. Another reason for their popularity is their flexibility. ARMA models can be easily modified to work with seasonal and multivariate data or to incorporate external regressor variables. Autoregressive integrated moving average (ARIMA) models, which are adjusted to accommodate non-stationary data, have especially seen widespread use by many researchers (Morales et al. 2009; Conejo et al. 2010).  Several researchers have used ARMA based models to create wind data for spatially correlated wind farms. Soder (2004) proposed a multi-variate ARMA model to simulate outcomes of wind speeds for correlated farms based on the available forecast. Actual wind speed forecast data was not available for the sites considered in the study thus persistent forecasting was assumed. The spatial correlation of different wind farm generations was replicated by adding correlated random variables to the various ARMA models built to represent individual wind farms. Instead of simulating the possible outcomes of the wind speed directly, it was indirectly calculated by first creating WFE scenarios and then adding it to the available forecast. The author claims that the method can also be applied to create scenario trees by using the same simulation outcome up to a certain point in time and then simulating different outcomes again from that point on (for example, using one outcome for 1-5 hours then simulating three outcomes at the 5th hour and so on).  The authors did not discuss the how the validity of the 27  models was investigated, and the relationship between the forecast magnitude and the wind forecast error was also not discussed. Boone (2005) extended the model presented in Soder (2004) to work with real wind speed and forecast data. He applied the ARMA (1,1) process to simulate the wind forecast error in the German and Nordic power systems. In this project, wind speeds were first produced and then converted to wind power predictions. The study’s results showed that ARMA process can accurately simulate WFE and that the accuracy of the simulation depends on the source of wind data. The author has found that the WFE variance is dependent on the wind speed, but its relationship with the wind speed was not quantified in the study. The conclusion made regarding this relationship was that assuming a constant variance is adequate for simulations runs carried out under the most common wind speed settings.  Synthetic wind speed generation without consideration for forecast was discussed in Morales et al. (2010). Similar to Soder (2004) and Boone (2005), the authors used an ARMA model with spatially correlated white noise terms to create a dataset of spatially and temporally correlated wind speeds. However, in this study, a probability transformation was carried out to improve the accuracy of the model outputs. ARMA models assume that the data used for model fitting is normally distributed. To correct for this invalid assumption, the wind speeds, which were assumed to follow Weibull distribution, were transformed into normal distribution before model fitting. The outputs of the ARMA model, which are also normally distributed, were then transformed back to Weibull distribution. The authors claimed that this technique generated statistically consistent scenarios that preserved the spatial and temporal autocorrelation of the different wind farms. The marginal distributions of the wind speed stochastic process were also satisfactorily reproduced. Sharma et al. (2013) also used a similar approach to generate wind speed scenarios and concluded that ARMA models preserve the distribution parameters and ACF. Both papers focused on the generation of scenarios without a forecast, so they did not discuss the property of the WFE. The flexibility of the ARMA model is showcased in Chen et al. (2010) where a new type of ARMA model, named LARIMA (limited-ARIMA) is proposed. The ‘Limited’ term in the name LARIMA indicates that the model outputs are restricted by the physical limitations of the wind power generators. In this study, wind speed was considered to be a non-stationary stochastic process, which 28  necessitated the use of a variance stabilizing transformation and differencing. The outputs of the ordinary ARMA and ARIMA models match the mean of the observed data but have overestimated its variance. The authors claim that this deficiency is the result of using unbounded ARMA models to model a bounded stochastic process. To solve this issue, they suggested a LARIMA model, where the autoregressive parameter is limited to a maximum value of 1. Since this limit would remove the nonstationary property of the process, a deterministic trend was used to reintroduce it. Using this approach, the LARIMA parameters that resulted in the best fit to the observed probability distribution were chosen. The authors then compare the output of this model to outputs of a first-order MC model in terms of ACF, partial autocorrelation (PACF), cumulative distribution function (CDF), number of model parameters and the ability to capture seasonal properties. The comparison showed that the LARIMA model is superior in reproducing the statistical properties except for the ACF, where it performed comparably with the MC model. The LARIMA model also has a much lower number of model parameters than the MC model, which suggests that less data is needed to estimate its parameters accurately. The LARIMA model was also found to provide a better fit in CDF than the ordinary ARMA/ARIMA models, especially at the tails of the distribution.  Other types of scenario generation techniques have also been applied successfully. Ma et al. (2013) presented a scenario generation method that can be used when a forecast is available. Realizing that the property of the WFE is dependent on the forecast, the authors categorized the normalized forecast values into 50 different bins and independently analyzed the WFE for each bin. The analysis showed that the WFE distributions vary from one forecast bin to other. The authors then created an empirical cumulative distribution (ECDF) for each forecast bin to approximate the property of the WFE. The ECDF was chosen instead of a CDF to avoid making any assumptions about the distribution. The scenarios were generated by sampling a multivariate normal distribution and transforming the sampling outputs back into wind power using the ECDF for the forecast bin. The structure of the multivariate normal distribution’s covariance matrix was controlled to ensure that the hour to hour variation of the generated wind power scenarios is close to the observed variation. However, the authors did not explicitly model the autoregressive property of the WFE and the ACF of the model outputs were not validated.  29  c. Other Types of Models Moment matching algorithms have also been applied by many researchers for scenario generation. The algorithm proposed by Høyland and Wallace (2001) has especially seen widespread use. This algorithm is based on non-linear optimization and restart-heuristics that minimizes the distance measure between a specified set of statistical moments and the generated scenarios’ statistical moments. Moment matching can guarantee the preservation of the statistical moments, but it does not explicitly model the ACF. Therefore, special attention needs to paid to the model outputs to ensure that the ACF is accurately represented. Lowery and O’Malley (2013) used the technique to create a stochastic tree tool that generates WFE and LFE scenarios. In the proposed technique, the user defines four statistical moments (i.e., mean, variance, skewness, and kurtosis) for each time step and the algorithm creates scenario sets that will have moments that are close to the provided ones. To enable the creation of realistic and representative scenarios, the ACF of the forecast error was also incorporated into the scenario generation process. The authors state that they used the moment matching technique because it enables the user to have greater control of the statistical properties of the outcomes. This is possible in moment matching because the user can control on the moment of the scenario fan that matters the most in the modeling environment. This property has been used to assess the impact of WFE statistics on unit commitment in Lowery and O’Malley (2012) (see Section 2.3).  Greenhall (2013) also used moment matching to generate wind power scenarios. The scenarios were created with a nonlinear optimization algorithm that minimizes the difference between the four observed and generated moments. In the objective function, different weights were assigned to the four moments, with mean and variance having the highest weights. The author found the procedure to have the limitation of creating duplicate scenarios, which might be resolved by varying the weights given to the four moments or by adding more constraints in the optimization. However, even with these limitations, the moment matching algorithm performed better than the other scenario generation technique used, namely the analog method. In the analog method, if for example the DA wind scenarios are needed, days in the past that have similar weather parameters like temperature and air pressure for the current day are mined and the wind power outputs for the following day were taken as the DA scenarios. This approach requires substantial amounts of archived data for the algorithm to perform satisfactorily.   30  Machine learning based algorithms have also been utilized to create wind speed scenarios. Vagropoulos et al. (2016) proposed a wind speed scenario generation technique that is based on artificial neural networks (ANN) and other exogenous variables. The proposed technique uses an ANN that takes in several explanatory variables (like the wind speed from the previous n timesteps and the wind speed forecast for time t) as an input and outputs wind speed prediction for the next timestep. A normally distributed white noise that was derived from the ANN model’s residual was added to the ANN’s prediction to get the realization for that times step. This Gaussian white noise is random, so it ensures that all simulations create unique trajectories. The predicted wind speed from this process at time t was then used as an input for the next time step (t+1), and this rolling update approach was followed until the desired forecast horizon is reached. After generating one scenario for a specified time horizon, the process is repeated starting from the first timestep until the desired number of scenarios are generated. A scenario reduction technique based on probability distance was used to reduce the numerous scenarios generated by the ANN-based method. The authors compared the quality of the scenarios obtained using this method with scenarios obtained from ARMA models and found it to be superior, both in terms of root sum square error and optimization results. The main disadvantage of this approach is the complexity associated with training the ANN, which requires ample experience and good engineering judgment. In this section, several scenario generation algorithms have been reviewed. The literature also shows the prominence of MC and ARMA models for creating wind scenarios. The strengths and weaknesses of MC, ARMA, and moment matching models, as well as the development of innovative scenario generation techniques,  have been discussed. Most of the available literature discusses the generation of short-term wind speed or power scenarios based on the observed statistics without incorporating DA wind forecasts. This indicates that methodologies of scenario generation that appropriately capture the WFE properties are not abundant.     2.4.2 Scenario Reduction Most of the scenario generation methods discussed in Section 2.4.1 can produce as many scenarios as desired. Some of the techniques like the ARMA and Monte Carlo based methods require the generation of numerous scenarios before the property of the underlying probability distribution can be captured confidently. As discussed earlier, the size of the SO problem and its computation burden 31  increases with the number of scenarios considered. This is especially problematic for mixed-integer problems since the size of the problem increases non-linearly with an increase in the number of scenarios.  Therefore, to avoid computationally demanding, or even intractable SO problems, a rational way of reducing the number of generated scenarios is needed.  The main goal of scenario reduction is to reduce the number of available scenarios to a desired fixed cardinality or accuracy without losing much of the relevant stochastic information they contain (Conejo et al. 2010). Several researchers have investigated approaches to achieve this result. Although there is no consistent way of classifying scenario generation methods, Bruninx (2016) classified them as methods based on probability metrics, clustering, importance sampling and moment matching.  Some of the most used scenario reduction methods for electric power related problems are the two algorithms introduced by Dupačová et al. (2003). In the paper, the authors presented two algorithms for scenario reduction that are based on the Kantorovich probability distance. Kantorovich distance (also called Wasserstein distance or Earth mover’s distance) is a metric that is used to compare how close two probability distributions (or sets of scenarios in this case) are. The best way to measure the quality of the scenario reduction would be to compare the optimal solutions from the original and reduced scenario set, but that is impractical because it requires solving a deterministic optimization problem for all of the scenarios. For this reason, the probability distance metric is used as a proxy to measure the quality of the scenario reduction algorithm (Conejo et al. 2010; Dupačová et al. 2003).  The two algorithms proposed by Dupačová et al. (2003) are the forward selection and backward reduction algorithms. In the forward selection algorithm, a single scenario that results in the lowest Kantorovich distance when compared with the original scenario set (i.e., a scenario that is the closest to the original set in terms of Kantorovich distances) is first selected from the original scenario set. This divides the original scenario set into sets of selected and rejected scenarios. The scenario reduction process continues by iteratively selecting one scenario at a time from the set of rejected scenarios and adding it to the set of selected scenarios until the desired accuracy metric or cardinality is reached. The probabilities from the rejected scenarios are then transferred into the selected scenarios. This method works best when a drastic reduction in the number of scenarios is desired or when the desired cardinality or the selected scenario set is small. The backward reduction method 32  works by iteratively deleting a single scenario from the original set that has the highest similarity in terms of Kantorovich distance with another scenario until the desired accuracy is reached or the remaining scenarios match the desired cardinality. In the final stages of both algorithms, probabilities from each of the rejected/deleted scenarios are transferred to the preserved scenarios that are closest to them. Numerical case studies are also presented in the paper which show that both algorithms give satisfactory results. Heitsch and Römisch (2003) presented a faster and more efficient algorithm to perform the calculations suggested in Dupačová et al. (2003). Gröwe-Kuska et al. (2003) applied the algorithm discussed in Dupačová et al. (2003), and Heitsch and Römisch (2003) to solve a stochastic portfolio management problem for a German power utility. Their conclusion was that the optimal solution of the optimization problem could be satisfactorily approximated by using the smaller number of scenarios obtained from this scenario reduction method. It is worth noting that the scenario reduction techniques discussed above are heuristic and that they do not guarantee the most accurate scenario sets. For example, there might be another reduced set of scenarios with the same cardinality that might offer a better representation of the original scenario set.  In Dupačová et al. (2003), the distance between scenario pairs is calculated by the Euclidean norm of the difference between the scenario pairs. Bruninx (2016) and Morales et al. (2009) argue that this is not an appropriate metric to compare scenario sets. They suggest running a deterministic optimization for every single scenario and using the optimal value of each optimization as a comparison metric. They claim this approach captures the effect of the scenarios on the optimal value and enable the selection of the most representative scenarios. This approach is especially beneficial for wind power as the norm of the difference between scenarios will not take into account the variability of wind power in the individual scenarios, which can have a significant impact on the solution of the SO problem. The disadvantage of this approach is the computational burden of having to solve an optimization problem for each scenario.  Sumaili et al. (2011) presented a clustering approach to solve the scenario reduction problem. The algorithm takes numerous wind power scenarios generated from a Monte Carlo method and applies an evolutionary particle swarm optimization to cluster scenarios that are similar. The similarity criterion between two wind power scenarios was defined by the maximum deviation between the 33  two scenarios in a 24-hour period. One focal scenario from the cluster is used as a representative element of the cluster, and a probability that is proportional to the number of elements in the cluster is assigned to it. A test performed on a simplified unit commitment model showed that the method could give a realistic representation of the wind forecast uncertainty. The use of importance sampling as a means for scenario reduction can also be found in the literature. For example, Papavasiliou et al. (2011) and Papavasiliou and Oren (2013) used importance sampling in the context of wind power to select the most appropriate and representative scenarios from a scenario set. In these two publications, the authors claim that the scenario reduction algorithm proposed by Dupačová et al. (2003) and Heitsch and Römisch (2003) perform poorly in real life unit commitment problems without transmission constraint and contingencies. They suggest that this is because the algorithms do not guarantee the preservation of the wind power generation statistical moments. They also add that the algorithms do not enable the user to select and use some specific scenarios that are deemed crucial for the system stability (e.g., scenarios with the minimum wind production or maximum error). The importance sampling algorithm proposed in Papavasiliou et al. (2011) works in two steps. First, the scenarios deemed important for the modeling are selected according to 11 selection criteria from an original set of scenarios generated by an AR model. Then, the selected scenarios are assigned weighted probabilities. The probabilities are assigned in a way that matches the means of the selected and original scenarios. Papavasiliou and Oren (2013) improved the algorithm by fashioning a new procedure that selects the ‘important’ scenarios based on the optimal value of each individual scenario. Four numerical experiments showed that the new procedure has a better performance than other rival algorithms. An added benefit of the new procedure is that it removes the subjectivity of selecting the ‘importance criteria’ that is found in Papavasiliou et al. (2011).   Moment matching, in addition to its use as a scenario generation tool, can also be used as scenario reduction tool. Li et al. (2016) proposed an algorithm that is partially based on moment matching to reduce scenarios of WPTS. The algorithm attempts to simultaneously minimize both probability and moment distances between the reduced and originally sampled scenarios. The optimization algorithm that achieves this reduction is a highly nonlinear, difficult-to-solve problem, which necessitated a heuristic solution approach. The authors proposed a heuristic search algorithm that iteratively and randomly selects scenarios from the original set and then calculates the probability and moment distances until a predefined accuracy threshold is achieved. This scenario reduction method was 34  compared with the backward reduction algorithm and a particle swarm optimization based scenario reduction method. The result showed that the heuristics search algorithm is superior at minimizing the probability and moment distances. The computation time is also shown to be significantly smaller than the backward reduction algorithm, especially when the number of original scenarios is large.   Dvorkin et al. (2014) compared four scenario reduction algorithms: fast forward selection formulated in Dupačová et al. (2003), backward reduction, k-means clustering and importance sampling. The algorithms were compared by running an optimization model for DA electricity market bidding. The outcomes of the comparison showed that the forward selection algorithm has the best accuracy and the least computational burden.  2.5 Stochastic Optimization of Wind Power Optimization models have been applied to solve hydropower management, planning, and operation problems. Reservoir operation and hydropower optimization models, like other optimization models, are defined by their objective functions, constraints, parameters, and variables. Typical objective functions of reservoir operation problems include the minimization of operating cost or the maximization of net trade benefits, energy generation, reliability, sustainability, or other types of benefits. More than one objective can also be considered within the multi-objective optimization framework. Typical constraints in the hydropower optimization field include constraints that impose limits on storage, or reservoir elevations, turbine discharge, power production, spill discharges, transmission capacity and ramps. Other types of constraints that are considered include the continuity equation, load-resource balance, hydraulic equations, and turbine efficiency curves. The common variables are power generation, discharge through turbines, spills, trade schedules and units to be committed. An extensive overview of the various types of algorithms that are used to solve reservoir operation problems is given in Yeh (1985) and Labadie (2004), where multiple algorithms including linear programming, nonlinear programming, dynamic programming, networking flow optimization, stochastic optimization, multi-objective optimization, heuristics and machine learning are discussed.  Yeh (1985) states that when important parameters are uncertain, it is crucial to evaluate the expected performance of the system and its reliability. The evaluation can be achieved by using SO. As discussed in Chapter 1, SO are optimization algorithms where at least some of the parameters’ uncertainty is 35  modeled explicitly, and perfect knowledge of the future is not assumed. SO has been applied to solve reservoir optimization problems since the early 1960s. Most of the parameters that were considered uncertain in the earlier models were inflows, load, market prices, and unit availability. Recently, with their growing prominence, highly uncertain renewable energy resources such as wind are also treated as stochastic parameters.  SO is recommended in an electric system with significant wind penetration because of the high level of uncertainty that is associated with wind power at different timeframes (Conejo et al. 2010; Juan M. Morales et al. 2010). Several researchers have studied the utility of stochastic optimization in reducing the cost of integrating wind power that arises from its variability and uncertainty. Some of the relevant research projects are discussed below. Note that all the projects discussed below are not dealing with DA uncertainty. Some of them address the management of the within-the-day wind power uncertainty, but the methodologies can be extended to manage DA uncertainty. Bathurst et al. (2002) developed a SO model to trade wind power in a short-term (i.e., less than 24 hours) deregulated electricity market while minimizing the imbalance cost. The imbalance costs in this study were the ‘top up’ cost which was paid when wind generation is lower than the contracted amount and the cost of spilling water when wind generation is higher than contracted amount. A persistent forecast was assumed for the wind energy in the short term. A Markov model was developed that divides the wind energy generation into 11 states and creates a state transition probability matrix for different forecast lead times. Then, a bidding strategy that minimizes the expected imbalance cost or maximizes the revenue for the predicted wind generation was devised. The model weighs the imbalance costs against the contract price for all possibilities in a given time delay and determines the optimum bid. This approach reduces the imbalance cost as compared to a bidding strategy based on a deterministic persistence forecast. The reduction in imbalance cost was found to be more significant when the imbalance cost is higher than the contract price.  Matevosyan and Söder (2006) presented a mixed-integer SO model that maximizes the profit (i.e., minimizes the imbalance cost) for a wind farm operator bidding in a DA market. In this study, the wind speed forecast error (WSFE) is represented as a stochastic process using scenarios. Several WSFE scenarios were created using an ARMA(1,1) model and added to the original wind speed forecast to produce several wind speed scenarios. These scenarios were then converted into wind 36  power scenarios using the power curve. The cardinality of these scenarios is then reduced by the scenario reduction method presented in Heitsch and Römisch (2003). This paper assumed a market environment that includes both a DA market and a regulation (balancing) market. The authors assumed the wind to be a price taker, which is an assumption that is common when the wind penetration in the electricity market is low. The objective function included terms for the revenue from selling power and costs incurred for balancing. The optimization was run independently for each hour since the hour to hour coupling of operation is not expected in a wind farm. The bidding strategy was applied with real data obtained from Denmark. The data was used to generate 1000 WSFE scenarios which were then reduced to 50 so that they can be used in the optimization. The results of the optimization showed that depending on the market prices, the profit from the stochastic model is higher or equal to the profit from a deterministic model where the operator bids all the forecasted energy and pays the imbalance cost. For the test case of March 2003, the profit from the stochastic model was found to be 5.45% higher than the deterministic one.  Morales et al. (2010) further developed the approach proposed in Matevosyan and Söder (2006) by improving the mathematical formulation of the optimization problem and by including adjustment markets. They also added risk control into their formulation. The optimization is formulated assuming a pool-based market environment that has DA, adjustment and balance markets. In this market environment, adjustment transactions are carried out closer to power delivery time (usually few hours), and the balancing transactions are carried out few minutes before power delivery. The objective was to maximize the profit from the DA and adjustment markets while minimizing the penalties incurred in the balancing markets. ARMA models were utilized to characterize the uncertainties in the market prices and wind power production. The optimization problem was formulated as a multistage stochastic linear programming problem, with the stages referring to the DA, adjustment and balancing market stages. Risk control was incorporated into the model using the ‘conditional value at risk’ risk measure. Through a case study, the benefits of the adjustment markets, risk control, and a stochastic formulation were demonstrated. The stochastic formulation was shown to have a higher benefit when the volatility of the stochastic processes is higher.  The papers reviewed above discussed offering strategies for a power producer whose product was only wind power. The combined operation and optimization of wind power with hydropower has also been investigated by many researchers. Angarita et al. (2009) presented a SO technique that 37  optimizes the combined hydro-wind operation in a pool-based market. In this paper, the authors argue that the combined bid can mitigate the substantial penalties incurred by wind power producers in the electric market. To test this theory, they presented a mixed-integer SO algorithm with the objective function of maximizing the combined profit of the two plants. Detailed modeling of the hydro operation was undertaken, with the highly nonlinear hydropower-related functions presented as concave piecewise linear functions. The wind power was considered to be a stochastic process, with three power generation scenarios for each hour. The results of the optimization showed that up to 45% and 15% of the imbalance penalties could be saved by the joint bid when the wind forecast mean absolute errors are at about 28% and 4.3% respectively. One of the conclusions of the study was that offering a combined bid is more attractive when the imbalance prices are high.  Tuohy et al. (2008) assessed the benefit of stochastic scheduling of electric systems with significant amounts of wind power. The benefit was assessed by running a stochastic unit commitment and comparing the results with two other approaches: a deterministic unit commitment model with reserves and a deterministic unit commitment model that assumes a perfect forecast (i.e., no reserves). The unit commitment problem was formulated as a multi-stage mixed-integer stochastic optimization problem with hourly rolling planning that minimizes the system cost. The system cost comprises of import/export transaction as well as fuel, start-up, emission, and reserve costs. The model was tested on a future resource portfolio of the Irish power system, which has 6000 MW of installed wind capacity, about 42% penetration.  Numerous scenarios for wind generation, load, unit availability, and reserve were generated with a Monte-Carlo method and were subsequently reduced. The results from an optimization run for a duration of one year showed that stochastic optimization could lower the system cost by 0.6%. The results also showed that perfect forecast would lower the system cost by about 1.7%. This suggests that stochastic unit commitment can reduce the cost of uncertainty in the Irish power system by more than one-third. Ruiz et al. (2009) investigated the benefit of a unit commitment approach that uses SO and reserves to manage the DA uncertainty of wind power. The traditional approach to deal with the DA wind uncertainty has been to run a deterministic unit commitment while keeping operating reserves to ensure reliability. The other approach has been the use of SO without explicitly committing operating reserves. In this study, the authors argue that including the reserve constraint in a stochastic formulation results in better policies that lower the chance of wind power curtailments. The 38  optimization problem was formulated as a mixed-integer problem with an objective function that minimizes the system cost while satisfying the constraints. One of the constraints is the required reserve amount, which is made unique for each scenario. The parameters that were considered uncertain were system load, wind generation, and unit availability. The performance of the model was examined by implementing it on the Public Service Company of Colorado (PSCo) system. The stochastic solution was compared with the traditional modeling approach and a perfect foresight benchmark policy. The traditional policy (i.e., deterministic optimization with reserves) has the highest system cost, which was lowered by 1.7% when perfect foresight (i.e., deterministic, no reserves) was assumed. The stochastic policy had a saving of 0.2% over the traditional policy, which is lower than the saving reported by Tuohy et al. (2008) and some other researchers. The authors state that this lower stochastic saving is probably the result of the flexibility of the PSCo system, which has a large number of fast start turbines and a pumped-hydro unit. An additional benefit of the stochastic policy is that it reduced the duration of wind power curtailment. Greenhall (2013) used two scenario generation methods, the moment matching technique and an ‘analogs’ technique (see Section 2.4.1) to study the effect of stochastic optimization in the Electric Reliability Council of Texas (ERCOT) system. Scenarios created from both techniques were implemented in a stochastic optimization problem for the ERCOT system. The results indicated that around 0.2% cost saving could be obtained from the stochastic formulation when the wind penetration is around 25%. This result is considered ‘good' since the maximum saving that can be expected from a perfect forecast in the ERCOT system at 25% penetration is 0.85%. The author notes that although the percentages appear small, they are significant in dollar value ($6 million in 9 months). Surprisingly, the stochastic savings reduced to 0.06% for both scenario generation methods when the penetration was set at 30%. The reason given for this phenomenon was that the reduced net load at higher wind penetration provides less opportunity for cost savings since surplus capacity is available. The moment matching technique was shown to provide better cost reductions than the analogs technique. Stochastic savings were also seen to increase with the number of scenarios included in the optimization.   Zima-Bočkarjova et al. (2010) presented a stochastic optimization-based algorithm for the coordinated operation and DA bidding of wind and hydropower. A novel approach that is based on game theory was also proposed to enable fair and transparent sharing of the extra profit gained from 39  the coordinated operation. The model was formulated as a two-stage SO with recourse. The wind power and the market prices were modeled as stochastic parameters with three scenarios each. A case study that was carried out to assess the performance of the algorithm showed that coordinated operation increases the profit for both wind power producers and the hydropower producers. The effectiveness of the profit sharing scheme was also demonstrated.  Garcia-Gonzalez et al. (2008) discussed the joint operation and bidding of pumped-hydro and wind power for an electricity market. The optimization was formulated as a two-stage mixed-integer optimization problem, with the DA bids as the first stage decisions and the spot market transactions as the recourse (second stage) decisions. The wind energy and market prices were considered as stochastic processes. Rivas Guzman (2010) applied a similar model to assess the value of a pumped-storage hydropower plant for integrating wind to the BC Hydro system. In the model, load scenarios were created using the moment matching technique, and wind scenarios were generated with a Markov model. Both wind and load scenarios were reduced based on the method discussed in Dupačová et al. (2003).  PandŽić et al. (2013) similarly presented a stochastic offering model for a virtual power plant that has the generating resources of wind, pumped-hydro, and combined cycle gas turbines. The focus of these studies was not to quantify the benefit of stochastic optimization, so stochastic savings were not reported. However, their adoption of stochastic models indicates that the value of stochastic optimization has been recognized by researchers.  In general, this section of the literature review shows that SO is a widely used approach to deal with the electric system. The publications also show that the stochastic savings vary significantly from one electric system to another and that they depend on several factors such as types of generators available, the amount of uncertainty in the system, the market structure, the geographic distribution of wind farms, and season of the year.  40   Modeling Methodology In this chapter, the modeling methodology that is applied to achieve the goals of the thesis is detailed. The main goals of the thesis will be discussed in detail under different sections. The first section describes the approaches used to statistically analyze the WFE data. The second and third sections lay out the methodology used to generate scenarios using Markov chain and ARMA models. The fourth section explains the process of validation and evaluation of the scenarios generated using the two models. The fifth section outlines the scenario reduction algorithm used. The sixth and seventh sections outline the details of the simulation and optimization model used in this thesis.  The eighth section explains how the results of the optimization are analyzed.  3.1 Statistical Analysis of the Day Ahead Wind Forecast Error One of the challenges of wind power is its uncertainty. As mentioned in the first two chapters, the WFE in the DA timeframe is one of the primary sources of the wind integration cost. Therefore, identifying the relevant statistical properties of the WFE is a crucial step in the development of mitigation strategies. In this thesis, the WFE value for a particular period t is defined as:  𝑊𝐹𝐸𝑡 =  𝐷𝐴 𝑊𝑖𝑛𝑑𝐺𝑒𝑛𝐹𝑟𝑐𝑡 − 𝐴𝑐𝑡𝑊𝑖𝑛𝑑𝐺𝑒𝑛𝑡   𝑀𝑎𝑥𝑖𝑚𝑢𝑚 𝑊𝑖𝑛𝑑 𝐺𝑒𝑛 1 where 𝐷𝐴 𝑊𝑖𝑛𝑑𝐺𝑒𝑛𝐹𝑟𝑐𝑡 is the DA wind generation forecast in MW, 𝐴𝑐𝑡𝑊𝑖𝑛𝑑𝐺𝑒𝑛𝑡 is the actual wind generation (realization) in MW and 𝑀𝑎𝑥𝑖𝑚𝑢𝑚 𝑊𝑖𝑛𝑑 𝐺𝑒𝑛 is the maximum possible wind generation (i.e. the installed capacity) in MW. The forecast error is normalized by the 𝑀𝑎𝑥𝑖𝑚𝑢𝑚 𝑊𝑖𝑛𝑑 𝐺𝑒𝑛 so that the WFE value can only vary from -1 to +1. In this definition, negative error represents underforecasting of wind generation, and positive error represents overforecasting of wind generation. All subsequent analysis is done on these normalized values. After defining the WFE, the next step is to identify the relevant statistical properties, namely, the statistical moments, seasonality, stationarity, probability distributions, autocorrelation, and partial autocorrelation.  3.1.1 Statistical Moments  The statistical moments that are considered in this study are the mean WFE, mean absolute WFE (MAE), variance, and kurtosis. These statistical moments are defined as: 41   𝑀𝑒𝑎𝑛 = ?̅? =∑ 𝑊𝐹𝐸𝑡𝑛𝑡=1𝑛 2   𝑀𝑒𝑎𝑛 𝐴𝑏𝑠𝑜𝑙𝑢𝑡𝑒 𝑊𝐹𝐸 = 𝑀𝐴𝐸 =  ∑ |𝑊𝐹𝐸𝑡|𝑛𝑡=1𝑛 3  𝑉𝑎𝑟𝑖𝑎𝑛𝑐𝑒 = (𝑆𝑡𝑎𝑛𝑑𝑎𝑟𝑑 𝐷𝑒𝑣𝑖𝑎𝑡𝑖𝑜𝑛)2 = 𝑠𝑡𝑑2 =∑ (𝑊𝐹𝐸𝑡 − ?̅?)2𝑛𝑡=1𝑛 − 1 4     𝐾𝑢𝑟𝑡𝑜𝑠𝑖𝑠 = 𝐾 =  𝑛(𝑛 + 1)(𝑛 − 1)(𝑛 − 2)(𝑛 − 3)∗  ∑ (𝑊𝐹𝐸𝑡 − ?̅?)4𝑛𝑡=1  𝑠𝑡𝑑2   5 where t is the time step and n is the number of total data points. Since the WFE takes both positive and negative value that cancel each other out when added, the mean only shows if there is bias in the forecast. The MAE, which is calculated by averaging the absolute value of the WFE as given in Equation 3, shows the magnitude of the average forecast error irrespective of the error sign. The standard deviation (std, σ), which is the square root of the variance given in Equation 4, shows the variation in the WFE. Kurtosis (K) is a measure of the combined weight of the tails relative to the rest of the distribution and is sometimes used to measure the ‘peakedness’ of a distribution. It indicates whether most of the variation is caused by few large errors (i.e., outliers) or many small errors. As mentioned in Hodge et al. (2011) and Hodge and Milligan (2011), higher kurtosis suggests better forecast accuracy. In general, high kurtosis values (i.e., greater than the K = 3 for the normal distribution) suggest a ‘sharper’ distribution where most of the variance is caused by few large errors. The other commonly used moment, the skewness, was not included in this analysis because it did not provide a lot of additional information (i.e., asymmetry was not a major feature of the WFE time series).  3.1.2 Stationarity and Seasonality Non-stationarity is the condition where the statistical parameters of a time series like the mean, variance and ACF change with time. Seasonality, which is a special case of non-stationarity, is a condition in a time series where the data shows a predictable and reoccurring pattern in a given period of time. For example, demand for electric energy, wholesale energy market prices, and reservoir inflows show seasonal trends. Non-stationarity and seasonality affect the model identification and development process; thus, they have to be investigated.  (Box et al. 2008; Conejo et al. 2010). In this thesis, stationarity is used to refer to the behavior of the WFE with respect to the forecast lead time (i.e., the time between the forecast issuance and the realization). Forecast accuracy 42  is considered to progressively decrease as the forecast lead time increases and as a result, many researchers like Ma et al. (2013) and Pinson et al. (2009) consider the WFE to be influenced by the forecast lead time.  In this study, the stationary of the WFE with regards to the forecast lead time is investigated using bar graphs, trend analysis, and ACF statistics. A similar approach is also used to examine the seasonality. The ‘auto.arima’ tool in the ‘forecast’  statistical package (Hyndman 2017) in the ‘R’ programming language is also used to assess if the assumption of stationary and seasonality improves model performance. 3.1.3 Probability Distributions The joint, marginal and conditional probability distribution functions (PDFs) of the WFE are important statistical properties that display the ranges of possible WFE outcomes and their probabilities. The identification of correct distributions to represent the WFE is vital since the misrepresentation of the characteristics of the WFE can lead to reduced reliability and efficiency of system operations (B.-M. Hodge et al. 2012; Hodge & Milligan 2011; Archer et al. 2017).  The joint PDF considered in this thesis is a bivariate function of the WFE and an exogenous variable. The exogenous variable considered in this thesis is the wind forecast. Many researchers like Bludszuweit et al. (2008), Bruninx, (2016), Ma et al., (2013), and Tewari et al., (2011) have indicated that the WFE distribution and property is influenced by the wind forecast level. Therefore, this association is also investigated in this thesis. The joint PDF is one of the tools used to investigate the influence of the forecast on the WFE. To define the joint PDF formally, let X be the WFE random variable and Y to be the wind forecast random variable. For these two continuous random variables, the joint probability distribution fXY(x,y) is defined such that  𝑃[(𝑋, 𝑌) ∈ 𝐴] =  ∬ 𝑓𝑋𝑌(𝑥, 𝑦)𝑑𝑥𝑑𝑦 𝐴    𝑎𝑛𝑑   ∬ 𝑓𝑋𝑌(𝑥, 𝑦)𝑑𝑥𝑑𝑦 = 1∞−∞ 6 where (X,Y) є A is an event in the xy-plane. The marginal PDF fX(x) shows the univariate WFE distribution without consideration for the forecast. It is given by: 43   𝑓𝑋(𝑥) =  ∫ 𝑓𝑋𝑌(𝑥, 𝑦)𝑑𝑦∞−∞  7 The conditional PDF of X given Y, 𝑓𝑋|𝑌(𝑥|𝑦), shows the WFE distribution conditioned upon the wind forecast. It is given by:  𝑓𝑋|𝑌(𝑥|𝑦) =𝑓𝑋𝑌(𝑥, 𝑦)𝑓𝑌(𝑥) 8 If the WFE distributions that are conditioned on the wind forecast are dissimilar, it is appropriate to use the wind forecast as an exogenous variable. If the conditional distributions are not very different from each other or from the marginal distribution, then the marginal distribution is adequate for the modeling process and an exogenous variable is not needed.  The empirical joint, marginal and conditional distributions of the WFE are first studied using histograms and scatter plots. Then, distribution fitting is carried out to determine if any parametric distributions offer an acceptable fit to the given data. The distribution fitting is done using the ‘distribution fitting’ application on MATLAB(2017a) together with allfitdist MATLAB function (Sheppard 2012) which is obtained from the MathWorks website. The distribution fitting application allows the user to test several distributions manually while the allfitdist function automatically fits the appropriate distribution to the data and lists the ones that offer better fits. Both approaches use the maximum likelihood estimation to estimate the distribution parameters and fit the models. If suitable parametric distributions that accurately capture the properties of the WFE are not found, then empirical distributions will be adopted.  3.1.4 Autocorrelation and Partial Autocorrelation Functions As mentioned in Section 2.4.1, the ACF illustrates the similarity between a time series and a time-lagged version of itself. It is a metric that indicates the persistence level of a time series. The ACF contains information about the autocorrelation coefficient (ACC) at different time lags. The ACC is a metric that measures the linear correlation of a time series among observations made at specific lags (lag-k, with k being 1hr, 2hr, etc.). The lag-k ACC is estimated by shifting the time series data k steps and calculating the linear correlation with the unshifted time series. More formally, for measurements, 𝑥1, 𝑥2, ..., 𝑥𝑡at time t1, t2, ..., tn and mean of 𝜇, the lag-k ACC 𝜌𝑘 is defined as: 44   𝜌𝑘 =𝐸[(𝑥𝑡 − 𝜇)(𝑥𝑡+𝑘 − 𝜇)]√(𝐸[(𝑥𝑡 − 𝜇)2]𝐸[(𝑥𝑡+𝑘 − 𝜇)2])=  𝐸[(𝑥𝑡 − 𝜇)(𝑥𝑡+𝑘 − 𝜇)]𝜎𝑥2 9 ACCs can take values between -1 and 1. High ACC values (i.e., close to positive 1) indicate a strong persistence property in the time series while low ACC values (i.e., close to -1) indicate high-frequency oscillations in the time series. ACC values close to zero suggest that there is little or no persistence in the time series. The ACF plot, which shows these ACC values at various time lags, is often plotted with a confidence interval for the ACC estimates. If the ACC lies outside of the no-autocorrelation confidence interval which is centered at zero ACC, the time series is assumed to have autocorrelation (Chen et al. 2010; Box et al. 2008).  The partial autocorrelation coefficient (PACC) is the measure of the correlation of the realizations at time t and time t+k when the effects of the realizations at t+1, t+2,.., t+k-1 are removed.  For example, the lag-2 PACC shows the correlation between the realization at time t and t+2 after removing the effect of the realization at t+1 on t+2.  The k-lag PACC for an ith order AR process, 𝜙𝑘𝑖, can given by:  𝜙𝑘,𝑖 =𝜌𝑘 − ∑ 𝜙𝑘−1,𝑖 𝜌𝑘−1,𝑖𝑘−1𝑖=11 − ∑ 𝜙𝑘−1,𝑖 𝜌,𝑖𝑘−1𝑖=1 10 where 𝜌𝑘 is the lag-k ACC, 𝜙0,0 is 0 and  𝜙𝑘−1,𝑖 =  𝜙𝑘−2,𝑖 − 𝜙𝑘−1𝜙𝑘−2,𝑘−𝑖. The partial autocorrelation function (PACF), like the ACF, contains information about the PACC at different time lags. Its plot is obtained by plotting the PACCs of various time lags along with the confidence bands for their estimates. The PACF usually determines the order of an MC model and the p parameters of an ARMA(p,q) model (Chen et al. 2010; Box et al. 2008). 3.2 A Markov Chain Based Scenario Generation Model Once the relevant properties of the WFE are identified, a scenario generation model that can create scenarios that capture the identified properties is developed. The first model proposed in this thesis is based on MCs. As mentioned in the literature review, MCs are stochastic processes that can be characterized by the transition probabilities between discretized system states (Shamshad et al. 2005). The MC scenario generation model proposed in this thesis is a higher-order conditional MC model. 45  The order of the MC model, which is the number of previous timesteps that influence the current timestep (as mentioned in Section 2.3.1), is decided on based on the PACF of the WFE time series. An exogenous variable (i.e., wind forecast) is used to condition the transition probabilities. A Markov chain modeling process starts with the selection of the appropriate number of states. There are no standard guidelines on how to select the number of states or their upper and lower bounds. This decision is influenced by many factors such as the input data size and accuracy desired. The decision usually requires a compromise between model accuracy and model simplicity. A large state number improves the accuracy of model outputs but requires more input data. A large number of states also means that the number of model parameters is higher, which is undesirable because the parameters become too dependent on the input data (i.e., a slight change in the input data leads to a significant change in the parameter values). The estimation of a large number of parameters also introduces errors and widens the confidence intervals around the parameters. Therefore, the principle of parsimony is applied in the model selection. Under this principle, the purpose of model selection is to obtain a model with the smallest number of parameters that will provide satisfactory results (Box et al. 2008; Chen et al. 2010; Papaefthymiou & Klockl 2008).  The method proposed in this thesis for selecting the appropriate number of states is to discretize the WFE into different number of states and compare the ACF and PACF of the discretized and continuous WFEs. The smallest number of states that provides ‘good enough’ result is then selected. Once the number of states is decided, the next step is the discretization of the WFE into different mutually exclusive states. The step is described following the notation used by Tang et al. (2015). For k number of states, the range containing the minimum and maximum WFE (εmin, εmax) is divided into k disjoint intervals that define a state, with the ith state holding all the WFE values that fall between the intervals 𝜉𝑙,𝑖  (lower bound) and 𝜉𝑢𝑖  (upper bound). More formally, the intervals [𝜉𝑙,𝑖 𝜉𝑢𝑖 ] for i,j=1,..,k are created such that:   [𝜉𝑙,𝑖 𝜉𝑢𝑖 ]  ∩  [𝜉𝑙,𝑗𝜉𝑢𝑗] =  ∅ 𝑓𝑜𝑟 𝑖 ≠ 𝑗 and  [𝜀𝑚𝑖𝑛, 𝜀𝑚𝑎𝑥]  ⊂ ⋃ [𝜉𝑙,𝑖 𝜉𝑢𝑖 ] 𝑘𝑖=1  11 Then, the WFE values are assigned WFE states using the following rule: 46   𝐸𝑡 = 𝑖 𝑖𝑓 𝜉𝑙,𝑖 ≤ 𝑊𝐹𝐸𝑡 < 𝜉𝑢,𝑖  12 where 𝐸𝑡  is the WFE state at time t. By following this procedure, k WFE states are created. A similar method is used to divide the wind forecast into various states.  After WFE is discretized, the next step is to calculate the state transition probabilities. The state transition probability 𝑃𝑗𝑘𝑓,𝑖 is defined as the probability that the state will transition to state i given the WFE states j,..,k from the previous t-1,.., t-n time steps and the state of the wind forecast f at time step t. For S number of WFE states and F number of wind forecast states, 𝑃𝑗𝑘𝑓,𝑖 can be mathematically expressed as:  𝑃𝑗𝑘𝑓,𝑖 = 𝑃(𝐸𝑡 = 𝑖|𝐸𝑡−1 = 𝑗, . . , 𝐸𝑡−𝑛 = 𝑘, 𝜙𝑡 = 𝑓) 13 where Φ is an exogenous variable, i, j, k = {1, 2,.., S} are the WFE states, f = {1, 2,..., F} is the wind forecast state, and n is the order of the MC model. The transition probability 𝑃𝑗𝑘𝑓,𝑖 is estimated by empirically counting the number of transitions that occur to a specific state i given the n+1 conditions (Ε t-1 =j,.. Ε t-n =k, Φ t =f) and dividing it by the total number of transitions to all states i = 1,2….S given the same n+1 conditions. This means the transition probabilities are calculated using maximum likelihood estimation (Shamshad et al. 2005; Tang et al. 2015). Formally, if N is the frequency of the transitions made to state i given the three conditions and S is the total number of states, the transition probability 𝑃𝑗𝑘𝑓,𝑖is estimated as:   𝑃𝑗𝑘𝑓,𝑖 =  𝑁𝑗𝑘𝑓,𝑖∑ 𝑁𝑗𝑘𝑓,𝑖S𝑖=1 14 For sampling purposes, this matrix is first converted to a cumulative state transition probability matrix (CSTPM) 𝐶 by:   𝐶𝑗𝑘𝑓,i =  ∑ 𝑃𝑗𝑘𝑓,m𝑖m=1 15 where j, k, i = {1, 2, …S} are the WFE states and f = {1, 2, …F} is the wind forecast state. 47  The WFE time series is then generated from the CSTPM using an inverse transform method. A pseudo-random number generator that gives a uniform distribution between 0 and 1 (U ~ Unif[0,1] ) is used to sample the required number of points and the CSTPM 𝐶𝑗𝑘𝑓,𝑖 is used to convert the sampled points to WFE. For t = {3 ..T}, where T is the length of simulation, the WFE state at time t (Et) is created from the random number at time t (U[t]) using:  𝐸𝑡 = 𝑖 if   𝐶𝑗𝑘𝑓,𝑖−1  < 𝑈[𝑡] ≤  𝐶𝑗𝑘𝑓,𝑖 16 The first two hours of the simulation do not have the CSTPM 𝐶𝑗𝑘𝑓,𝑖 since they do not have the error states at time t and t-1. Therefore, cumulative transition probability matrices based on the matrices Pfi and Pjfi are used for those two hours. The outputs of the simulation are WFE states, which are then converted to WFE values by replacing them with representative values for the states.  This simulation approach can be used to sample as many WFE scenarios as desired. The number of scenarios required to confidently capture the underlying WFE distribution is determined by creating numerous scenarios for a single day and checking how many scenarios it takes before the mean and standard deviation start to stabilize. By the law of large numbers (LLN), the sample statistical parameters are expected to approach their population counterparts as the number of samples increase. When the statistical parameters stop varying as new scenarios are added, it can be safely assumed that adequate points have been sampled to accurately capture the underlying distribution. This investigation is done on various days with varying wind forecast levels so that the number of scenarios selected is adequate to capture the WFE distribution regardless of the wind forecast or season of the year.  3.3 An ARMA Based Scenario Generation Model The second approach that was investigated to model the WFE is the application of an autoregressive moving average (ARMA) model. ARMA models are a combination of autoregressive (AR) and moving average (MA) models that have been used extensively to model autocorrelated stochastic processes. There are several variations of ARMA models that are well suited for the modeling of some specific types of stochastic time series such as ARIMA for non-stationary process, SARIMA for seasonal time series, MARIMA for multivariate analysis, ARIMAX for analysis using an external regressor, etc. An 48  ARMA(p,q) stochastic process X with p number of autoregressive parameters and q number of moving average parameters can mathematically be expressed as:  𝑋𝑡 = 𝑐 +  ∑ 𝜙𝑖𝑋𝑡−𝑖 + ∑ 𝜃𝑖𝜀𝑡−𝑖𝑞𝑖=1+ 𝜀𝑡𝑝𝑖=1 17 where 𝑐 is the constant term; 𝜙1, 𝜙2… 𝜙𝑝 are the autoregressive parameters; 𝜃1, 𝜃2 ... 𝜃𝑞  are the moving average parameters, and εt is an error term with zero mean and constant variance (i.e., white noise). The different combination of 𝑐, 𝜙, 𝜃, and 𝜀  are able to represent most types of stochastic time series. For a time series being modeled by an ARMA model, the 𝜙 and θ parameters define its ACF and PACF, the c term is related to its mean, and variance of εt is related to its variance (Box et al. 2008).  The initial estimates of the ARMA parameters p and q are determined from the ACF and PACF. Different reasonable combinations of p and q parameters are then used to fit an ARMA model using the R programming language. The ‘arima’ function in R fits the model through maximum likelihood estimation (R Core Team 2017). After fitting the models, the Akaike Information Criteria (AIC) values that were associated with the models is used to compare the quality of the models. AIC compares various models based on their ‘goodness of fit’ and number of model parameters. AIC penalizes the number of parameters and therefore promotes more parsimonious models (Box et al. 2008). The selection is further verified by using the ‘auto.arima’ function in the ‘forecast’ package for the programming language R (Hyndman 2017). The function selects the most suitable number of parameters and provides estimates for those parameters. Once the parameters are estimated, the ‘arma.sim’ function in the R programming language is used to simulate from the ARMA model.  ARMA processes assume that the input data is a normally distributed stochastic process. The outputs of ARMA model simulations are also normally distributed. Therefore, modeling non-normally distributed stochastic processes without violating the normality assumption of ARMA models requires probability transformation (Conejo et al. 2010). In this study, the WFE’s marginal distribution is transformed to a new stochastic process Z with a marginal standard normal distribution using the cumulative distribution function (CDF) according to: 49   𝑍 =  𝜙−1[𝐹𝑌(𝑌)] 18 where 𝐹𝑌 is the marginal CDF of the WFE and 𝜙 is the CDF of the standard normal random variable. The process is shown graphically in Figure 6 (the graph is only for illustration purposes, the CDF of the WFE in the graph is not the actual CDF). Following the arrows labeled ‘d’ indicates the transformation of the WFE to the standard normal distribution. In this process, the CDF value for the WFE is first read and then the same value from a standard normal distribution for the corresponding CDF value is read. For example, assume that the yi in Figure 6 is a WFE value of 0.1 and its CDF is 0.82. Then, its equivalent in the standard normal distribution is obtained by reading the value that corresponds to CDF of 0.82 in the standard normal CDF. This will give a value of around 0.5. This means the WFE value of 0.1 is transformed to its standard normal equivalent value of 0.5. This is done for every WFE value to get the new transformed time series Z.   Figure 6 - Probability transformation  (based on Conejo et al. (2010)) The output of an ARMA model is a normally distributed time series that is not affected by any exogenous variable. To create representative WFE scenarios that are conditioned on an exogenous variable (i.e., the wind forecast in this case), the model output is reverse transformed using CDFs that are conditioned on the forecast. These conditional ECDFs are derived first by dividing the forecast into various disjoint forecast bins and then creating a CDF for the WFE values that reside in each bin. 50  These ECDFs approximate the joint probability distribution of the wind forecast and the WFE. After obtaining CDFs, the ARMA output at time t is reverse transformed back into WFE using the CDF for the forecast bin at time t. The transformation procedure is shown in Figure 6 by the arrows labeled ‘i’. This inverse transformation can be written as:  𝑌 =  𝐹𝑌𝑓−1[𝜙(𝑍)] 19 where 𝐹𝑌𝑓−1 is the conditional CDF of the WFE for forecast bin f. This method was selected to represent the marginal distribution because testing showed that it performs better at capturing the statistical properties of the WFE as compared to a bivariate ARMA model and an ARMAX model that uses the wind forecast as an external regressor.  Like the MC-based scenario generation method, this ARMA-based approach can also be used to sample as many WFE paths as desired. The same LLN-based methodology used on the MC model to determine the required number of scenarios is also applied for the ARMA model.  3.4 Evaluation of Scenario Quality Once the scenarios are generated using the MC and ARMA based models, the validity of the generated scenarios is investigated. The validation is done by creating a single WFE time series from both models that is equal in length to the original WFE time series (i.e., three years of hourly values) and comparing the statistical properties of the synthetic time series to the properties of the original time series. The statistical properties that are compared this way include the statistical moments (specifically the mean, MAE, standard deviation, and kurtosis), the probability distributions, the ACF, and the PACF. It was not possible to conduct a split sample validation because all the available WFE data is used to develop the two scenario generation models, both of which required large input dataset. The disadvantage of using all the available data for model development without cross-validation is that it might lead to an undetected overfitting problem.   The stability of the scenario generation techniques (i.e., their ability to create scenarios with similar properties when the simulations are rerun) is also examined. The stability is first assessed by simulating the three-year-long synthetic time series ten times and verifying that the statistical properties of all simulations are similar. The analysis is then conducted for a shorter time horizon by 51  creating several scenarios sets for a 24-hr period and comparing their statistical properties. Different scenario sets that are generated from the same stochastic model for the same time interval should have similar statistical properties. Otherwise, each scenario set can result in significantly different unit commitments, which leads to in-sample instability (Rahimiyan et al. 2011; Xu et al. 2015). The similarity of statistical parameters calculated for the different scenarios sets also ensures the underlying distribution is satisfactorily captured.  3.5 Scenario Reduction Approach Once an adequate number of scenarios are generated, and a satisfactory representation of the underlying WFE distribution is guaranteed, an algorithm is devised to reduce the cardinality of these scenarios. As mentioned in Section 2.4.2, scenario reduction is needed because large cardinality of scenarios unnecessarily increases the size and computational burden of the SO problem, potentially rendering it intractable. The scenario reduction methodology applied in this thesis is an approach based on the Kantorovich distance as presented Dupačová et al. (2003). This approach has a proven performance and has found widespread use in the SO problems as discussed in Chapter 2. The algorithm chosen for the scenario reduction is the improved, fast forward scenario selection algorithm presented by Heitsch and Römisch (2003). This heuristic algorithm has a low computational burden and works best when the cardinality of the preserved scenarios is significantly less than the cardinality of the original scenarios.  The scenario reduction methodology used is presented below following the notation found in Morales et al. (2009) and Conejo et al. (2010).  For the distribution 𝑄 defined by the scenario set Ω ,  the goal of scenario reduction is to obtain a reduced scenario subset Ω𝑠 ⊂ Ω  of a desired cardinality and to assign new probabilities to the reduced scenarios so that the reduced distribution 𝑄′ defined by the subset Ω𝑠 provides the best approximation to the original distribution 𝑄 in terms of a probability distance. In this thesis, the probability distance metric used is the Kantorovich distance 𝐷𝑘(. ), which can be defined for the two probability distributions 𝑄 and 𝑄′ as: 52   𝐷𝑘(𝑄, 𝑄′) =  ∑ 𝜋𝜔  min𝜔′∈ Ωs𝜈(𝜔, 𝜔′)𝜔 ∈ Ω\Ω𝑠  20 where πω and τω stand for the probabilities of scenarios ω and ω’ in scenario sets Ω and ΩS respectively. The above problem is known as the Monge-Kantorovich mass transportation problem. 𝜈(𝜔, 𝜔′) is usually called the cost function, and it is used to quantify the closeness or similarity of two scenarios. Given these definitions, the heuristics algorithm used to find the optimal reduced set is started by computing the cost function. The cost function ν(ω, ω′)  is computed for each pair of scenarios ω and ω′ in Ω. The cost function can be calculated in many ways as mentioned in Section 2.4.2. It is calculated here as the Euclidean norm of the difference between pairs of scenarios. More formally, for a stochastic process 𝝀 with initial sets of scenarios 𝜆(𝜔) such that 𝝀 = {𝝀(𝜔) , 𝜔 = 1,2, . . 𝑁Ω}, the cost function is calculated as:  𝜈(ω, ω′) = ‖𝝀(𝜔) − 𝝀 (𝜔′)‖ =  √∑((𝝀(𝜔) − 𝝀 (𝜔′))2 21 where || . || represents the norm operator. This cost function formulation compares the similarities between scenarios using the Euclidean distance between them. The cost function suggested by Bruninx (2016) and Morales et al. (2009), which requires an optimization problem to be solved for each scenario, is not feasible for this study because of the size of optimization problem that is involved, which is large even in its deterministic form. After the cost function is calculated, the forward selection procedure that is followed to produce a reduced scenario set with N number of scenarios can be summarized as: Step 1 – The Kantorovich distance is computed for each scenario, and the scenario that has the lowest Kantorovich distance is selected. This divides the original scenario set into sets of selected and rejected scenarios. Step 2…N – The cost function is updated and the Kantorovich distance are recomputed for the rejected scenarios. From the sets of rejected scenarios, the scenario that has the lowest Kantorovich 53  distance is selected and added to the set of selected scenarios. The process is repeated until the desired number of scenarios are selected.  Step N+1- Finally, after the last scenario is selected, the probability of each rejected scenario is transferred to the selected scenario that is closest to it in terms of Kantorovich distance. The algorithm can be found in detail in Conejo et al. (2010). The compromise between the number of reduced scenarios and information contained in the scenarios is assessed using the relative distance metric that was used in Dupačová et al. (2003), Heitsch and Römisch (2003), and Gröwe-Kuska et al. (2003). The relative distance drel(n) for a reduced scenario set containing n number of elements is estimated by:  𝑑𝑟𝑒𝑙(𝑛) =  𝑑(𝑛)𝑑(a) 22 where d(n) is the Kantorovich distance between original set of scenarios and n number of reduced scenarios while d(a) is the distance between the original set of scenarios and the single optimal scenario that minimizes the Kantorovich distance (i.e., the best deterministic approximation of the scenario set or the first scenario selected in the fast forward algorithm). A drel(n) value of 0 is obtained when d(n) is also 0, which happens when n is the same as the cardinality of the original scenario set (i.e., when no scenarios are removed from the original set). A drel(n) value of 1 indicates that n is 1 (i.e., only one scenario is preserved from the original set).  Thus, 1- drel(n)  reflects the amount of information retained in terms of Kantorovich distance in a range of 0 to 1 (Follestad et al. 2011; Dupačová et al. 2003). A curve showing the number of retained scenario as the independent variable and the 1- drel(n) value as the dependent variable is used to decide the minimum number of scenarios to retain.   A simplified flowchart showing the process of model development and scenario generation is given in Figure 7. In the figure, the boxes with solid lines show some of the steps followed in the model development process and the boxes with broken lines show some of the steps followed to get the reduced wind generation scenarios.  54   Figure 7 - A simplified flowchart showing the process of scenario generation Yes ARMA-based model MC-based model    Statistically analyze the WFE data Decide the number of states to use  Discretize the WFE into WFE states  Estimate state transition probabilities  Create a cumulative state transition probability matrix Sample numerous scenarios for a single day and for three years  Give WFE values for the sampled WFE states Transform the WFE distribution to the standard normal using the marginal ECDF  Decide the proper p and q parameters to use  Fit an ARMA model   Reverse transform the model outputs to WFE values using conditional ECDFs  Reduce the number of scenarios  Add the reduced WFE scenarios to the DA wind forecast to get the DA wind generation scenarios  Are the model outputs acceptable? Check statistical moments, marginal and joint probability distributions, the ACF, and the PACF of the simulation outputs Check stability of the simulation outputs  No No Sample numerous scenarios for a single day and for three years  Use the MC-based model to generate several WFE scenarios for a single day   Use the ARMA-based model to generate several WFE scenarios for a single day   55  3.6 Optimization Models Once the WFE scenarios are generated and reduced, the preserved scenarios are converted into DA wind generation scenarios by adding them to the single DA wind generation forecast. These DA wind generation scenarios are used in a SO model. The final SO model used in this thesis is based on an existing deterministic model at BC Hydro called the generalized optimization model (GOM). GOM uses many inputs from another deterministic modeling tool called the hydraulic simulation model (HYSIM). These two models are discussed below. 3.6.1 Hydraulic Simulation Model (HYSIM) HYSIM is a deterministic simulation tool with a monthly time-step that is used to model the BC Hydro system. HYSIM determines the most economical dispatch of the resource mix under a range of historical streamflow sequences. The resource mix may include hydropower plants, thermal plants, wind farms, import/exports as well as firm energy contracts. HYSIM dispatches the large hydropower plants and makes use of energy import/export to meet residual load (i.e., total load – fixed energy resources) while minimizing operation cost in long-term operation (BC Hydro 2013). The primary purpose of HYSIM is for the estimation of energy capabilities of proposed hydro developments, the estimation of future operating costs, and preparation of inputs for more detailed and granular models like GOM.  HYSIM uses the marginal value of storage ($/MWh) tables of individual hydropower plants along with electricity market prices to make decisions on power plant operations and energy market transactions. The marginal value of storage is generally negatively correlated with the amount of water stored in the reservoir (i.e., marginal value is high when storage is low and vice versa). Plant operation and market transactions decisions are made for each month of the year serially without the knowledge of future information.  In HYSIM, a study for one load year is usually paired with multiple historical streamflow sequences and the operation of the resource portfolio is simulated for each streamflow year individually. The use of historical streamflow sequences ensures that some level of risk analysis is incorporated in the modeling process. Some of the outputs from HYSIM including month end forebay elevation targets 56  for the reservoirs, energy production from thermal plants, Columbia River Treaty operation rule curves, etc. are used as inputs to GOM. 3.6.2 The Generalized Optimization Model (GOM) GOM is one of the primary tools used for operation planning studies at BC Hydro. It is used to assess the trading and operational benefits/impacts that result from changes in system operations, unit outages, efficiency upgrades, incorporation of new generation resources, transmission system upgrades, changes in constraints, etc. In this thesis, GOM is used to assess the impact of a change in the market bidding procedure. GOM is a deterministic, multi-period, variable time step linear programming (LP) optimization model. It includes a detailed hydraulic description of the hydropower plants in the BC Hydro system. GOM produces optimal generation and trading schedules that maximize the value of BC Hydro resources while balancing system load with resources and meeting operational constraints. In other words, it dispatches resources in the most economical manner to meet load without violating constraints. GOM is an extension of the short-term optimization model (STOM), which was developed by Shawwash (2000) for a similar purpose but a shorter planning horizon. STOM is used for optimization runs that are carried out at hourly timesteps with a planning horizon of less than seven days. GOM was extended to allow more flexibility in the time step and planning time horizon. Both GOM and STOM are written in the programming language AMPL (A Mathematical Programming Language), and they use the commercial optimization solver CPLEX. GOM usually optimizes the major hydropower plants while the rest of the generation in the system is considered as a fixed generation asset. The fixed generation resources include the smaller hydropower plants, thermal plants, and energy from IPPs. The smaller dispatchable hydropower plants can also be optimized in GOM if they are specifically being studied, but they are often treated as fixed generation since their flexibility and capacity is limited (Guzman 2010; BC Hydro 2013). Power from IPP’s including wind power and run-of-river hydro (ROR), are treated as a non-dispatchable but curtailable resource. Curtailing wind power and ROR generation might be required when dispatchable plants are operating at their minimum allowed generation, but there is still overgeneration (i.e., generation higher than demand), when transmission capacity is exhausted, or when releasing the wind 57  reserves (see Section 3.7.5  for detail) is more beneficial than using the wind power generated. For example, Evans (2009) discovered that curtailing wind for economic reasons can increase the value of wind energy in the BC Hydro system by up to 4%.  GOM requires many inputs such as market prices, system load, transmission capacity, non-dispatchable generation (i.e., ROR and wind), natural reservoir inflows, environmental water release requirements, power generation curves, generation limits, flow limits, etc. These inputs are obtained from various sources that include other models, specific studies, corporate forecast and historical data. GOM usually uses an hourly time step to reflect the hourly dispatch of energy and capture the variability of load, prices and other time-of-day dependent parameters. The forebay levels for storage reservoirs are fixed for the first and last hour of the optimization run. These forebay limits are obtained from long-term HYSIM studies that consider multiple water years. Other reservoir boundary conditions, such as the Columbia River Treaty operation, are also derived from the outputs of HYSIM.  The value of the GOM objective function is used to assess the economic impact of a change in system operation. The objective function of GOM is set to maximize the trade benefits from market transactions with the US and Alberta. The trade benefits are computed by multiplying the market prices for the Alberta and US trading hubs by export/import amounts. The short version of the objective function is given as:  𝑀𝑎𝑥 𝑅𝑒𝑣 ∶ ∑(𝐸𝑥𝑃𝑡  ∗ 𝐸𝑥𝑃𝑟𝑡 − 𝐼𝑚𝑃𝑡 ∗ 𝐼𝑚𝑃𝑟𝑡  − 𝐶𝑟𝑡𝑡 ∗ (𝑅𝐸𝐶 + 𝐸𝑥𝑃𝑟𝑡) − 𝑆𝑝𝑙𝑡 ∗ 𝑃𝑛𝑡𝑡)𝑇𝑡=1 23 where t is the time step, which is one hour in this thesis and T is the length of the planning horizon in hours. The first term in the equation represents the revenue generated from exporting energy (𝐸𝑥𝑃𝑡, in MWh) to the US and Alberta markets at the export price (𝐸𝑥𝑃𝑟𝑡, in $/MWh). The second term represents the cost of importing energy (𝐼𝑚𝑃𝑡, in MWh) from the US and Alberta at the import price (𝐼𝑚𝑃𝑟𝑡, in $/MWh). The import and export prices are different because of transmission (wheeling) charges and losses. The last two terms are penalties. The third term penalizes curtailment of the non-dispatchable renewable resource like wind power. The penalty is calculated by multiplying 58  the curtailed power amount (𝐶𝑟𝑡𝑡, in MWh) by the sum of the export price (𝐸𝑥𝑃𝑟𝑡) and the renewable energy credit (𝑅𝐸𝐶, in $/MWh) that is lost when curtailing wind power. The fourth term penalizes water spilled from the reservoirs (𝑆𝑝𝑙𝑡, in m3/s]). The spill penalty cost (𝑃𝑛𝑡𝑡, in $/m3/s) can be related to the marginal value of water at the marginal reservoir in the system such as the Williston reservoir behind W.A.C Bennett dam or the energy price at the time of the spill (BC Hydro 2013). The objective function value is the sum of these four hourly terms over the planning horizon (T, in hours). The planning horizon for GOM is usually one year but since this thesis deals with DA planning, a 24-hr planning horizon is used. GOM also has the flexibility to include a water value function in the objective function to represent the optimal tradeoff between energy traded in the market and the value of water held in storage. Some of the variables in GOM include the amount of power to be generated from the large hydropower plants, the amount of energy to be exported/imported and the amount of wind to curtail. The water released through turbines is calculated from the power generations using power generation production curves. The non-linear relationship between the water released and the power produced is expressed as piecewise linear functions. The relationship between storage and elevation is also defined this way. The main constraints of GOM include the load-resource balance, water balance, turbine and spillway flow limits, plant generation limits, forebay limits, available transmission capacity, and operating reserve requirements. Some of the important constraints are discussed below. The load resource balance puts forward the condition that at every time step in the planning horizon, the sum of the power generated plus imports should be greater than or equal to the domestic load plus exports. For J optimized plants, the constraint is written as:  (∑ 𝑃𝑗,𝑡𝐽𝑗=1) + 𝐹𝑥𝑑𝑃𝑡 + 𝑊𝑖𝑛𝑑𝑣𝑡 + 𝑅𝑜𝑅𝑡 + 𝐼𝑚𝑃𝑡 − 𝐸𝑥𝑃𝑡 ≥ 𝐿𝑜𝑎𝑑𝑡 , ∀𝑡 24 where 𝑃𝑗,𝑡 is the generation from the hydropower plant j, 𝐹𝑥𝑑𝑃𝑡  is the fixed generation that includes the small hydro and thermal, 𝑊𝑖𝑛𝑑𝑣𝑡  is the wind generation, 𝑅𝑜𝑅𝑡 is the ROR generation, 𝐼𝑚𝑃𝑡 is 59  the energy imported from the US and/or Alberta and  𝐸𝑥𝑃𝑡  is the export to US/Alberta. All terms in Equation 24 are in units of MWh.   The simplified version of the mass balance constraint is:  𝑉𝑡 =  𝑉𝑡−1 + (𝑄𝑛𝑎𝑡𝑡 + 𝑄𝑢𝑝𝑠𝑡 − 𝑄𝑡𝑢𝑟𝑡 − 𝑄𝑠𝑝𝑖𝑙𝑙𝑡) ∗ 𝑡𝑖𝑚𝑒, ∀𝑡 25 where 𝑉𝑡 is the storage in a reservoir at the current hour, 𝑉𝑡−1 is the storage from the previous hour, 𝑄𝑛𝑎𝑡𝑡 is the natural inflow, 𝑄𝑢𝑝𝑠𝑡 is the inflow that is obtained from upstream reservoirs, 𝑄𝑡𝑢𝑟𝑡 is the turbine release, 𝑄𝑠𝑝𝑖𝑙𝑙𝑡 is the water released over the spillway, and 𝑡𝑖𝑚𝑒 is the length of time step, which is one hour in this study. Different units can be used for components of Equation 25 but m3/s is usually used for the flow rates and the storage is usually converted to m3/s-day. The transmission constraints, which limit the amount of exported and imported energy, are given as:  𝐸𝑥𝑃𝑡 ≤ 𝑇𝑟𝑛𝑠. 𝐸𝑥𝑀𝑎𝑥𝑡  𝑎𝑛𝑑 𝐼𝑚𝑃𝑡 ≤ 𝑇𝑟𝑛𝑠. 𝐼𝑚𝑀𝑎𝑥𝑡 , ∀𝑡 26 where 𝑇𝑟𝑛𝑠. 𝐸𝑥𝑀𝑎𝑥𝑡  and 𝑇𝑟𝑛𝑠. 𝐸𝑥𝑀𝑎𝑥𝑡  are the available export/import transmission capacity in MWh or MW. Not that since the time step considered in this thesis is one hour, both units are sometimes used interchangeably for some of the parameters.  Most of the variables like power generation, turbine flows, reservoir storage, forebay, etc. have explicitly stated upper and lower bounds in GOM. A typical constraint that bounds a generic variable 𝑉𝑎𝑟𝑡  is given as:   𝑃𝑎𝑟𝑚𝑖𝑛,𝑡 ≤ 𝑉𝑎𝑟𝑡  ≤ 𝑃𝑎𝑟𝑚𝑎𝑥,𝑡 , ∀𝑡 27 where 𝑃𝑎𝑟𝑚𝑖𝑛,𝑡 is the minimum limit for the variable and 𝑃𝑎𝑟𝑚𝑎𝑥,𝑡 is the maximum limit for the variable.  The wind generation is also bounded as:  0 ≤ 𝑊𝑖𝑛𝑑𝑣𝑡 ≤ 𝑊𝑖𝑛𝑑𝑡 , ∀𝑡 28 60  where 𝑊𝑖𝑛𝑑𝑣𝑡 is the wind generation at time t after curtailment (in MWh)  and 𝑊𝑖𝑛𝑑𝑡 is the total wind power generated from the wind farms (in MWh). This constraint gives the model flexibility to curtail wind when it is not needed.  Equations 23 to 28 describe some of the main components of the GOM model. These equations, along with many other variables, parameters and constraints are assumed to characterize the workings of the BC Hydro system.  3.7 Modifications Made to the Generalized Optimization Model (GOM) The GOM model currently in use at BC Hydro was modified to perform the required studies in this thesis. As mentioned in Chapter 1, one of the goals of the thesis is to use SO to reduce the DA wind integration cost. To achieve this goal, different types of models have to be compared. The types of models that are compared to estimate the value of SO can generally be divided into two models: the Two-stage deterministic models and the Stochastic model, as discussed below: i. Two-stage deterministic models: in this type of model, wind power generation in the DA is represented by single forecast and is assumed to be deterministic. The DA decisions are made based on the wind generation forecast, but the RT operation is carried out based on the actual wind power generation.  ii. Stochastic models: in these types of models, the wind is considered uncertain and is represented by different wind generation scenarios instead of a single deterministic forecast. The DA decisions are made based on the wind generation scenarios, and RT operations are based on the actual wind power generation.  Note that in both models it is assumed that wind is known with certainty in RT system operations. This is a reasonable assumption as the system operator typically has full knowledge of the magnitude of current wind generation in the system in real-time.  These two types of models require different modifications. Since most of the needed modifications are similar, they are discussed below together in the same subsections. When the difference in the modifications between models is significant, they are discussed separately. Some of the modifications made to GOM are discussed below under three sections.  61  3.7.1 Modification to Optimization Algorithm For type (i) models, GOM was converted from a one-stage deterministic optimization model into a two-stage deterministic optimization model. This was done by fixing the DA offers based on an optimization run that considers the wind generation forecast and running RT operations with the actual wind power realization and the fixed DA bids. This approach correctly represents the day to day operational procedures used at BC Hydro.  One-stage deterministic linear programming models can be represented in the following canonical form:  𝑀𝑎𝑥𝑖𝑚𝑖𝑧𝑒 𝑍 = 𝒄𝑇𝒙     29  𝑨𝒙 ≤ 𝑩, 𝒙 ≥ 0 30 where 𝑍 is the objective function, 𝒄𝑇 is the coefficient matrix for the decision variable matrix 𝒙, and 𝑨 and 𝑩 are known matrices for the corresponding constraints. Meanwhile, two-stage deterministic optimization problems can be represented in the following canonical form:  𝑀𝑎𝑥𝑖𝑚𝑖𝑧𝑒 𝑍2 = 𝒄𝑇𝒙  + 𝒅𝑇𝒚  31  𝑨𝒙 ≤ 𝑩 32  𝑻𝒙 + 𝑾𝒚 ≤ 𝑯 33   𝒙 ≥ 0, 𝒚 ≥ 0 34 where Z2 is the objective function, c is the coefficient of the first stage decision matrix x, d is the coefficient of the second stage decision matrix y, and A, B, T, W, H are all known matrices. This formulation works as a two-stage optimization problem using the following steps: Step 1. First, the optimization is run with forecasted W and H values. An optimal value based on the forecast is obtained along with values for x and y. Step 2. The values for x obtained from step 1 are fixed. These are the first stage decisions. The values for y are discarded. Step 3. The optimization is rerun with the fixed x values and actual W and H values. This gives the values for y (i.e., the second stage decisions) and the actual optimal value. 62  For type (ii) models, GOM is adapted into a two-stage stochastic LP model. As discussed in the literature review, two-stage stochastic optimization algorithms are often used to solve DA scheduling and market clearing problems. In two-stage stochastic optimization problems, the first stage decisions are also called here and now decisions while the second stage decisions are called recourse decisions. The difference between the deterministic type (i) models and the stochastic type (ii) models is that the first stage decisions in type (ii) models are fixed based on the wind generation scenarios instead of the single deterministic forecast. Two-stage stochastic linear programming problems can be represented in the following canonical form:  𝑚𝑎𝑥𝑖𝑚𝑖𝑧𝑒 𝑍𝑠𝑡𝑜𝑐 = 𝒄𝑇𝑥 + ∑ 𝜋(𝑤)𝒒(𝑤)𝑇𝒚(𝑤)𝑤∈𝛺 35  𝑨𝑥 ≤ 𝑩 36  𝑻(𝑤)𝒙 + 𝑾(𝑤)𝒚(𝑤) ≤ 𝑯(𝑤), ∀𝑤 ∈ 𝛺 37  𝒙 ∈ 𝑋, 𝒚(𝑤) ∈ 𝑌, ∀𝑤 ∈ 𝛺 38 where c is as discussed above, x is the first stage decision variable matrix, w is the scenario in the scenario set Ω, π(w) is the probability of scenario w, q(w) is the coefficient matrix of the second stage variables and y(w) is the second stage decision variable matrix that is dependent on scenario w. The matrices T(w) and W(w) and H(w) are known matrices that may or may not depend on the scenario w (Conejo et al. 2010). Equation 35 specifies that the objective function of a two-stage optimization problem is the sum of two parts: the revenue from the first stage decisions and the expectation of the revenues from the second stage decisions. Equation 37  indicates that the constraints are specified for each scenario.  Stochastic LP models have many advantages over other stochastic models such as mixed integer programming, dynamic programming or nonlinear programming models. Some of these advantages include faster solution time, less computational burden, guaranteed global optimal solution, it does not need initial solutions, the availability of readily available solvers with proven performance, and sensitivity analysis information. However, stochastic LP models also have some disadvantages such as its inability to accurately capture non-linear relationships and its difficulty with modeling multi-stage optimization problems  (Goodarzi et al. 2014; Chong & Zak 2013; Labadie 2004). 63  In this thesis, the first stage decisions (x) are the DA market transactions, which have to be made before the wind generation is known with certainty. Meanwhile, the second stage decisions (y) are the RT market transactions and system operation decisions, which are made after the wind power generation is known or can be accurately predicted. In addition, when the distinction is important, the optimal values from the first stage runs are called ‘expected optimal values’ and the optimal values from the second stage runs are called ‘actual optimal values’. 3.7.2 Added Variables and Parameters The version of GOM currently in use does not distinguish between DA market and the real-time (RT) (also called spot) market. This is because the purpose of GOM is to conduct planning studies that assume perfect forecast and do not require this distinction. Currently in GOM, all the market transactions are assumed to take place in the real-time markets. However, only a small portion of the energy is traded in the RT markets. BC Hydro, through its energy trading subsidiary Powerex, engages in different types of market transactions (including long-term bilateral (forward) contracts, DA transactions, and RT transactions) with various trading partners in US and Canada. Most of the transactions modeled in GOM are assumed to be finalized at the Mid-C trading hub in the forward and DA markets. The price volatility and the difficulty of securing transmission access in real-time restrict Powerex’s ability to heavily trade in the RT markets.  Although GOM currently does not distinguish between DA and RT markets, it can be modified to accommodate both when the distinction between the two is important. For example, Guzman (2010) separated the markets into DA and RT markets. When looking at the actual scheduling process and the uncertainties involved, the distinction between the DA market and RT market is important. Although it is desirable to do most of the transactions in the DA market, the level of uncertainty encountered on the DA time frame necessitates a secondary RT market where differences between the DA commitments and RT load/generations can be balanced. The RT market transactions are usually concluded 20 minutes before the start of the scheduled hour and therefore face much lower levels of uncertainty. Any deviations from the hour-ahead schedules are balanced using RT trades and ancillary services which include regulation and spinning reserves and other capacity contracts such as dynamic schedules.  64  Since the focus of this thesis is on improving the DA scheduling process, the distinction between the DA and RT markets is vital. To reflect this in GOM, new variables that represent the DA and RT markets are defined. A new parameter representing the RT market prices was also added. The RT market prices are subject to high price volatility and are influenced by many factors such as the deviation of the RT load from DA commitments, unit outages, high forecast errors of non-dispatchable generation, etc. For this reason, modeling the RT market price is a complicated process. Because of this complication, and given that the focus of this thesis to isolate the value of SO in reducing the cost of DA wind uncertainty, the RT market prices were considered to be the same as the DA market prices. This is a reasonable assumption given that RT average prices tend to converge to DA in the longterm (i.e., tend to revert to the mean DA prices) as can be seen in the 2017 California Independent System Operator (CAISO) quarterly market reports. The ancillary service (AS) market is also modeled in this thesis. This is especially appropriate given that BC Hydro and Powerex have started participating in the AS market. Two types of AS are considered in this thesis: regulation and spinning reserves. Reserves, in general, are generation capacity withheld to ensure reliable operations. Regulating reserves address minute-to-minute imbalances while spinning reserves address imbalances in longer, within-the-hour time spans. AS can be provided in the up directions (ramping up generation) and in the down directions (ramping down generation) (BC Hydro 2013). The net energy used while providing up and down reserves is expected to net out to zero, so it was not considered in this study. Considering the AS market is important because the amount of wind reserves that are withheld limits BC Hydro’s ability to offer AS in the market (see Section 3.7.5 for more detail).  Parameters regarding the marginal value of water (MVW) are also introduced. Consideration for MVW is necessary because fixing the forebay levels at the first and last hour of the optimization, as is currently done in GOM for runs with longer time horizons, is not feasible for runs with 24-hr planning horizon. This is because the daily end forebay limits excessively restrict the operation of the system and can cause infeasibilities. Therefore, MVW is used in conjunction with daily storage targets to guide the operation of the hydropower plants. The MVW is estimated from the HYSIM system marginal cost and is made unique for each plant by factoring in the generation capability of each plant (e.g., water in a high head plant is more valuable than water in a low head plant). The daily storage targets are derived from HYSIM monthly forebay targets using the following steps: 65  1. The targets from HYSIM were used in the deterministic GOM runs with one-year time horizon. 2. The ending forebay levels of the GOM runs (i.e., the forebay at the 24th hour) for all days of the year were recorded and made the new daily forebay targets. 3. The daily storage targets were derived from the daily forebay targets by using the piecewise linear elevation-storage curves. Any deviation from these daily storage targets was penalized if the ending storage was less than the target storage or rewarded if the ending storage was higher than the target storage. The cost of deviation was the marginal value of water in the reservoir where deviation occurred (see Equation 39 below).   Some of the existing parameters and variables that are dependent on the outcome of the stochastic process are also indexed by the scenarios. In other words, some of the parameters and variables are made unique for each realization of the stochastic process, i.e., the wind generation. Some of the indexed variables include the wind generation, the hydropower plant generation, forebay levels, and storage. The indexing of parameters and variables also leads to the changes in the objective function and constraints.  3.7.3 Modification to the Objective Function The objective function was modified to include the trade benefits from both the DA and RT markets. The DA market transactions are taken as the first stage decisions and hence are not indexed over the scenarios. The RT market transactions are considered as the second stage decisions and are therefore indexed over the scenarios.  The revenue in the SO’s objective function is a random variable because it depends on the wind generation random variable. Some wind generation scenarios (e.g., lower than forecasted wind generation) might lead to more imports and negative revenue while other (e.g., scenarios with higher wind generation than forecasted) might result in more exports and positive revenue. The expectation of this random variable is the objective function. For the planning horizon T, with the subscripts t, s, and j representing the hour, scenario, and hydropower plant respectively, the stochastic objective function is stated as: 66  𝑀𝑎𝑥𝑖𝑚𝑖𝑧𝑒 𝑅𝑒𝑣_𝑆𝑡𝑜𝑐: ∑(𝐷𝐴. 𝐸𝑥𝑃𝑡 ∗ 𝐷𝐴. 𝐸𝑥𝑃𝑟𝑡 − 𝐷𝐴. 𝐼𝑚𝑃𝑡 ∗ 𝐷𝐴. 𝐼𝑚𝑃𝑟𝑡 + 𝐷𝐴. 𝐴𝑆𝑢𝑝𝑡𝑇𝑡=1∗ 𝐷𝐴. 𝐴𝑆𝑢𝑝𝑃𝑟𝑡 + 𝐷𝐴. 𝐴𝑆𝑑𝑛𝑡 ∗ 𝐷𝐴. 𝐴𝑆𝑑𝑛𝑃𝑟𝑡) + ∑ 𝑝𝑠 ∗ ∑(𝑇𝑡=1𝑆𝑠=1𝑅𝑇. 𝐸𝑥𝑃𝑡,𝑠∗ 𝑅𝑇. 𝐸𝑥𝑃𝑟𝑡 − 𝑅𝑇. 𝐼𝑚𝑃𝑡,𝑠 ∗ 𝑅𝑇. 𝐼𝑚𝑃𝑟𝑡 + 𝑅𝑇. 𝐴𝑆𝑢𝑝𝑡,𝑠 ∗ 𝑅𝑇. 𝐴𝑆𝑢𝑝𝑃𝑟𝑡+ 𝑅𝑇. 𝐴𝑆𝑑𝑛𝑡,𝑠 ∗ 𝑅𝑇. 𝐴𝑆𝑑𝑛𝑃𝑟𝑡 − 𝐶𝑟𝑡𝑡,𝑠 ∗ (𝑅𝐸𝐶 + 𝑅𝑇. 𝐸𝑥𝑃𝑟𝑡) − 𝑆𝑝𝑙𝑡,𝑠 ∗ 𝑃𝑛𝑙𝑡− 𝐿𝑜𝑎𝑑𝑆ℎ𝑒𝑑𝑡,𝑠 ∗ 𝑃𝑛𝑙𝑆ℎ𝑒𝑑𝑡 + ∑(𝑀𝑉𝑊𝑗,𝑠 ∗ (𝑉𝑒𝑛𝑑𝑗,𝑠 − 𝑉𝑡𝑎𝑟𝑔𝑒𝑡𝑗,𝑠𝐽𝑗=1)) 39 where, • 𝐷𝐴. 𝐸𝑥𝑃𝑡  is the DA export, in MWh,  • 𝐷𝐴. 𝐸𝑥𝑃𝑟𝑡 is the DA export price (in $/MWh,  • 𝐷𝐴. 𝐼𝑚𝑃𝑡 is the DA import,  in MWh, • 𝐷𝐴. 𝐼𝑚𝑃𝑟𝑡  is the DA import price, in $/MWh, • 𝐷𝐴. 𝐴𝑆𝑢𝑝𝑡  is the up AS sold in the DA market, in MW, • 𝐷𝐴. 𝐴𝑆𝑢𝑝𝑃𝑟𝑡  is the DA up AS price, in $/MW,  • 𝐷𝐴. 𝐴𝑆𝑑𝑛𝑡 is the down reserves sold in the DA market, in MW,  • 𝐷𝐴. 𝐴𝑆𝑑𝑛𝑃𝑟𝑡is the down AS price, in $/MW, • 𝑝𝑠 is the probability of the scenario s,  • 𝑅𝑇. 𝐸𝑥𝑃𝑡,𝑠 is the RT export, in MWh, • 𝑅𝑇. 𝐸𝑥𝑃𝑟𝑡  is the RT export price, in $/MWh, • 𝑅𝑇. 𝐼𝑚𝑃𝑡,𝑠  is the RT import, in MWh, • 𝑅𝑇. 𝐼𝑚𝑃𝑟𝑡  is the RT import prices, in $/MWh, • 𝑅𝑇. 𝐴𝑆𝑢𝑝𝑡,𝑠 is the up AS sold in the RT market, in MW,  • 𝑅𝑇. 𝐴𝑆𝑢𝑝𝑃𝑟𝑡  is the RT up AS price, in $/MW,  • 𝑅𝑇. 𝐴𝑆𝑑𝑛𝑡,𝑠 is the down reserves sold in the RT market, in MW,  • 𝑅𝑇. 𝐴𝑆𝑑𝑛𝑃𝑟𝑡 is the down AS price in the RT market, in $/MW,  • 𝐿𝑜𝑎𝑑𝑆ℎ𝑒𝑑𝑡,𝑠 is the load shedding, in MWh (curtailed load amount),  67  • 𝑃𝑛𝑙𝑆ℎ𝑒𝑑𝑡  is the penalty for curtailing load, in $/MWh (see Section 3.7.4 below for details about load shedding),  • 𝑀𝑉𝑊𝑗,𝑠 is the Marginal Value of Water in ($/(m3/s))-day,  • 𝑉𝑒𝑛𝑑𝑗,𝑠 is the storage, in m3/s-day at the last hour of the panning horizon,  • 𝑉𝑡𝑎𝑟𝑔𝑒𝑡𝑗,𝑠 is the storage target in m3/s-day for the last hour of the planning horizon.  Other variables and parameters are as described in sections 3.6.2 and 3.7.2.  Since the objective function presented in Equation 39 maximizes the expectation of the revenue random variable, the first stage decisions that result from it are suboptimal for individual wind scenarios but are the best decisions for the wind generation uncertainty captured in the scenarios. The formulation of the objective function is risk neutral since it only considers the expected value of the revenue random variable. This formulation is expected to yield higher benefits in the long run compared to other formulations that take risk into consideration (Conejo et al. 2010). For the deterministic type (i) models, Equation 39 is applied to a single scenario (i.e., the wind forecast) that is assigned a probability of 1.  3.7.4 Modifications to the Constraints For type (ii) models, the constraints that contain stochastic parameters or variables are modified to ensure that all constraints are respected for all scenarios. For example, the load-resource balance is restated as:  (∑ 𝑃𝑗,𝑡,𝑠𝐽𝑗=1 ) + 𝐹𝑥𝑑𝑃𝑡 + 𝑊𝑖𝑛𝑑𝑣𝑡,𝑠 + 𝑅𝑜𝑅𝑡 + 𝐷𝐴. 𝐼𝑚𝑃𝑡 + 𝑅𝑇. 𝐼𝑚𝑃𝑡,𝑠 −𝐷𝐴. 𝐸𝑥𝑃𝑡 −  𝑅𝑇. 𝐸𝑥𝑃𝑡,𝑠 + 𝐿𝑜𝑎𝑑𝑆ℎ𝑒𝑑𝑡,𝑠 ≥ 𝐿𝑜𝑎𝑑𝑡 , ∀𝑠, ∀𝑡  40 This constraint ensures that no matter what the wind power output in scenario s, the load is served using the combination of available resources. The load shed (𝐿𝑜𝑎𝑑𝑆ℎ𝑒𝑑𝑡,𝑠) variable is included to avoid infeasibilities and is penalized at a high cost in the objective function so that it only occurs when load curtailment is the only way of maintaining feasibility.  The water balance equation is redefined as: 68   𝑉𝑡,𝑠 =  𝑉𝑡−1,𝑠 + 𝑄𝑛𝑎𝑡𝑡 + 𝑄𝑢𝑝𝑠𝑡,𝑠 − 𝑄𝑡𝑢𝑟𝑡,𝑠 − 𝑄𝑠𝑝𝑖𝑙𝑙𝑡,𝑠, ∀𝑠, ∀𝑡 41 This constraint shows that the storage is dependent on the wind power scenarios. The transmission constraints are restated to include the DA and RT transactions as:  𝐷𝐴. 𝐸𝑥𝑃𝑡 +   𝑅𝑇. 𝐸𝑥𝑃𝑡,𝑠 + 𝐷𝐴. 𝐴𝑆𝑢𝑝𝑡 + 𝑅𝑇. 𝐴𝑆𝑢𝑝𝑡,𝑠 ≤ 𝑇𝑟𝑛𝑠. 𝐸𝑥𝑀𝑎𝑥𝑡  , ∀𝑠, ∀𝑡 𝑎𝑛𝑑  𝐷𝐴. 𝐼𝑚𝑃𝑡 + 𝑅𝑇. 𝐼𝑚𝑃𝑡,𝑠 + 𝐷𝐴. 𝐴𝑆𝑑𝑛𝑡 + 𝑅𝑇. 𝐴𝑆𝑑𝑛𝑡,𝑠 ≤ 𝑇𝑟𝑛𝑠. 𝐼𝑚𝑀𝑎𝑥𝑡 , ∀𝑠, ∀𝑡 42 These constraints ensure that there is enough transmission capacity to trade electricity in the DA and RT markets in addition to having the capability to serve AS.  Constraints that limit the RT transactions by the market liquidity are also added to the model. These constraints are presented as:                            𝑅𝑇. 𝐸𝑥𝑃𝑡,𝑠 + 𝑅𝑇. 𝐴𝑆𝑢𝑝𝑡,𝑠 ≤ 𝑅𝑇. 𝐸𝑥𝑀𝑎𝑥𝑡   , ∀𝑠, ∀𝑡 𝑎𝑛𝑑  𝑅𝑇. 𝐼𝑚𝑃𝑡,𝑠 + 𝑅𝑇. 𝐴𝑆𝑑𝑛𝑡,𝑠 ≤ 𝑅𝑇. 𝐼𝑚𝑀𝑎𝑥𝑡 , ∀𝑠, ∀𝑡 43 Since these constraints put upper bound on the possible RT market transactions, they can be used to control the magnitude of trading in the RT market. Upper and lower bound pertaining to several of other variables are presented in their stochastic form as:   𝑃𝑎𝑟𝑚𝑖𝑛,𝑡 ≤ 𝑉𝑎𝑟𝑡,𝑠  ≤ 𝑃𝑎𝑟𝑚𝑎𝑥,𝑡  ∀𝑠, ∀𝑡 44 These types of constraints ensure that the upper and lower bound constraints are not violated for all scenarios.   The wind generation is bounded by:  0 ≤ 𝑊𝐼𝑁𝐷𝑉𝑡,𝑠 ≤ 𝑊𝑖𝑛𝑑𝑡,𝑠∀𝑠, ∀𝑡 45 For type (i) models, Equations 40 to 45 are stated for only one scenario as discussed above.  69  3.7.5 Modification to Operating Reserves Operating reserves (ORs) are generating capacity withheld to ensure reliable operation of an electric system. ORs are called upon whenever there is an imbalance between planned generation and demand. Up reserves are needed when the demand is higher than the generation and down reserves are required when the generation is higher than the demand. It is noted here that the AS service discussed in Section 3.7.2 are reserves offered for the AS market while the reserves discussed in this section are the AS required to balance BC Hydro’s system. GOM currently accounts for contingency and load reserves as part of ‘operating reserve obligations’ (ORO). The ORO is the percentage of the generating capacity that is reserved from various generation resources to handle forced outages and large load forecast errors. The ORO does not account for wind reserves, so they are considered separately. The current GOM only considers within-the-hour reserves, which are calculated following the methodology detailed in BC Hydro (2013). Since they were not the focus of this thesis, they were not modified.  In the current model, reserve constraints are formulated to ensure that there is enough generation capacity to manage the variability of wind power. The wind up reserve constraint is given as:  ∑(𝑃𝑚𝑎𝑥𝑗,𝑡 − 𝑃𝑗,𝑡 − 𝑃𝑗,𝑡 ∗ 𝑂𝑅𝑂)𝐽𝑗=1≥ (1 − 𝐶𝑟𝑡. 𝐹𝑟𝑐𝑛𝑡) ∗ 𝑊𝑡𝑜𝑡𝑎𝑙. 𝑢𝑝. 𝑟𝑒𝑠𝑡 , ∀𝑡 46 where 𝑃𝑚𝑎𝑥𝑗,𝑡 is the maximum possible generation of plant j at time t, 𝑃𝑗,𝑡 is the power generated by plant j at time t, 𝑂𝑅𝑂 is the operating reserve obligation, 𝐶𝑟𝑡. 𝐹𝑟𝑐𝑛𝑡  is the fraction of wind curtailed calculated as (𝑊𝑖𝑛𝑑𝑡,𝑠 − 𝑊𝐼𝑁𝐷𝑉𝑡,𝑠)/𝑊𝑖𝑛𝑑𝑡,𝑠, and 𝑊𝑡𝑜𝑡𝑎𝑙. 𝑢𝑝. 𝑟𝑒𝑠t is the total within-the-hour up wind reserve needed at time t. This constraint ensures that the capacity left at the plant after generation and load reserves are deducted is enough to provide the wind up reserves. The inclusion of 𝐶𝑟𝑡. 𝐹𝑟𝑐𝑛𝑡 in the constraint is to ensure that if the wind power is curtailed, the wind reserve that has to be kept is also reduced.  The wind down reserves constraint is provided by: 70   ∑(𝑃𝑗,𝑡 − 𝑃𝑚𝑖𝑛𝑗,𝑡)𝐽𝑗=1 ≥  (1 − 𝐶𝑟𝑡. 𝐹𝑟𝑐𝑛𝑡)   ∗ 𝑊𝑡𝑜𝑡𝑎𝑙. 𝑑𝑜𝑤𝑛. 𝑟𝑒𝑠𝑡 , ∀𝑡 47 where 𝑃𝑚𝑖𝑛𝑗,𝑡 is the minimum allowed generation for plant j at time t and 𝑊𝑡𝑜𝑡𝑎𝑙. 𝑑𝑜𝑤𝑛. 𝑟𝑒𝑠𝑡 is the total within-the-hour down wind reserve. This constraint ensures that there is enough room for plant generation to move down in case the wind power output is higher than the forecast.  In this thesis, DA wind reserves are also required in addition to the within-the-hour wind reserves. The DA wind reserves were calculated independently for each month of the year. For each month, the amount of reserve that is expected to ensure a reliable operation of the BC Hydro system is set to be the amount that covers 99.7% (3 stds) of the observed DA forecast errors. This is in accordance with the methodology used in BC Hydro (2013). A liquidity allowance for 300 MW is made since that amount of power can be purchased in the real-time market without incurring significant financial losses. The resulting reserve amount from this calculation is considered to be the incremental wind reserve required in addition to the load and contingency reserves captured in the ORO. Once the reserves are calculated, they are incorporated into Equations 46 and 47 by making the terms 𝑊𝑡𝑜𝑡𝑎𝑙. 𝑢𝑝. 𝑟𝑒𝑠t  and 𝑊𝑡𝑜𝑡𝑎𝑙. 𝑑𝑜𝑤𝑛. 𝑟𝑒𝑠𝑡  in the equations equal to the sum of the within-the-hour and DA reserves. When the DA reserves are not considered in the model, the terms will only refer to the within-the-hour reserves. Equations 46 and 47 are modified for the stochastic type (ii) models by indexing 𝑃𝑗,𝑡 and 𝐶𝑟𝑡. 𝐹𝑟𝑐𝑛𝑡 by the scenarios. 3.8 Estimation of Stochastic Savings As mentioned in Chapter 1, stochastic savings are defined as the percentage of the cost that can be saved by using SO instead of deterministic optimization. The stochastic savings that resulted from using a bidding strategy based on the two-stage stochastic optimization algorithm is investigated following the methodology detailed in Tuohy et al. (2008), Greenhall (2013), and Ruiz et al. (2009). Four types of models are used to compare the savings due to SO: a. Deterministic model without reserves:  in these models, a perfect forecast is assumed, and no DA reserves are kept. These models represent an ideal condition where the wind has no uncertainty and therefore provide the lower bound in the cost of wind integration.  71  b. Deterministic model with reserves:  here the wind forecast is assumed to be uncertain and DA reserves are kept to manage the WFE. The DA trade decisions are made using the wind forecast with DA reserves. These DA reserves are then released in RT operations so that the system will have the capacity needed to manage the WFE. This model approximates the current operation practices at BC Hydro.  c. Stochastic model with MC generated scenarios and without reserves:  here the wind is considered to be uncertain, but no reserves are explicitly kept to manage the WFE. The model decides on DA market transactions after assessing the WFE scenarios created by the MC model. In RT operations, the wind generation scenarios are replaced by the actual wind generation. The cost of this model is expected to be between models discussed in (a) and (b). d. Stochastic model with ARMA generated scenarios and without reserves: same as in model (c) but the scenarios for the SO are obtained from the ARMA model. The model in (b) is expected to have the lowest optimal value (i.e., maximum cost) while the one in (a) is expected to have the highest optimal value (i.e., lowest cost). Given that model (b) resembles current BC Hydro operations the most, the DA wind integration cost in this thesis is defined as the difference between the cost of model (b) and the cost of model (a). A similar approach model was used to determine the DA wind integration cost by BC Hydro (2013). Comparing the optimal values of the models in (c) and (d) to the optimal value of the model in (b) gives the value of SO in the BC Hydro system. When analyzing the value of SO, the term ‘relative cost’ is used in addition to stochastic saving. Relative cost is defined as the cost of a stochastic model relative to the model in (b). Relative cost is related to stochastic saving since stochastic saving is equal to 100% minus relative cost. The effect of the scenario generation algorithm and the quality of scenarios on stochastic savings is investigated by comparing the optimal values of the model in (c) to the model in (d). Some of the equations that are used to estimate the stochastic savings are given below.    𝐶𝑜𝑠𝑡𝑖 = ∑(𝑅𝑒𝑣𝑎,𝑗𝐷𝑑=1− 𝑅𝑒𝑣𝑖.𝑗), 𝑖 ∈ {𝑎, 𝑏, 𝑐, 𝑑} 48 72   𝑊𝑖𝑛𝑑𝐼𝑛𝑡𝑒𝑔𝑟𝑎𝑡𝑖𝑜𝑛𝐶𝑜𝑠𝑡𝑖 =𝐶𝑜𝑠𝑡𝑖𝑇𝑜𝑡𝑎𝑙𝐴𝑛𝑛𝑢𝑎𝑙𝑊𝑖𝑛𝑑𝐺𝑒𝑛 , 𝑖 ∈ {𝑏, 𝑐, 𝑑} 49  𝑅𝑒𝑙𝑎𝑡𝑖𝑣𝑒𝐶𝑜𝑠𝑡𝑖 =∑ (𝑅𝑒𝑣𝑎,𝑗𝐷𝑗=1 − 𝑅𝑒𝑣𝑖.𝑗)𝐶𝑜𝑠𝑡𝑏 =𝐶𝑜𝑠𝑡𝑖𝐶𝑜𝑠𝑡𝑏, 𝑖 ∈ {𝑐, 𝑑} 50  𝑆𝑡𝑜𝑐𝑎𝑠𝑡𝑖𝑐𝑆𝑎𝑣𝑖𝑛𝑔𝑖 = 100% − 𝑅𝑒𝑙𝑎𝑡𝑖𝑣𝑒𝐶𝑜𝑠𝑡𝑖, 𝑖 ∈ {𝑐, 𝑑}  51 In the equations given above, the subscript j represents the day, D is the number of days in the study period, 𝐶𝑜𝑠𝑡𝑖 is the cost of model (i), and 𝑅𝑒𝑣𝑖.𝑗 is the revenue obtained from model (i) on day j. Both the cost and revenue in the above equations are in units of Canadian dollars ($) (CAD). The wind integration cost (𝑊𝑖𝑛𝑑𝐼𝑛𝑡𝑒𝑔𝑟𝑎𝑡𝑖𝑜𝑛𝐶𝑜𝑠𝑡𝑖 in $/MWh) is calculated by dividing 𝐶𝑜𝑠𝑡𝑖 by the total annual wind generation (𝑇𝑜𝑡𝑎𝑙𝐴𝑛𝑛𝑢𝑎𝑙𝑊𝑖𝑛𝑑𝐺𝑒𝑛, in MWh). Note that the cost of model (b) (𝐶𝑜𝑠𝑡𝑏), is used to calculate the relative costs (𝑅𝑒𝑙𝑎𝑡𝑖𝑣𝑒𝐶𝑜𝑠𝑡𝑖) of models (c) and (d). The stochastic saving of model (i) (𝑆𝑡𝑜𝑐𝑎𝑠𝑡𝑖𝑐𝑆𝑎𝑣𝑖𝑛𝑔𝑖) is calculated as 100% minus the relative cost of model (i). As discussed in Section 3.7.1, the models in (a) to (d) have two optimal values, the expected optimal values, and the actual optimal values. Most of the published papers, including the ones discussed in Chapter 2, use the expected optimal values to assess the performance of stochastic models. Consequently, an assessment that uses the expected optimal values is conducted in this thesis. However, more focus is given to the assessment that is based on the actual optimal values. According to Conejo et al. (2010), an assessment that uses the actual optimal values is called out-of-sample assessment. In out-of-sample assessments, the performance of the decisions obtained from different types of models is compared in terms of actual real-world results. Most of the assessment is done after the daily optimal values obtained from models (a) to (d) are summed for the study year. Note that the term ‘out-of-sample assessment’ is used in this thesis to refer to the performance assessment of the SO model based on the actual optimal values and not the cross-validation of the scenario generation models with split samples. The term ‘out-of-sample’ is used in the sense that the actual wind generation is thought to be independent of the wind generation scenarios used in the SO (it is outside of the sample of wind generation scenarios used in the SO). Sensitivity analysis is also conducted to assess the effects reservoir inflow conditions, per unit load shedding cost, and AS market depth assumption will have on the stochastic savings. 73   Case Study and Inputs In this chapter, a case study which is used to investigate the performance of the proposed models outlined in the previous chapter is presented. In the first section, a general overview of the BC Hydro system is presented.  The second and third sections then discuss the various modeling assumptions and data that were used as inputs in the optimization models. The fourth section describes the different cases investigated in this thesis. 4.1 The BC Hydro System BC Hydro is crown cooperation and power utility owned by the people and government of British Columbia. BC Hydro generates and delivers electric power to more than 1.9 million homes, businesses, and industries, serving more than 95% of the province’s population. BC Hydro is a vertically integrated utility, which means it owns the generation, transmission, and distribution assets (BC Hydro 2017).  Currently, the BC Hydro system is comprised of 32 generating stations, 30 hydropower plants and two thermal plants. These generating stations, which have a total capacity of more than 12000 MW, produce on average 46000 GWh of energy. More than 98% of this energy is generated from renewable sources. The Columbia River region of BC is responsible for 58% of BC Hydro’s capacity and 48% of its energy production while the Peace Region is responsible for 29% of capacity and 38% of generation. Lower Mainland provides about 9% of capacity and 10% of generation. BC Hydro purchases about 21000 GWh of energy annually from IPPs that operate more than 110 hydropower, biomass, and wind plants (BC Hydro 2017). As discussed in Section 2.1 wind IPPs have a total capacity of 672 MW, and generate more than 2000 GWh of energy. Figure 8 shows the location of some of the large renewable power generating plants in BC.  BC Hydro’s generation from its main generation areas in the Peace and Columbia regions is delivered to its customers that are mostly concentrated in the Lower Mainland through 18000 km of high voltage transmission lines and 56000 km of lower voltage distribution lines. Voltage conversion is carried out in over 300 substations spread across the province. BC Hydro also has transmission interties that connect it to the US and Alberta for the trading of electric energy and AS (BC Hydro 2017).  74   Figure 8 - Renewable generation resources in British Columbia Source: National Energy Board (2016) As discussed in Section 3.6.2, GOM is used to model the operation of the BC Hydro system. Some of the model inputs are discussed below.  4.2 Model Inputs and Assumptions Some of the input data required for the optimization models include load data, inflows, market prices, various limits on generation, transmission, flow, storage, reservoir elevation, flow ramps, storage ramps, generation ramps, etc. The sources and properties of some of these data are discussed below. 4.2.1 Load The study year considered for this thesis is the year 2035. It is noted here that the term ‘year’, when used by BC Hydro in the context of load year, water year, wind year or study year, refers to a period spanning from Oct. 1 to Sep. 30 (e.g., load year 2035 is defined as the period from Oct 1, 2035 to 75  Sep 30, 2036). This definition of a year is selected to maximize the correlation between annual precipitation and annual reservoir inflow values.  The load year 2035 is used because it is the year when the system is expected to be in load-resource balance. The hourly load data for the study year 2035 was derived by first using a monthly load forecast for the year 2035 and shaping it into hourly values with a typical BC Hydro hourly load shape. The BC Hydro system load data resulting from this process is shown in Figure 9. The load is highest in the winter months and lowest in the summer months.   Figure 9 - BC Hydro load forecast for 2035 4.2.2 Reservoir Inflows The natural inflow records for the various reservoirs were obtained from historical records. To determine the value of SO in the BC Hydro system under different inflow conditions, three specific water years with normal, wet, and dry inflow conditions were chosen for this thesis. Figure 10 shows the cumulative inflows into the large hydropower plant’s reservoirs for the water years 1964 to 1973. These years are considered because they show the variation in inflow that can be expected in the BC Hydro system (BC Hydro 2013; Evans 2009; Guzman 2010). The three water years for this thesis were selected from these 10 years; 1969/70 was considered as a dry year, 1968/69 a normal water year, and 1966/1967 a wet year.  020004000600080001000012000140001600018000Oct Nov Dec Jan Feb Mar Apr May Jun Jul Aug SepLoad (MW)76  The average inflows into the largest reservoir in the BC Hydro system, Williston Reservoir, is shown in Figure 11. The pattern, which is typical of for most of BC Hydro’s reservoirs, indicates that most of the inflows are generated when the accumulated snow melts in the freshet period (i.e., late spring/early summer months).    Figure 10 - Cumulative inflows for BC Hydro's large reservoirs (1964-1973)  Source: Evans (2009)  Figure 11 - Average inflows into the Williston Reservoir (GMS) (1930-2009) 01000200030004000500060007000Oct Nov Dec Jan Feb Mar Apr May Jun Jul Aug SepInflow (m3/s)77  4.2.3 Price  Prices for electricity, natural gas and AS are needed for this study. The electricity prices are obtained from a corporate price forecast for the Mid-C trading hub. The import and export electricity prices are differentiated by wheeling charges and transmission loss factors. The prices used are shown in Figure 12. The blue line in the figure shows the 24-hour moving average market price.  The typical price shape shows higher prices in the late summer (where air conditioning in California increases demand) and winter months (where heating in the Pacific Northwest increases demand). The prices are the lowest in the freshet period when water is abundant in the Pacific Northwest and temperatures are mild in the west coast. The natural gas prices were also obtained from corporate forecasts. Electricity market prices are highly correlated to natural gas prices, and they are used to determine the cost of running thermal plants. Market prices for the types of ASs considered in this thesis, namely the regulation up, regulation down and spinning reserves prices were obtained from CAISO_EXP region for the water year 2016. The CAISO_EXP region is an appropriate data source for this study given that BC Hydro participates in the California AS market. Figure 13 shows the AS prices used for this thesis. The AS prices are very volatile and exceed $100/MW for brief periods.   Figure 12 - Price data for the Mid-C trading hub (2035/36) -20020406080100Oct Nov Dec Jan Feb Mar Apr May Jun Jul Aug SepPrice ($/MWh)78   Figure 13 - AS prices (2016/2017) 4.2.4 Wind Data The data for this study was obtained from DNV Global Energy Concepts, Inc. (DNV-GEC), which was contracted by BC Hydro to provide synthetic wind data for several hypothetical wind farms across four BC regions with high potential for wind power development. DNV-GEC subcontracted the data generating work to 3TIER, a mesoscale atmospheric modeling and forecasting services firm. 3TIER created the synthetic data using the Weather Research and Forecasting (WRF) model, a mesoscale numerical weather prediction (NWP) model. As inputs, the model took global weather archive data and high-resolution topographic data. The historical global weather archives, which are archived representations of the historic state of the atmosphere over the planet, were used to initialize the WRF model. A higher resolution representation of the terrain, land use, and vegetation, which is required to accurately model the individual wind farms, was gathered from several sources and fed to the WRF model (DNV-GEC 2009).  The output from the WRF model was validated and calibrated using actual wind speed measurements at several Environment Canada weather stations, independent power producer (IPP) meteorological towers and BC Hydro meteorological towers. Since the power output from this procedure produces unrealistically smooth wind power time series, a SCORE (statistical correction to output from a record Extension) process was used to statically correct this behavior. Using this methodology, a 10-year long wind data spanning between Sep 1998 and Aug 2008 was produced for 104 hypothetical 020406080100Oct Nov Dec Jan Feb Mar Apr May Jun Jul Aug SepPrice ($/MW)Reg Up Prices Reg Down PricesSpin Prices 24 per. Mov. Avg. (Reg Up Prices)24 per. Mov. Avg. (Reg Down Prices) 24 per. Mov. Avg. (Spin Prices)79  wind farms in BC at a 10-minute time resolution. This was the primary output provided by DNV-GEC, and it is treated in this thesis as the actual wind generation.  For assessing the DA wind forecast error, DNV-GEC also provided 3-year long NWP forecast data spanning between Jan 2005 – Dec 2007 for each wind farm at an hourly resolution. To produce this wind forecast data, 3TIER ran another WRF model that was configured to provide once-a-day forecasts with a forecast horizon of 96 hours. The forecast simulations were initialized using archived global forecast grids (DNV-GEC 2009).  The NWP forecast is expected to give a level of accuracy and WFE statistical characterization that is comparable to the accuracy and WFE characterization of the actual wind forecasts used by BC Hydro.  Several hypothetical wind farms were selected to arrive at 15% wind penetration level, which means that the wind farms were added until the total installed wind capacity reached about 15% of the forecasted peak BC Hydro load for the study year 2035. The subset of selected hypothetical 18 wind farms 13 of which are located in the BC’s Peace Region, 3 in the Southern Interior, 1 in the North Coast and 1 in Vancouver Island. Some of these wind farms are currently operating in the BC Hydro system as discussed in 2.1. The wind farms were selected based on their UECs. This selection approach leads to limited geographic diversity, with 76% of the wind energy being generated by the 13 wind farms in the Peace region. The NWP forecast provided by DNV-GEC is for 96 hours initiated at midnight, so the DA forecast was taken to be the forecast for hour-ending 25 to hour-ending 48. A 3-year 10-minute wind data for the period Jan 2005 to Dec 2007 was extracted from the 10-year synthetic data discussed above and were converted into hourly averages to get the ‘observed’ data. This process resulted in 26280 data points for both forecasts and observations. In this thesis, it was decided to use the aggregate wind power forecast instead of the wind speed forecasts for the individual wind farms. This is more practical given the number and variety of wind farms considered in this study. Using the wind speed forecasts would be less practical since the 18 wind farms have to be considered individually while accounting for the spatial and temporal autocorrelation between them. To avoid this issue, the wind generation forecast provided by 3TIER for each wind farm was gathered and added to get one aggregate wind generation forecast. 3TIER derived the wind power forecasts from the wind speed forecasts by using the appropriate power curves that considered the turbine characteristics and the weather variables at each wind farm. 80  Analyzing the aggregate wind power forecast makes the analysis simpler, more practical, and more holistic. This also enables the development of better MC model as discussed by Papaefthymiou and Klockl (2008).   4.2.5 Forebay and Storage Targets The starting and ending forebay elevations are derived from HYSIM runs that consider multiple water years. Other boundary conditions and operations constraints, such as the Columbia River Treaty operation rule curves, were also taken from HYSIM.  4.2.6 Optimized Plants GOM usually optimizes the five major hydropower plants: Mica (MCA), Revelstoke (REV), Arrow Lakes (ARD), GM Shrum (GMS), and Peace Canyon (PCN). These five plants currently generate more than 72% of BC Hydro’s annual energy generation. The study year for this thesis is beyond the completion date of the Site-C dam, so Site-C (STC) is optimized together with the five major plants. The rest of the generation in the system is considered a fixed generation and therefore is not optimized. Some of the properties of the large hydropower plants included in this case study are given in Table 3. Table 3 - Some of BC Hydro's generation resources modeled Facility Reservoir name Live Storage (km3) Capacity (MW) Avg. Annual Gen (GWh) WAC Bennett dam - GM Shrum plant Williston 39.4 3000 13,230 Peace Canyon Dinosaur <1.0      700 3,110 Site-C Site-C <1.0 1,100 5,100 Mica Kinbasket 14.8 2,832 7,000 Revelstoke Revelstoke 1.85 3100 7,820 Keenleyside Arrow Lake 8.77 192 2,480  In addition to providing energy, the large hydropower plants are also capable of supplying AS. The AS sold in the market are divided into regulation and spinning reserves. Regulation reserves, which are 81  required to balance the minute to minute variations in load and generation, can only be provided by generation units equipped with Automatic Generation Control (AGC). The plants that were considered to provide regulation reserves were GMS, MCA, and REV. Spinning reserves, which manage imbalances created in longer time span than 10 minutes, can be provided by any of the optimized plants. 4.2.7 Reserves The wind reserves calculated according to the methodology discussed in Section 3.7.5 are given in Figure 14. The graph shows that the DA reserves are orders of magnitude higher than the within-the-hour reserves.   Figure 14 - Wind reserves 4.3 Other Inputs and Assumptions Some of the other inputs into the optimization model and their sources are presented in Table 4. The inputs listed in the table were kept the same for all optimization models to enable a fair comparison.    02004006008001000120014001600  Oct   Nov   Dec   Jan   Feb   Mar   Apr   May   Jun   Jul   Aug   SepReserves (MW)DA Up Reserves DA Down ReservesWithin-the-Hour Up Reserves Within-the-Hour Down Reserves82  Table 4 – Other inputs and assumptions that went into GOM Input Type Description (Data Source) Fixed generation Taken from HYSIM run that considers historical data. These fixed generations include small hydropower plants, Seven Mile generating plant, thermal plants, and ROR plants IPP generation The historical and contracted energy generation profiles were used to represent the generation from various IPPs Outage schedules The default planned outage schedules for the various units were used Operating reserve obligations The default value of 0.14 was used. This is used to account for load forecast error, transmission line outages, and major unit outages Transmission limits Current tie limits used by hydro used for all the models. It ranges from 1575 to 2350 MW for the BC-US intertie and from 300 to 500 MW for the BC-Alberta interties    Flow limits Physical and environmental restrictions that influence the plant’s generation capacity and spill releases were considered Generation limits Limited based on physical generator limits, number of units available, flow limits, reservoir storage, etc.  Market depth Both energy and AS market was limited by the available transmission capacity Load shedding penalty A per unit load shedding penalty of $10000/MWh was imposed to discourage load shedding   4.4 Cases Investigated Several model runs were conducted to estimate the stochastic savings and their sensitivity to input assumptions. The normal operating condition was assumed to be the water year 1968 (i.e., the average water year) coupled with the wind year 2006. Another two cases with dry and wet water years, coupled with wind years 2005 and 2006 were investigated to assess the effect of the annual water and wind variability on the wind integration cost and the stochastic savings. The drawback of coupling wind and water years that do not match is that it misrepresents the correlation between 83  the two weather parameters. This representation might affect the estimation of the wind integration cost and the stochastic savings. The combinations of water, wind and load year investigated for this thesis are shown below. Table 5 - Water, wind and load year combination  Water Year Wind Year Load Year Normal water year 1968 2006 2035 Dry year 1969 2005 2035 Wet year 1966 2006 2035  The cases investigated in this thesis are shown below in Figure 15. In the figure, the black boxes show the model runs used for the normal operating condition case while the grey boxes show the runs used for the sensitivity analysis. The sensitivity analysis for the water year was examined without changing the assumptions regarding the load shedding penalty cost or the AS market depth assumption. The sensitivity analysis for the load shedding penalty cost and the AS market depth assumption was conducted for the normal year only. An additional model was also run for the dry year with a low load shedding penalty cost assumption.  Figure 15 - Cases Investigated 84   Results and Discussion In this chapter, the results of the case study are presented and discussed. The first section presents the statistical characteristics of the WFE. The second section presents the results of scenario generation using the MC-based model. The third section discusses the results of scenario generation using the ARMA-based model. The fourth section presents the results of the evaluation of the scenarios generated by the MC and the ARMA based models. The fifth section presents the results of the scenario reduction process. The sixth section presents the results of the SO, where the stochastic savings, DA wind integration cost, load shedding and wind power curtailment are discussed. Finally, the last section discusses the sensitivity of the stochastic saving to varying input data and modeling assumptions.  5.1 Day-Ahead Wind Forecast Error Properties The DA WFE was calculated from the data discussed in Section 4.2.4 using Equation 1. For convenience, the WFE equation is repeated here.  𝑊𝐹𝐸𝑡 =  𝐷𝐴 𝑊𝑖𝑛𝑑𝐺𝑒𝑛𝐹𝑟𝑐𝑡 − 𝐴𝑐𝑡𝑊𝑖𝑛𝑑𝐺𝑒𝑛𝑡   𝑀𝑎𝑥𝑖𝑚𝑢𝑚 𝑊𝑖𝑛𝑑 𝐺𝑒𝑛 1 According to the definition given in Equation 1, a negative WFE is higher than expected wind generation that needs to be balanced by ramping down dispatchable hydropower generation and a positive WFE is lower than expected wind generation that needs to be balanced by ramping up dispatchable hydropower generation. The entire WFE time series is shown in Figure 16.  Figure 16 – The DA wind forecast error time series -0.8-0.6-0.4-0.200.20.40.60.80 5000 10000 15000 20000 25000WFEHour85  5.1.1 Statistical Moments Figure 16 shows the fluctuation of the WFE, which is close to the zero most of the time. However, its mean value is -2.37%, not 0. This indicates that the forecast is biased and that on average, the forecast is lower than the observed wind. Table 6 lists some of the WFE statistics (the WFE is written as a percentage in the table for convenience).  Table 6 – WFE statistics  MAE Mean Std Min Max Kurtosis Value 9.77% -2.37% 13.31% -63.27% 70.64% 4.88  Table 6 shows that the absolute mean average hourly WFE is 9.77% and that the WFE values range from -63.27% to 70.64%. These show that the WFE in BC is significant in the DA time frame since a WFE equal to about 10% of the installed capacity is expected for each hour and that a WFE up to 70% of the installed capacity can be witnessed. The standard deviation of 13.31%, which is higher than the MAE, suggests that the WFE can vary considerably. The Kurtosis value of 4.88 indicates that the WFE is a leptokurtic distribution, which implies that most of the variance comes from few large errors (outliers) instead of many small errors.  5.1.2 Stationarity and Seasonality  The dependence of the WFE time series on the forecast lead time is investigated following the methodology discussed in Section 3.1.2. Figure 17 and Figure 18 show the MAE and std for the forecast lead times from 25 to 48 hours. The figures show some increase in the MAE and std as the lead time increases. However, the increase is not substantial.  The average WFE for the last five hours is only about 9% higher than the WFE of the first five hours. This suggests that the forecast accuracy does not significantly worsen past the 24-hr forecast lead time. Therefore, the forecast lead time was not used as an input in the scenario generation process. Including the forecast lead time as an input would complicate the model without adding much value. Seasonality of the WFE was also investigated and summarized in Figure 19 and Figure 20. The figures show that the forecasted wind generation (and the actual generation), varies significantly by season. It can be seen that the average wind power generation in the winter is more than twice of that in the 86  summer. However, the various seasons have comparable levels of average WFE as seen in Figure 20, which indicates that the seasonality of the MAE is not that significant. Similar to the forecast lead time, adding the season of the year in the scenario generation process will complicate the model without adding much value. Therefore, the season of the year was not taken as input into the scenario generation process.  In conclusion, the WFE was considered to be independent of both the forecast lead time and the season of the year.  Figure 17 – Hourly MAE   Figure 18 – Hourly std of the WFE   Figure 19 - Monthly average wind generation  Figure 20 - Monthly MAE 00.020.040.060.080.10.1225 27 29 31 33 35 37 39 41 43 45 47MAEHour (hr ending)00.020.040.060.080.10.120.140.1625 27 29 31 33 35 37 39 41 43 45 47StdHour (hr ending)00.10.20.30.40.50.6JanFebMarAprMay JunJulAugSepOctNovDecWind ForecastMonth00.020.040.060.080.10.12JanFebMarAprMay JunJulAugSepOctNovDecMAEMonth87  5.1.3 Probability Distribution The probability distribution of the WFE was first visualized using the histogram given in Figure 21. The figure shows that the WFE values between -0.02 and 0.02 occur most frequently and that frequency of occurrence decreases as the absolute value of the WFE increases.  The scatter plot given in Figure 22 shows the joint probability distribution of the forecast level and the WFE. In this thesis, a scatter plot is used to show the joint probability distribution instead of a 3D-histogram because it was found to better for visualization. The figure shows that as many researchers like Bludszuweit et al. (2008), Bruninx, (2016), Ma et al. (2013), and Tewari et al. (2011) have suggested, the forecast level is an important variable that influences the shape of the WFE distribution. It can be seen from the figure that the WFE variation is high for intermediate forecast (i.e., between 25% and 80% forecast percentiles) but low for high and medium forecasts. This property was further investigated and shown in Figure 23 and Figure 24 Figure 23 shows the MAE of the WFE varies significantly depending on wind forecast level. The figure suggests that on average, the WFE is highest for intermediate wind forecast levels. It shows that the expected WFE is usually low for very low and very high wind forecasts. Figure 24, which shows the standard deviation for each forecast bin also has a similar trend with WFE variation being the highest for intermediate forecasts levels and lowest for low and high forecast levels. This can be partially explained by the nonlinear property of the wind power curve. As Tewari et al. (2011) pointed out, the nonlinear steep slope of the wind power curve at intermediate wind speeds can magnify small errors in the wind speed forecast. At low and high wind speed levels (i.e., lower than the cut-in speed and above the rated output speed as seen in Figure 4), the power curve is horizontal, which means small errors in wind speed forecasts have no effect on the accuracy of the wind generation forecast. Attempts were made to fit parametric probability distributions to WFE data. Before fitting distributions, the WFE data was shifted by +0.65 (i.e., 0.65 was added to the WFE data) so that even the lowest WFE value becomes positive. This was done so that probability distributions that do not accept negative values, like the lognormal distribution, can be used in the fitting process. Figure 25 shows the marginal WFE histogram and the best fit distributions. It shows that the marginal distribution of the WFE is a leptokurtic distribution with marked peaks and slimmer shoulders. This 88  property makes most of the parametric distributions a poor fit since they are not able to reflect the high kurtosis of the WFE distribution. The Beta distribution suggested by some researchers as discussed in Section 2.3 was not found to be a good fit for the WFE data set used in this thesis. A comparison based on the AIC values showed that the parametric distributions that fit the WFE data the best are the t location-scale distribution, followed by the logistic, normal and Rician distributions.   Figure 21 - WFE histogram  Figure 22 – Scatter plot of the wind forecast and WFE  Figure 23 - MAE for each forecast bin  Figure 24 - Std for each forecast bin 00.020.040.060.080.10.120.140.165 15 25 35 45 55 65 75 85 95MAEWind Forecast Percentile00.020.040.060.080.10.120.140.160.185 15 25 35 45 55 65 75 85 95StdWind Forecast Percentile89  As discussed above, the forecast level influences the distribution of the WFE. For example, a WFE of 0.2 is not possible when the forecast is only 0.1, and WFE of -0.2 is not possible when the forecast is 0.9. This is because it is impossible to have wind generation less than 0 or higher than the installed capacity. These bounds on the WFE cause it to have different distribution depending on the forecast levels. To show this property, the forecast was divided into 20 bins in increments of 0.05, and conditional histograms along with their best fit distributions were plotted. Figure 26, Figure 27, and Figure 28 show the conditional histograms and best fit distributions for typical low, medium, and high forecast bins. The figures show that the conditional WFE distributions vary significantly based on the forecast bins. The distribution is skewed left for low forecasts (e.g., Figure 26) and skewed right for high forecasts (e.g., Figure 28). The conditional PDF for the WFE of intermediate forecasts (e.g., Figure 27) is more symmetrical and resembles that of the normal distribution. The different shape of the distribution based on the forecast level indicates that the forecast level needs to be considered for an accurate statistical representation of the WFE.  Figure 25 – WFE marginal histogram and best fit distributions  Figure 26 – WFE conditional histogram and best fit distributions for the second forecast bin [0.05,0.1) 90   Figure 27 - WFE histogram and best fit distributions for the 10th forecast bin [0.45,0.5)  Figure 28 - WFE histogram and best fit distributions for the 19th forecast bin [0.9,0.95) In summary, it has been shown that the MAE and WFE variance are high for intermediate forecasts while they are low for high and low forecasts. It has also been shown that the distribution of the WFE changes depending on the forecast level. A proper scenario generation algorithm should have outputs that reflect this property. Attempts of distribution fitting have also shown that finding best fit parametric distributions for the WFE is challenging due to the WFE’s high kurtosis. 5.1.4 Autocorrelation and Partial Autocorrelation Functions The autocorrelation function (ACF) plotted in Figure 29 shows that the ACC decreases exponentially, although there is a slight uptick around 24 hours. The plot shows that the WFE is highly autocorrelated, suggesting that there is a high level of persistence.  A partial autocorrelation diagram was used to isolate the order of the autoregressive process. The PACF plot, Figure 30, shows that the PACC is statistically significant for the first 2 hours. This suggests that the WFE is a second-order autoregressive stochastic process. The negative PACC at lag 2 shows that there is some anti-persistence property incorporated into the forecast. This means that if the WFE is likely to reverse sign between t and t+2. This property reduces the overall persistence in the WFE. 91   Figure 29 - Sample ACF of the WFE time series  Figure 30 - Sample PACF of the WFE time series  5.2 Scenario Generation Using Markov Chains The goal of the scenario generation is to create a set of probabilistic scenarios that properly capture the properties of the WFE presented in section 5.1. The forecast lead time and season of the year will not be considered in the scenario generation process because they were found not to affect the WFE significantly as discussed in section 5.1.2. The conditional and joint probability distributions have also shown that the WFE is dependent on the forecast level, which suggests that the wind forecast should be used as an exogenous variable. Therefore, the transition probabilities will be conditioned on the forecast level. As seen in Figure 30 and discussed in Section 5.1, the PACF for the WFE is significant for the first 2 lags, showing the presence of a second-order autoregressive process. Thus, the proposed MC-based model in this study is based on a second-order MC. Therefore, in summary, the MC based model used in this thesis is a second-order MC with a wind forecast exogenous variable. After determining the fundamental properties of the MC model, the next step is selection of the number of states and state size. As discussed in Section 3.2, it is desirable to have the least number of states that will provide satisfactory results. This is especially important here since the model under consideration is a second-order MC-based model with an exogenous variable whose number of parameters increases rapidly with the number of states considered (i.e., the number of parameters in 92  the proposed model are n3 * f, where n is the number of WFE states, and f is the number of wind forecast states).  As discussed in Section 3.2, the appropriate number of states was chosen by discretizing the WFE into various numbers of states and by comparing the properties of the resulting time series. This is done for states numbers of 11, 20, 32 and 82. The states were created by using Equations 11 and 12. The comparison of the ACF and PACF for various levels of discretization is shown in Figure 31 and Figure 32.  Figure 31 - ACF comparison   Figure 32 -PACF comparison Figure 31 suggests that the ACF of the WFE is not considerably affected by the discretization. However, Figure 32 shows that the PACF is affected by discretization. The PACF of discretized WFE time series is different from that of the continuous WFE time series, especially when the WFE is discretized into a smaller number of states. The lag-2 PACC of the discretized WFE time series is much lower than the continuous one for the 11-state discretization but is almost the same as the continuous one for the 82-state discretization. However, the lag-3 and lag-4 PACCs of the 11-state discretization are higher than that of the continuous time series. This suggests that the anti-persistence property of the continuous WFE time series (i.e., negative lag-2 PACC) has been distributed into lags of 2, 3, and 4 hours for the discretized time series. This has the effect of slightly lowering the ACC of the 11-state WFE time series as seen in Figure 31. Yet, given the small impact of discretization on the ACF and the fact that the number of parameters that have to be estimated is 93  much lower for the 11-states discretization (11*11*11*f =1331f) than even the 20-state discretization (20*20*20*f = 8000f), 11 states are deemed sufficient for this study. The intervals that defined the 11 states were discretized by 0.1 except for the three intervals that reside in the middle (i.e., between -0.1 and 0.1) where the frequency of the WFEs is higher and finer resolution was desired. The column headers in Table 7 show the selected intervals and their corresponding states.  A similar discretization process was applied to the wind forecast to create forecast states. As a compromise between model accuracy and parsimony, a total of 10 forecast states were found to be ideal for this study. The ranges used for the wind forecast state discretization were values ranging from 0 to 1 at increments of 0.1. This means that there are a total of 13310 (i.e., 1331*10) parameters in the proposed model. Most of the parameters (i.e., the transition probabilities) are zero since some state transitions are almost impossible. For example, the probability that the WFE will transition from WFE state 4 at t-1 to WFE state 1 at time t is zero. Since the model treats such unobserved events as zero, estimating 13310 parameters from the available data does not make the accuracy of the estimates very questionable. After deciding on the order of the MC, the exogenous variable to use, and the number of states to consider, the transition probability matrixes are estimated following the steps outlined in Section 3.2. The state transition probability Pjkf,i  which is defined for the proposed MC-based model as P(Εt = i | Εt-1 = j, Εt-2 = k, Φt = f), where i, j, and k represent the WFE states and range from 1 to 11 while f represents the forecast states and ranges from 1 to 10. The transition probability Pjkf,i is estimated by counting the number of transitions to a specific state i given the three conditions (Εt-1, Εt-2, Φt) and dividing it to the total number of transitions to all states (Ε t = 1 to 11) given the same three conditions. For example, the number of transitions to state 6 given the three conditions Εt-1 = 5, Εt-2 = 6, Φt = 2 is 32 while the total number of transitions with these three conditions to all states, including state 6, is 159. Therefore, the probability that the WFE will transition to state 6 given the three conditions is 32/159=20.13%. The state transition probability matrix is constructed by repeating this calculation for every combination of i, j, k, and f.  The resulting transition probability is a 4-dimensional, 11 x 11 x 11 x 10 matrix which is difficult to illustrate. Hence, the 2-dimensional transition probability matrixes Pj,i, Pk,i, and Pf,i are given in Table 94  7, Table 8, and Table 9 respectively. The matrixes show the individual effects that the three conditions Ε t-1, Ε t-2, Φt have on Ε t.  The Pj,i transition matrix is given in Table 7. In the table, the rows show the starting states, the columns show the ending states, and the number in the intersecting cells shows the probability of transition from the starting state to the ending state. For example, if the starting state (WFE state at t-1) is state 4, there is 67.3% chance that the ending state (the WFE state at time t) is state 4, 19.9% chance that the ending state is state 5, and so on. It is worth noting that the sum of probabilities across the rows equals 1. The table shows that the diagonal elements of the matrix have the highest value for all the starting states, indicating that the WFE at time t is most likely to stay at the same state as the WFE at t-1. This is true even more so for the starting states 1, 6 and 11, which represent the states with the most negative, lowest, and most positive WFEs respectively. Another notable property is that state number 5 is the only state that can be transitioned to from all starting WFE states at t-1. Table 8 shows that at time t, the probability of staying at the WFE state of t-2 is lower than the probability of staying at the WFE state of t-1. This table also shows that it is possible to start at any state and end up in states 4 to 7 after two hours. The anti-persistence property of the WFE can also be seen in Table 8. For example, P9,9 decreased from 0.597 in Table 7 to 0.399 in Table 8 while P9,8 increased from to 0.240 in Table 7 to 0.324 in Table 8. This shows that the error has been reduced since state 8 has less error than state 9. Table 9 shows the probability of the ending WFE state for a given forecast state. The table shows that WFE is dependent on the forecast. The states with the highest and lowest forecasts (i.e., starting states 1 and 10) have the lowest WFE, which can be deduced from the fact that their ending states are more likely to be state 6. The forecast states 3-7 have the broadest WFE state range, likely to end up in more than 5 WFE states with similar probabilities. This table reflects the joint PDF of the wind forecast and WFE shown in Figure 22.   95  Table 7 - Pj,i transition probability matrix – ending WFE state given WFE state at t-1  Table 8 - Pk,i transition probability matrix– ending WFE state given WFE state at t-2  Table 9 - Pf,i transition probability matrix– ending WFE state given Forecast state at t  96  The joint state transition matrix Pjkfi combines the effect of the three transition probability matrices discussed above to capture the relevant properties observed in the WFE. The sampling process was carried out using the CSTPM and Equation 15. The equation is repeated here for convenience:   𝐶𝑗𝑘𝑓,i =  ∑ 𝑃𝑗𝑘𝑓,m𝑖m=1 15 for j, k, n = {1, 2, …, 11} are the WFE states and f = {1, 2, …, 12} is the wind forecast state. A sample of the CSTPM Cj,i that corresponds to Pj,I (Table 7) is shown in Table 10. Table 10 - Cumulative transition probability matrix Cj,i  As described in Section 3.2, the inverse transform method was used to sample WFE values from the CSTPM. As a simple illustration of the sampling process, consider Table 10. Given Εt-1 is 5 and the sampled random number is 0.67, the ending state Εt is 5 since 0.67 is between 0.14 and 0.79. The values on row 5, which correspond to the starting state Ε t-1 of 5, ensure that for the starting state 5, 1%(0.01-0), 13% (0.14-0.01), 65% (0.79-0.14), 21% (0.98-0.79) and 2%(1-0.98) of the ending states Ε t will be states 3, 4, 5, 6, and 7 respectively. This also makes it impossible to start from the starting state 5 and end at states 1, 2, 8, 9, 10 and 11. This ensures that the probability distribution of the WFE is preserved.  To create N number of scenarios for a single day, random numbers with 24 rows and N columns are first sampled. These random numbers are then transferred into WFE states by using the CSTPM. For 97  the first forecast hour, Εt-1 and Εt-2 are not available. This is because the forecast is issued every 24 hours and the WFE error state from the final hours of the previous day are not applicable. Therefore, Εt-1 and Εt-2 are not used to estimate Ε1, which means that the only information used when creating Ε1 is Φ1. Thus, a CSTPM produced from Pfi is used to sample the first hour of the day. Similarly, for the second hour (Ε2), the Εt-2 information is not available. Therefore, a CSTPM matrix produced from Pjf,i (excluding k) is used for sampling. For the rest of the hours, the complete CSTPM matrix is used since all the required information to calculate Pjkf,i is available.  The outputs of the simulation, the WFE state time series, is converted to WFE values by replacing the states with the mean WFE values of the states. For example, WFE state 9, which contain WFE values that range from 0.2 to 0.3, is converted to the WFE value of 0.242 (0.242 is the mean value for WFEs that reside in state 9). The other popular approach of converting states to values, namely the sampling of uniformly distributed random values from within the states (for example, sampling between 0.2 and 0.3 for state 9) was also tested. This approach has the advantage of giving continuous WFE values as outputs but the scenarios it generated were found to be less representative of the observed data as compared to the chosen approach. The number of scenarios needed to capture the WFE distribution adequately is determined in accordance with the methodology discussed in Section 3.2. For this process, 2000 scenarios for a single day were first generated. Then the mean and std of the scenario sets were calculated for the first n scenarios (n ranging from 1 to 2000) and plotted. This plot is then used to decide the number of adequate scenarios to sample. Figure 33 is one such plot for a typical day. It shows how the mean and std of the WFE evolve for the 24 hours of the day when the number of scenarios considered in the calculation of the statistical moments is increased. In the figure, the computed sample mean and std vary as the number of scenarios increases and they start to stabilize and approach the population mean and std after hundreds of scenarios are sampled. From the plot, 1000 scenarios were deemed to be adequate to sufficiently capture the WFE distribution since the statistical parameters do not vary significantly beyond this point. The fact that the WFE stabilizes at different mean and standard deviation values for the different hours of the day shows that all forecast levels have their own unique statistical moments. Similar figures plotted for days with other forecast levels and generation characteristics also showed that 1000 samples are adequate. 98  Using MATLAB to generate the 1000 scenarios on a personal computer with 8GB of RAM and 7th generation Core i5 processor for one single day takes less than 0.6 seconds on average. This indicates that this scenario generation approach has a low computational burden.   Figure 33 - Mean (bottom half) and std (upper half) computed after considering n number of scenarios for a typical day (MC-based model). The various lines represent the 24 hours in a day.   5.3 Scenario Generation Using ARMA The second type of model proposed in this thesis WFE is an ARMA based model. ARMA models take a different approach to modeling the WFE time series. No state discretization is needed, and the calculation of state transitions probabilities is also not required. All the properties of the stochastic process are captured in few model parameters. As mentioned in Section 3.3, An ARMA(p,q) process 𝑋𝑡  is expressed as:  𝑋𝑡 = 𝑐 + ∑ 𝜙𝑖𝑋𝑡−𝑖 + ∑ 𝜃𝑖𝜀𝑡−𝑖𝑞𝑖=1+ 𝜀𝑡𝑝𝑖=1 17 where 𝑐 is the constant term; 𝜙𝑖 are the autoregressive parameters; 𝜃𝑖 are the moving average parameters, and 𝜀𝑡  is the white noise term.  99  The inputs needed for the ARMA model include the p and q parameters, and the CDF that is used for probability distribution transformation. Initial values for these inputs are obtained from the analysis of the WFE conducted in sections 5.1.3 and 5.1.4. The initial ARMA parameters p and q (i.e., the number of time steps to consider) are determined from the ACF and PACF plots shown in Figure 29 and Figure 30. The statistically significant PACC for the first two lags suggest a second-order autoregressive process and thus, a p value of 2 is used. The nearly exponential drop in ACF indicates a first or second order moving average process and that a q value of 1 or 2 can be used. As mentioned in Section 3.3, ARMA models assume that the input data is normally distributed. Therefore, to accurately model a time series using an ARMA model, a probability distribution transformation is necessary. This transformation requires the CDF of the input data. Since it was concluded in Section 5.1.3 that no parametric probability distribution can accurately capture the properties of the marginal WFE distribution, the empirical CDF (ECDF) is used. Using the ECDF ensures that the WFE is represented accurately and that any erroneous assumptions are not made about the WFE distribution. This is in line with the methodology used by Ma et al. (2013). The marginal distribution used for transforming the WFE is shown in Figure 34. Once the WFE is transformed into a normally distributed data using the methodology described in Section 3.3, it is fit into an ARMA model. To ensure that the most appropriate numbers of parameters are selected, different ARMA models with varying combinations of p and q were fitted, and their goodness-of-fit was compared based on AIC and -log likelihood values. ARIMA models, which assume the WFE to be nonstationary, were also considered. The results of the AIC and -log likelihood comparisons are shown in Table 11 below.  Table 11 - AIC and -log likelihood comparisons ARMA Order AIC values -Log Likelihood ARMA(1,1) 16063.03 8027.52 ARMA(1,2) 15801.74 7895.87 ARMA(2,1) 15685.05 7837.52 ARMA(2,2) 15666.52 7827.26 ARMA(2,3) 15668.51 7827.26 ARMA(3,2) 15668.55 7827.27 ARMA(3,3) 15669.47 7826.73 100  ARMA Order AIC values -Log Likelihood ARIMA(1,1,1) 17156.77 8575.38 ARIMA(2,1,1) 15718.32   7855.16 ARIMA(2,1,2) 15672.84 7831.42  The selection was also verified by using the auto.arima function in the ‘forecast’ package (Hyndman 2017) for the R programming language. The function optimizes the number of AR and MA parameters to be used. In addition, the function investigates whether the inclusion of seasonal components or nonstationary improves the model. The function also indicated that ARMA(2,2) process is ideal to represent the WFE. The function also confirms that the WFE can be considered a stationary stochastic process.  The parameters estimated from the fitting process were: 𝜙1 = 1.5133, 𝜙2= -0.5678, 𝜃1 = -0.3026, 𝜃2 = -0.0634, 𝑐 = 0, and σ2  =  0.1062. That means, the WFE ARMA process can be written as: 𝑋𝑡 = 1.5133 ∗ 𝑋𝑡−1 − 0.5678 ∗ 𝑋𝑡−2  − 0.148 ∗  𝜀𝑡−1 − 0.0634 ∗ 𝜀𝑡−2  + 𝑁(0,0.1062) where 𝜀𝑡−1 and 𝜀𝑡−1 are the white noise terms and N(0,0.1062) represents a white noise with a mean of 0 and variance of 0.1062. Once the parameters are determined, the ARMA model is used to simulate a stochastic time series that is later to be transformed to WFE values.  As mentioned in Section 3.3, the output of the univariate ARMA model simulation is a normally distributed time series with zero mean and constant variance. It is also not conditioned on the forecast level. To ensure that that the probability distribution of the WFE is preserved and that the WFE statistics is conditioned on the forecast level, the time series that is simulated from the ARMA model is reverse transformed into WFE time series by using 20 conditional ECDFs. The 20 conditional ECDFs approximate the joint probability distribution of the wind forecast and the WFE, which ensures that the relationship between the WFE and the forecast are reflected in the scenario generation process. The different conditional ECDFs that are used for the reverse transformation are shown in Figure 34.  101  The 20 conditional ECDFs in Figure 34 all have different shapes which ensure that the WFE is conditioned on the forecast. All ECDFs have different statistical parameters like mean and variance. For example, bin 2, which represents wind forecast between 5 and 10% of the installed capacity, has a mean WFE of -0.042, and a relatively low standard deviation of 0.083 while bin 10 (i.e., forecast level between 45% and 50% of capacity) has a mean WFE of -0.017 and a relatively high standard deviation of 0.167. Meanwhile, bin 20 (i.e., forecast level between 95% and 100% of capacity) has a mean WFE of 0.04 and a relatively low standard deviation (0.074). These values are reflective of what is seen in Figure 22. The ECDF for bin 1 ensures that most of the WFEs are positive and the one for bin 20 ensures that the WFEs are mostly negative. Similarly, the properties of the WFE are accurately represented for the other forecast levels. One of the drawbacks of using ECDFs to reverse transform the outputs of the ARMA model into WFE values is that it bounds the WFE values by the minimum and maximum observed WFEs. A multivariate ARMA (i.e., MARMA) model that considers both the WFE and the forecast was also tested, but its outputs were less successful in capturing the joint distribution as compared to the approach used in this thesis.    Figure 34 – The marginal ECDF used for transforming the WFE to standard normal and the various conditional ECDFs used for reverse transforming normally distributed ARMA output back to WFE 102  The disadvantage of using 20 different ECDFs instead of the single marginal ECDF for the reverse transformation is that it lowers the ACF of the simulated WFE time series. This is because the ACF information contained in the untransformed time series is disturbed when carrying out reverse transformation with 20 ECDFs that have different characteristics. The small loss of ACF was the reason the single marginal ECDF was used instead of the 20 ECDFs to transform the WFE into normal distribution prior to model fitting. The univariate ARMA model only requires WFE data as an input so the forecast was not used as input. Therefore, only a single ECDF was used for the transformation since it ensures that the maximum amount of autocorrelation information is preserved. The use of the single marginal ECDF for the transformation and the multiple conditional ECDFs for the reverse transformation does not cause many issues since both the marginal and conditional ECDFs represent the same joint probability distribution. This ARMA model can also be used to sample as many WFE scenarios as desired. The ‘arima.sim’ function in R is utilized in simulation to generate the required number of scenarios. The number of ARMA-based scenarios that are required to capture the underlying WFE distribution reliably is determined using the same approach used for the MC-based model. Figure 35 shows that the sample mean and std of the WFE scenarios approach the population mean and std as the number of sampled scenarios increase. From Figure 35 and other similar plots made for other days with different forecast characteristics, 1000 scenarios were found to be adequate to capture the underlying characteristics of the WFE distribution. Creating 1000 scenario in R using a personal computer with 8GB of RAM and a 7th generation Core i5 CPU takes 0.3 seconds on average, which means that the ARMA-based scenario generation algorithm also has a low computational burden.  103   Figure 35 - Number of scenarios vs statistical parameters: standard deviation (upper half) and mean (lower half) for the 24 hours of a typical medium wind forecast day (ARMA model) 5.4 Scenario Quality Evaluation Various approaches were used to assess the performance of the MC and ARMA based scenario generation models in capturing the relevant properties of the WFE. First, to visualize the scenarios and compare the similarity of the outputs, 1000 equiprobable 24-hour scenarios for a typical day were generated using both models. Figure 36 and Figure 37 show the plots of these scenarios.  Figure 36 shows 1000 scenarios generated using the MC-based model along with regions of scenario concentrations, represented by shades. The outputs of the MC-based model are limited to 11 different WFE values since only 11 WFE states were considered. This leads to overlapping when the scenario are plotted on a figure. To help with the visualization, darker shading is applied to the areas with higher number of scenarios.  Figure 37 shows 50 scenarios generated by the ARMA-based model. Only 50 scenarios were shown because showing all 1000 scenarios would cover the entire graph and hide the shadings showing the scenario concentrations. In both figures, the dashed blue line with dots is the observed WFE and the thick black line is the mean of the 1000 scenarios  104   Figure 36 - Scenarios generated using the MC-based model.   Figure 37 – 50 Scenarios generated using the ARMA-based model.  Some general similarities can be seen in Figure 36 and Figure 37. For example, the WFE shows high variation in hours 1 to 6 and hours 18 to 24. This is because the forecast is intermediate for these hours. For hours 7 to 15, the WFE is more likely to be positive and its variation is lower (most of the 105  scenarios are in a tight band as shown by the dark shading in Figure 36 and Figure 37). This is because the forecast is high for hours 7 to 15, which means the error variation is low and it is more likely to be positive rather than negative. These graphs also show the effectiveness of the models in capturing the uncertainty bands of the WFE. Plots made for different days also revealed that the uncertainty bands are represented in the scenarios generated using both the MC-based and the ARMA-based models.  5.4.1 Statistical Moments One of the ways the quality of the two scenario generation models was assessed was by comparing the statistical moments of their outputs with the statistical moments of the original WFE time series. For this assessment, one scenario with a time horizon of three years was first generated using both the MC and the ARMA models to enable apple to apple comparison with the available WFE data. Then, the statistical parameters of the scenarios generated from the two models were computed and compared with the statistical moments of the WFE.  Table 12 shows the comparison of the statistical moments.  Table 12 - Statistical moments comparison  Mean MAE Std Kurtosis Original WFE -2.37% 9.77% 13.31% 4.88 MC-simulated -2.38% 9.45% 13.10% 4.68 ARMA-simulated -2.34% 9.75% 13.26% 4.83  The table shows that both the MC and ARMA based models are able to satisfactorily reproduce the statistical moments of the WFE. Although both models performed well, it can be seen that the ARMA model performed slightly better in recreating the MAE, standard deviation, and kurtosis values of the WFE, 5.4.2 Marginal and Conditional Distributions The three-year-long scenarios generated in section 5.4.1 were also used to investigate the ability of the ARMA-based and MC-based models to create WFE scenarios that preserves the underlying joint 106  and marginal distributions of the WFE. The comparisons based on the scatter plots and ECDFs are shown in Figure 38, Figure 39, and Figure 40. Figure 38 shows the discrete joint probability distribution obtained from the MC-based model. Small randomness was introduced into the plot to avoid displaying the WFE as a single line.  Although the WFE values shown on the graph are discrete, the figure still resembles the joint distribution shown in Figure 22. This indicates that the MC-based model is able to approximate the joint probability distribution of the wind forecast and the WFE. Figure 39, which shows the joint distribution obtained from the ARMA-based model’s outputs, resembles Figure 22 even more. This indicates that the ARMA-based model is better than the MC-based model at approximating the joint probability distribution. This improved approximation is mainly because of the lack of WFE discretization in the ARMA-based model. Yet, ‘stepping’ or stair-like shape is observed in the maximum positive and negative WFE bounds of Figure 39. This ‘stepping’ is the result of the reverse transformation that was carried out using 20 ECDFs derived from 20 discrete forecast bins that have different maximum and minimum WFE bounds.  The figures also show that both the ARMA and MC based models were effective in capturing the relationship between the forecast and the WFE. Both figures show that there is a low probability of having positive WFE when the wind forecast is low or having negative WFE when the forecast is high. The maximum variation of WFE is seen for intermediate forecasts while the minimum WFE variation is seen for low and high forecasts.  Figure 40 shows the ECDFs of the MC and ARMA based model outputs along with the ECDF of the actual WFE. It shows that both models accurately capture the marginal WFE distribution. The ECDF of the MC simulation is ‘stepped’ because of the discretization, but it still closely follows the ECDF of the observed stochastic process. The observed and ARMA simulated ECDFs are very close to each other and are only observable after magnification. This closeness is the result of the reverse transformation that uses the conditional ECDFs derived from the actual marginal ECDF.  107   Figure 38 – Scatter plot of WFE obtained from the MC-based model  Figure 39 - Scatter plot of WFE obtained from the ARMA-based model  Figure 40 - Comparison between the marginal ECDFs of the observed WFE time series, MC-based model simulated time series, and ARMA-based model simulated time series. 108  5.4.3 ACF and PACF Figure 41 shows the autocorrelation function comparison. The ACF of the MC-based model’s outputs is higher than both the ACF of the ARMA based model’s outputs and the ACF of the original WFE time series. The reason is that, as will be discussed in the following section, the MC-based model underestimates the lag-2 PACC, which means that the anti-persistence property is not fully captured in its outputs. The ACF of the ARMA-based model’s outputs and the ACF of the actual WFE are close to each other for the first 16 hours but start to diverge after the 17th hour. All the ACFs are statistically significant for the first 24 hours. In general, it can be seen from Figure 41 that the ARMA based model is better at recreating the ACF of the WFE.  Figure 42 shows that the lag-1 PACCs of the three WFE time series are close to each other. At lag 2, the PACC for the MC simulated WFE is significantly lower than the other two. This result is one of the side effects of the discretization. The original WFE’s PACC decreased at lag 2 when it was discretized into 11 states as shown in Figure 32. However, when the original WFE was discretized into 11 states, the third to five-hour lags had a negative PACC that was higher than the actual, which means the overall ACF was preserved. The MC-based model was only able to recreate a negative PACC for lag-2, so the anti-persistent property of the WFE was not adequately reflected in the simulated time series. The slightly lower PACC values of the ARMA-based model outputs for the first two lags is the result of the probability transformation that utilized different conditional ECDFs. In summary, the PACF plot shows that the ARMA-based model is better than the MC-based model at capturing the autoregressive property of the WFE.  5.4.4 Stability  The stability of the proposed scenario generation models was investigated by assessing their ability to create scenarios with similar properties in simulations. This is first done by creating the three-year-long synthetic time series discussed in section 5.4.1 ten times and verifying that the outputs of the simulations have consistent statistical properties. This comparison showed that the statistical parameters are consistent with the ones presented in Table 7. For shorter time scale, the investigation was conducted by using both MC and ARMA based models to create 25 sets of scenarios that each contain 1000 scenarios and by examining their properties.  109   Figure 41 - Autocorrelation function comparison  Figure 42 - Partial autocorrelation function comparison The properties of the various scenario sets generated from both models were found to be stable as different runs resulted in scenario sets that have consistent statistical parameters. Figure 43 and Figure 110  44 show the means of 25 WFE scenario sets that contain 1000 WFE scenarios each. It can be seen in the figures that both the MC based and ARMA based models have ‘stability’ when creating scenarios since the means of the WFEs estimated from several scenario simulations were very close to each other. Figure 43 and Figure 44 were plotted for a single day, but similar plots made for days with different wind forecast characteristics also showed stability.   Figure 43 - Scenario stability analysis – means (lower half) and standard deviations (upper half) of 25 simulations from the MC-based model for a typical day    Figure 44 - Scenario stability analysis – means (lower half) and standard deviations (upper half) of 25 simulations from the ARMA-based model for a typical day   In summary, the scenario quality evaluation has shown that both the MC and ARMA based models are effective in capturing the statistical moments and the probability distributions of the WFE stochastic process. The ARMA-based model was found to be better at capturing the ACF and PACF of the WFE time series. It was also observed that both models are able to generate stable scenarios that do not vary significantly when simulations were repeated.  5.5 Scenario Reduction  As discussed in Chapter 2 and Chapter 3, it is desirable to use a small number of scenarios in a SO problem in order to minimize the computational burden of the problem. The purpose of scenario reduction is to find a balance between the amount of information captured in the scenarios and size of the optimization problem. The 1000 equiprobable scenarios generated using the ARMA and MC 111  based models are too numerous for the large-scale optimization model considered in this thesis, so a reduced set of representative scenarios have to be produced.   The scenario reduction algorithm presented in the methodology section is implemented in both MATLAB and R programming languages to work with the MC and ARMA based models. The compromise between the information retained and the reduced number of scenarios is investigated following the methodology found in Dupačová et al. (2003) and Follestad et al. (2011). Figure 45 shows the number of scenarios (in a log scale) versus the percentage of information retained as measured in terms of Kantorovich distance. In the figure, 0% represents the information contained in a single scenario (i.e., the deterministic representation) while 100% represents the information contained in all scenarios. Both curves suggest that most of the scenarios contain similar information and that few scenarios can capture a significant amount of information about the underlying distribution. For example, the graph shows that more than 35% of the information contained in the 1000 equiprobable scenarios generated from both MC and ARMA based models can be captured by using just 10 scenarios with different probabilities. The curve is shown for one typical day but the curves for other days also show comparable properties.   Figure 45- Relationship between number of scenarios and percentage of information captured 112  Figure 45 and other similar plots for other days were used to decide on the number of scenarios to retain after the scenario reduction process. 40 scenarios were deemed adequate to capture the necessary set of information to make a reasonable stochastic decision as 40 scenarios contain more than 50% of the information contained in the 1000 scenarios. The size of the proposed SO problem was also found to be computationally feasible with 40 scenarios. For demonstration purposes, 40 scenarios reduced from the original 1000 are given in Figure 46 and Figure 47 for the MC-based and the ARMA-based models respectively. The figures show the reduced versions of the scenario sets presented in Figure 36 and Figure 37 respectively. It can be seen that most of the variation seen in Figure 36 and Figure 37 is still preserved after scenario reduction. Note that the scenarios shown in Figure 46 and Figure 47 do not have equal probabilities. The line thickness indicates the likelihood of the scenario’s occurrence. The thicker the line in these figures, the higher the probability. It can be seen from the graphs that the scenarios near the zero WFE value have higher probabilities.  The means and standard deviations of the full scenario sets and the weighted means and standard deviations of the reduced sets were also compared to verify that the statistical characteristics are preserved after the scenario reduction process. Figure 48 shows one such comparison for a typical day. It can be seen that the means of the full scenario sets and the reduced sets are very close to each other. This suggests that the statistical parameters of the WFE are adequately preserved after scenario reduction.  Once the scenario generation and reduction processes are completed, the next step is to add the WFE scenarios to the wind forecast to get the wind generation scenarios. When adding the WFE to the forecast, there is a small chance that the sum could result in wind generation that is either negative or higher than the installed capacity. This is because the outputs of both scenario generation algorithms do not perfectly capture the joint probability distribution of the wind forecast and the WFE. For example, it is seen in Figure 38 and Figure 39 that both models can give a positive WFE when the forecast is 0, which is impossible since it leads to negative wind generation. To avoid this, wind generation lower than zero was truncated at zero and generation higher than the installed capacity was truncated by the installed capacity. The final DA wind generation scenarios obtained from this process are used in the stochastic optimization.  113   Figure 46 - 40 reduced scenarios from 1000 (MC). The line thickness indicates relative probability  Figure 47 - 40 reduced scenarios from 1000 (ARMA). The line thickness indicates relative probability 114   Figure 48 - Statistical parameter comparison between the full and reduced scenario sets for a typical day 5.6 Optimization Results The results of the optimization studies and the comparison between the different models in terms of optimal value, load shedding, and wind power curtailment are presented in the following subsections.  For convenience, the four types of models described in Section 3.8 are repeated below a. A deterministic model without DA reserves b. A deterministic model with DA reserves c. A stochastic model with scenarios generated from the MC-based model and without DA reserves d. A stochastic model with scenarios generated from the ARMA-based model and without DA reserves The four models will be referred to as models (a), (b), (c), and (d) from now on. 5.6.1 Assessment Using the Expected Optimal Values First, an assessment is conducted for the normal water by comparing the expected optimal values of models (a) to (d) as discussed in Section 3.8. This provides an estimate of the savings that can be expected from using SO under normal operating conditions. Figure 40 shows the difference between -0.2-0.100.10.20.31 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24WFEHourMean of All MC Scenarios Weighted Mean of Reduced MC ScenariosMean of All ARMA Scenarios Weighted Mean of Reduced ARMA ScenariosStd of All MC Scenarios Weighted Std of Reduced MC ScenariosStd of All ARMA Scenarios Weighted Std of Reduced Arma Scenarios115  the optimal values of the four models and the optimal value of model (b), which is used as a reference. The graph shows that there are days where model (a) gives lower optimal values than the other models. This happens because model (a) shows the results of trading in markets with the actual wind generation while the other models show the results from trading in the markets with the wind forecast. Therefore, the optimal values of model (a) might be lower than the optimal values of the other models in days where the wind generation is lower than the wind forecast. The figures also show that models (c) and (d) rarely result in lower optimal values than model (b) since they have more freedom to trade in the markets, which is the result of not having to keep DA reserves explicitly.  Figure 49 – Difference in expected optimal values (the optimal value of model (b) is the reference) A more comprehensive comparison is achieved by summing the daily optimal values for the year and comparing the sums in accordance with the methodology outlined in Section 3.8. To reiterate the methodology, DA wind integration cost (also called cost in short) was defined as the difference between the sum of the optimal values of models (b) and (a). According to this definition, model (a) has zero integration cost while model (b) has maximum integration cost. Using this approach, the cost of wind integration for this water year was found to be $5.21/MWh. Figure 50 shows a comparison of the four models based on the sum of their expected optimal values. The figure shows that in general, models (c) and (d) results in significantly lower costs representing only 24% and 17% of the cost of model (b) respectively.  -$1,500,000-$1,000,000-$500,000 $- $500,000 $1,000,000 $1,500,000 $2,000,0000 50 100 150 200 250 300 350Difference in Optimal Values ($CAD)DayDeterministic Without Reserves Deterministic with ReservesStochastic with MC Based Scenarios Stochastic with ARMA Based Scenarios116   Figure 50 – Relative cost comparison based on the expected optimal values This suggests that for a normal year, SO is expected to reduce the wind integration cost by up to 83%. However, it is emphasized that these are only the expected savings, which assume that the DA reserves are adequate to manage all WFEs and that the WFE is perfectly represented by the wind generation scenarios sets. Therefore, these results may not reflect the performance of the models in real-world situations. The difference between the expected and actual performance can be significant, for example, when the WFE is very high and the amount of DA reserve is not sufficient to manage the WFE or when an extreme WFE is not captured by any of the WFE scenarios. Under these circumstances, a load shedding event is likely to occur which can significantly reduce the actual optimal values and the expected gains. Therefore, an ‘out-of-sample assessment’ as discussed in Conejo et al. (2010) is conducted to evaluate the performance of the stochastic models in real-world situations. 5.6.2 Out-of-Sample Assessment In this assessment, the DA market decisions obtained from the stochastic runs are fixed, and their performance is evaluated by running RT operations with the fixed DA decisions and the actual wind generation. The optimal values from this RT runs are then compared in a similar way as Section 5.6.1. Figure 51 shows the actual optimal values of the four types of models for the normal water year. There are significant deviations in optimal values in this figure that are not seen in Figure 49 (the Y-axis in Figure 51 is capped at ±1,000,000 to help with visibility). These deviations are the results of load shedding events. Figure 51 also shows that not all the models shed load on the same day.  Since 0%10%20%30%40%50%60%70%80%90%100%Deterministicwith ReservesStochastic withMC BasedScenariosStochastic withARMA BasedScenariosRelative Cost117  all models are affected by the WFE in the out-of-sample assessment, model (a) always performs as good as or better than the other models. This shows the value of perfect information.   Figure 51 - Difference in the daily actual optimal values When comparing the sum of the actual optimal values, it was found that the cost of DA wind integration has increased to $6.43/MWh. The increase in cost was the result of several load shedding events. Figure 52 shows that the costs of models (c) and (d) are now 30% and 52% of the cost of model (b) respectively. These equate to 70% stochastic saving for model (c) and 48% stochastic saving for model (d). While these cost savings are significant, they are less than what was found in the assessment that used the expected optimal values. Figure 53 shows how well the different components of the objective functions of the stochastic models (c) and (d) perform when compared to model (b). This comparison is helpful in determining the sources of the stochastic saving. The figure shows that most of the improvements of the stochastic models come from reduced DA imports, increased RT exports, increased RT AS transactions and reduced load shedding. The trading benefits of the stochastic models (c) and (d) as compared to model (b) are mostly the result of the DA capacity reserves that have to be kept for model (b). The reserves reduce the flexibility of the system and limits BC Hydro’s ability to participate in the energy and AS markets. The reduced load shedding in the stochastic models is the result of a better representation of the wind forecast uncertainty through the WFE scenarios. The main disadvantage -$1,000,000-$800,000-$600,000-$400,000-$200,000 $- $200,000 $400,000 $600,000 $800,000 $1,000,0000 50 100 150 200 250 300 350Difference in Objective Value (CAD)DayDeterministic Without Reserves Deterministic with ReservesStochastic with MC Based Scenarios Stochastic with ARMA Based Scenarios118  of the stochastic models appears to be that they miss their daily storage targets by using more water than targeted, creating negative storage deviations that are penalized in the objective function by the marginal value of water.   Figure 52 – Relative cost comparison based on the actual optimal values  Figure 53 - Comparison of the objective function components of models (c) and (d) with model (b).  Some interesting observations are made by studying the results in detail. As can be seen in Figure 51, most of the gains of one model over the others come from few days with large gains and many days 0%10%20%30%40%50%60%70%80%90%100%Deterministicwith ReservesStochastic withMC BasedScenariosStochastic withARMA BasedScenariosRelative Cost-20,000,000-15,000,000-10,000,000-5,000,00005,000,00010,000,00015,000,00020,000,00025,000,00030,000,000Difference in Objective Function ComponentsStochastic with MC Based Scenarios Stochastic with ARMA Based Scenarios119  with small gains. Model (a) results in large gains on days where models (b), (c) or (d) were forced to shed load. Model (a) resulted in minor gains over model (b) on days where the DA capacity reserves that have to be kept for model (b) limited the amount of energy and capacity that can be traded in the energy and AS markets. Model (a) showed minor gains over the stochastic models on days where the wind forecast was accurate, but the uncertainty captured in the WFE scenarios provided sufficient signals to the SO models to increase DA energy imports or to reduce DA exports. Model (b) had large gains over the stochastic models (c) and (d) on days where the stochastic models shed load because of an extreme WFE that was not captured in the WFE scenarios, but the DA capacity reserves kept in model (b) were sufficient to manage this extreme WFE and avoid load shedding. The stochastic models showed large gains over model (b) on days where the extreme WFE was higher than the DA capacity reserves but at least one of the scenarios used in the stochastic models captured this extreme WFE. In these types of days, the scenarios with the extreme WFE usually caused the stochastic models to import more energy in the DA market and avoid load shedding. Smaller gains were obtained from the stochastic models on days where the participation of model (b) in the AS or energy markets was limited because of the DA capacity reserves.  Model (c) showed large gains over model (d) when an extreme WFE was captured by the scenarios generated from the MC-based model but not by scenarios generated from the ARMA-based model. An example of this is day 110. Figure 54 and Figure 55 show the forecasted generation, the actual generation, and the 40 wind generation scenarios for day 110. On this day, the MC-based model’s scenarios were able to capture the extreme WFE (i.e., the WFE after the 19th hour of that day) and the ARMA-based model’s scenarios were not. Similarly, model (d) showed large improvements over model (c) on days where the ARMA model’s scenarios were able to capture the extreme WFE, but the MC-based model’s scenarios were not.  In general, as seen in Figure 52, model (c) performs better than model (d) in the out-of-sample assessment. This is because the scenarios from the MC-based model was, at least for the 2005 wind year, better at capturing extreme positive WFEs (as a reminder, positive WFE means less generation than forecasted). These extreme positive WFE scenarios signal the SO model to make conservative decisions and import high amounts of energy in the DA markets. These conservative decisions usually result in less load shedding events (see Figure 56 below for more details). The MC model performs 120  better at capturing extreme WFEs because it overestimates the ACF as shown in Figure 41. The overestimation of ACF by the MC-based model makes its simulation outputs more persistent than the actual WFE. The occurrence of this high persistence means that high WFE that is encountered at the start of the day is less likely to be significantly lowered before the end of the day. In other words, the MC-based model tends to generate scenarios that have lower/higher WFE that persist over extended periods. This increases the probability that scenarios with extreme WFEs that last the whole day are generated. This is not the correct behavior of the WFE but given the high per unit load shedding cost, the conservative trading decisions that arise from this behavior tend to result in better cost saving in the long run.   Figure 54 - Forecasted, observed and the 40 reduced wind generation scenarios for the 110th day (the MC-based model)  Figure 55 - Forecasted, observed and the 40 reduced wind generation scenarios for the 110th day (ARMA-based model)  Comparison between the expected optimal values and the actual optimal values was also conducted. The expected optimal values were much higher than the actual optimal values on days where the actual wind generation was much lower than forecasted. On these days, the model in RT operation was forced to use more water use than targeted, increase RT imports, or shed the load. The actual optimal values were higher than the expected ones on days where the actual wind realization was much higher than the forecast. This is because the extra generation is exported, and the large 050010001500200025001 4 7 10 13 16 19 22MWHourObserved Forecast050010001500200025001 4 7 10 13 16 19 22MWHourObserved Forecasted121  hydropower plants generated less energy than planned. By generating less energy than planned, the plants stored more water for future use and created a positive storage deviation, which is rewarded in the objective function by the marginal value of water. 5.6.3 Load Shedding The ability of the models to meet the system load in RT operations was also investigated. Figure 45 shows the total cumulative load shed by the different models in the normal study year. It can be seen that model (b), the deterministic model with DA capacity reserves, sheds the most load. This is because the DA capacity reserves are sometimes not sufficient to cover the largest WFEs. When comparing the outputs of models (c) and (d), it is apparent that model (c), on average, is better at avoiding load shedding. This is because the MC-based model’s scenarios are better at capturing extreme WFEs. Extreme WFE scenarios force the SO model to be conservative and to import energy in DA trading. It is also noted that there is no load shedding in model (a), indicating the wind forecast uncertainty is the cause of load shedding in the other models.   Figure 56 - Cumulative load shedding 0500100015002000250030003500400045005000Oct Nov Dec Jan Feb Mar Apr May Jun Jul Aug SepCumulative Load Shed (MW)Deterministic Without Reserves Deterministic with ReservesStochastic with MC Based Scenarios Stochastic with ARMA Based Scenarios122  5.6.4 Wind Power Curtailment Although the cost of wind power curtailment is not very significant, the performance of the four models in avoiding wind curtailment is still an important comparison metric because it shows the ability of the models to accommodate variable generation. Figure 57 shows the cumulative wind power curtailed for the normal water year. It shows that model (b) curtails the highest amount of wind. This is expected since curtailing wind also means reducing the capacity reserves needed to manage wind as discussed in Section 3.7.5. Under some circumstances, such as when the wind forecast is low, using the capacity gained from the reduced reserve is a more profitable decision than using the wind generation. The graph also shows that some amount of wind power is curtailed even in model (a). This curtailment occurred in the spring months when the load was low, ROR generation was high, dispatchable plants were operating near minimum allowable generation, and when electricity market prices were very low or even negative. The SO models performed comparably at minimizing wind power curtailment. Their overall performance at avoiding wind power curtailment is also similar to that of model (a).  Figure 57 - Cumulative wind power curtailment 050010001500200025003000Oct Nov Dec Jan Feb Mar Apr May Jun Jul Aug SepCumulative Wind Curtiled (MW)Deterministic Without Reserves Deterministic with ReservesStochastic with MC Based Scenarios Stochastic with ARMA Based Scenarios123  5.7 Sensitivity Analysis The results discussed so far are specific to the data and assumptions used in the study. Changes in some of the data and in assumptions are expected to affect the estimated DA wind integration cost and stochastic savings. Some of the input data and assumptions that are expected to affect the performance of the models include:  • Inflow conditions in the large hydropower reservoirs  • Unit cost of load shedding • AS market depth • Annual variability of wind • Energy and AS prices  • RT market liquidity • Difference between DA and RT prices  • Transmission limits • DA reserve calculation methodology, etc.  Sensitivity analysis is carried out for the first three inputs listed above and the results are discussed in the following sections. Sensitivity analysis for other inputs and assumptions can be carried out in the future. The sensitivity analysis discussed below is conducted using the out-of-sample results (i.e., actual optimal values).  5.7.1 Sensitivity to Hydrologic Conditions The results of the analysis discussed so far is for a representative normal water year, which was coupled with the 2006 wind year. The sensitivity of the stochastic savings to the water conditions was investigated by using the four models in the dry and wet years. The results from these model runs resulted in a DA wind integration cost to vary from $5.86/MWh for the wet year to $7.95/MWh for the dry year. The slight cost reduction seen in the wet years is the result of excess system capacity in the wet year, which minimized the impact of holding DA capacity reserves and lowered the probability of load shedding. The cost increase in the dry year is partially due to the lower levels of excess capacity in the system in the dry year. Another reason for the cost increase in the dry year is 124  that the wind year which was coupled with the dry year (i.e., wind year 2005, see Table 5) has more days with extreme WFEs that caused load shedding when compared to the wind year used for the normal and wet years (i.e., wind year 2006).  Figure 58 shows the relative cost for different inflow conditions. The stochastic savings are more pronounced in the normal and wet years. The stochastic models performed poorly in the dry year for the same two reasons which were responsible for increasing the wind integration cost in the dry year. The first reason was that the shortage of capacity in the dry year increased the likelihood of load shedding when the WFE was not captured by scenarios used in the SO optimization model. Almost all the increased cost from the model (c) shown in Figure 58 was due to load shedding. When a sensitivity analysis was conducted by reducing the load shedding penalty cost, the stochastic savings in the dry year for both models (c) and (d) increased to 77% (see Section 5.7.2 below for more details). The second reason for the poor performance of the stochastic models in the dry year was that the extreme WFEs of the 2005 wind year were not consistently captured in the set of reduced WFE scenarios. In the dry year, the scenarios generated by the ARMA-based model performed better than the MC-based model in capturing the extreme WFEs. Thus, model (d) resulted in lower costs than model (c). The stochastic savings for the wet year were found to be a bit higher than the savings in the normal year. This was mainly due to the excess capacity in the hydropower system in the wet year, which reduced the likelihood of load shedding events and their corresponding costs.   Figure 58 – Cost comparison for different water years 0%20%40%60%80%100%120%140%Normal Dry WetRelative CostStochastic with MC Based Scenarios Stochastic with ARMA Based Scenarios125  The sensitivity of the stochastic savings and the wind integration cost to the water and wind years suggests that more combinations of water and wind years should be considered in future studies to properly assess the impacts of DA wind integration and the value of using SO for DA scheduling.  5.7.2 Sensitivity to Load Shedding Penalty Cost The per unit load shedding penalty cost assumed for the analysis shown so far was $10000/MWh. This cost is quite high, and because of it, the objective function was heavily penalized on days when the model had no choice but to shed load. To investigate the effect of lower penalty cost, another case was prepared where the penalty cost was assumed to be the minimum of double the hourly energy import price, or $100/MWh. This minimum cost was imposed to discourage load shedding when the market prices were low or negative. This load shed pricing scheme reduced the DA wind integration cost to $1.51/MWh. The relative cost for this penalty cost assumption is shown in Figure 59. It can be seen that the relative cost of models (c) and (d) were reduced to 23%, which implies that the stochastic savings are 77%. The savings increased when the penalty cost was reduced because the stochastic models were able to trade in the energy and AS markets without the severe penalties associated with load shedding. The effect of low shed cost is even more significant in the dry years, where the relative cost decreased dramatically from 130% to 23% for model (c) and from 97% to 23% for model (d). This shows that value of stochastic optimization is even more significant if the penalty for not serving load is low.   Figure 59 – Cost comparison for the low load shedding penalty cost assumption 0%10%20%30%40%50%60%70%80%90%100%Deterministic withReservesStochastic withMC BasedScenariosStochastic withARMA BasedScenariosRelative Cost126  5.7.3 Sensitivity to AS Market Depth Assumption As identified in Section 5.6.2, one of the potential sources of the stochastic savings is increased AS market transactions. In the model used in that section, the AS market transactions were constrained by the available surplus generation capacity and transmission capacity. This is expected to overestimate the amount of AS that can be sold in the market since the AS offered to the market was not limited by market depth. In reality, the depth of the AS market might be much smaller than the AS available in the BC Hydro system. Therefore, a sensitivity analysis was conducted by limiting the maximum possible hourly DA AS transactions by 500 MW and the RT AS transaction by 100 MW. Other model inputs, including the load shedding cost, were left unchanged from the ones used in Section 5.6.2.  The result of this analysis shows that the DA wind integration cost was significantly reduced to $1.06/MWh when the market depth was limited. One of the reasons for this reduction in cost was that transmission capacity is not exhausted by DA AS transactions, which allowed more RT imports and lowered the probability of load shedding. Another reason was that for most periods (i.e., hours) in the study, the binding constraint that limited the amount of AS sold to the market was the market depth constraint instead of the DA reserve constraint or transmission constraint. This market depth constraint equally affected model (a) and model (b) for most periods of the year, so the AS transaction revenues of the two were not very different. This led to the reduced cost of model (b).  Comparing the results of the stochastic models shows a saving of 71% for model (c) and a saving of -69% for model (d). Saving of -69% indicates that model (d) actually increased the cost of wind integration as compared to model (b). It should be noted that most of the increased cost of model (d) is attributed to the cost of load shedding in a single day that was caused by an extreme WFE that was not captured by any of the reduced ARMA-based scenarios. In this case, a single day’s load shedding cost significantly skewed the results of this sensitivity analysis mainly due to the high cost of load shedding penalty which magnified the impacts of the lower market depth assumption. When the load shedding cost was reduced as in Section 5.7.2, the relative cost of model (d) was significantly reduced from 169% to 30%.  127   Figure 60 – Cost comparison for low AS market depth assumption The results of the sensitivity analysis are summarized in Table 13. It is noted from the table that the results of all sensitivity analysis conducted with low load shedding penalty cost resulted in stochastic savings of more than 70%. This shows the negative effect of a high load shedding penalty cost on the value of SO. Table 13 - Summary of sensitivity analysis results Water/Wind Year Modeling Assumptions Wind Integration Cost ($/MWh) Relative Cost Stochastic Saving SO model with MC-based scenarios SO model with MC-based scenarios SO model with MC-based scenarios SO model with MC-based scenarios 1968/2006 - Normal Base Case (PnlShed = $10000/MWh, no market depth assumption) 6.43 30% 52% 70% 48% Low Per Unit Load Shedding Cost (PnlShed = max(S100/MWh,2*ExPrt)) 1.51 23% 23% 77% 77% Low Market Depth (DA AS<=500, RT AS<=100) 1.06 29% 169% 71% -69% Low Per Unit Load Shedding Cost, Low Market Depth 1.06 29% 28% 71% 72% 1969/2005 - Dry Base Case 7.95 130% 96% -30% 4% Low Per Unit Load Shedding Cost 1.91 23% 23% 77% 77% 1966/2006 - Wet Base Case 5.86 38% 56% 62% 44%  0%20%40%60%80%100%120%140%160%180%Deterministic withReservesStochastic with MCBased ScenariosStochastic with ARMABased ScenariosRelative Cost128   Conclusion and Recommendations for Future Work This thesis had four goals: to statically characterize the WFE, to propose two scenario generation algorithms that create representative WFE scenarios, to quantify the value of SO in reducing the cost of wind uncertainty in the BC Hydro system, and to assess the effect of the scenario generation technique on the stochastic savings. To achieve these goals, first, synthetic wind forecast and wind generation data for various wind farms in BC were collected and statistically analyzed. Then, two scenario generation algorithms were developed to create WFE scenarios that accurately represent the relevant WFE statistical properties. Many scenarios were generated using these two algorithms. The numbers of the generated scenarios were then reduced by using scenario reduction algorithm to lower the computational burden of running the SO models. This reduced sets of scenarios were subsequently used in SO models that were developed by extending and enhancing an existing BC Hydro planning and operations model. Finally, by analyzing the results of the SO runs, the value of SO in reducing the cost of DA wind forecast uncertainty in the BC Hydro system was determined. This chapter summarizes the results of this process under four subsections that correspond to the four goals of the thesis. A future work section is also provided to provide recommendations for future research on this topic.  6.1 Statistical Properties of the WFE The behavior of the DA WFE in the BC Hydro system has been investigated in this thesis. The WFE was found to be significant in the DA time frame, with an MAE of about 10% of the installed capacity.  The magnitude of the WFE was found not to be significantly influenced by the forecast lead time or the season of the year. The DA wind forecast was found to be negatively biased, i.e., the wind power generation is underestimated by 2.39% of the installed capacity on average. The WFE has a high kurtosis value of 4.88, which made it difficult to find a parametric distribution that can accurately capture its properties. One of the important observations made about the property of the WFE was that it is significantly influenced by the wind forecast level. Low and high forecast levels tend to have low WFE and low variance while intermediate forecast levels were found to have high WFE and high variance. The WFE is also found to be a highly persistent second-order AR process.  129  6.2 Scenario Generation Two scenario generation algorithms were developed in this research to capture the statistical properties of the WFE. The first is a second order MC model with an exogenous variable. The second is an ARMA model which uses conditional probability distributions for probability distribution transformation. Both MC and ARMA based models were found to be effective in creating WFE scenarios that can capture the statistical properties of the WFE stochastic process, namely, the joint probability distribution, the ACF, and PACF. Both models were fast and efficient, and had low computational burdens. Both models required the generation of large number of scenario samples to confidently capture the underlying WFE probability distribution, which necessitated the use of a scenario reduction algorithm before the scenarios can be used in a SO problem. Although both models were found to be satisfactory in terms of capturing the statistical parameters of the WFE, the ARMA-based model was found to perform better at capturing the statistical properties, especially the PACF and the joint probability distribution of the WFE and wind forecast. Comparison of the two models’ outputs also highlighted the advantages and disadvantages of using the two scenario generation algorithms. The MC-based model has the advantage of being simpler to formulate as it does not require probability transformations or fitting a parametric probability distribution. Modeling a multivariate stochastic process using the MC-based model is also simpler. Contrary to what other researchers have discovered, the MC model used in this study overestimates the ACF. This is primarily due to the MC-based model’s inability to reproduce the lag-2 negative PACC of the WFE. This leads to the conclusion that the MC-based model used in this thesis is not appropriate to model a stochastic process that has a significant anti-persistence property.  Most of the underperformance of the MC-based model was caused by the discretization of the continuous WFE. Discretization masks some of the properties of the stochastic process. For example, the negative lag-2 PACC of the WFE is masked after discretization. Another drawback of the MC-based model was its inability to predict unobserved state transitions as the transition probabilities were estimated empirically based only on observed events. The MC-based model also had numerous model parameters, which necessitates a large input data set for accurate parameter estimation.  130  The ARMA-based model avoids the pitfalls of the MC model as it does not require discretization and its prediction are not limited by observed events. The ARMA-based model has much fewer parameters, which means the model needs less data for parameter estimation and is more parsimonious. It is found that the ARMA model is able to explicitly model the autocorrelation of the WFE, which makes it more successful at modeling the anti-persistence property of the WFE. These factors combined help the ARMA-based model to be better at representing the ACF, PACF, CDF, and statistical moments of the WFE. However, the ARMA-based model has some drawbacks. Working with ARMA model can be more complicated since it requires specialized statistical software for model fitting and simulation. In addition, the ARMA-based model required probability transformation, which causes a loss of autocorrelation information and adds another source of modeling difficulty. Modeling a multivariate stochastic process using ARMA models is also more challenging as compared to MC-based models. 6.3 The Value of Stochastic Optimization in Reducing the DA Wind Integration Cost Stochastic optimization has considerable value in reducing the cost of DA wind forecast uncertainty in the BC Hydro system. It was shown in this thesis that SO may be able to reduce the DA wind integration cost by up to 70% under normal operating conditions. Most of the cost savings come from reduced DA imports, increased AS market transactions and reduced load shedding events. The SO models were also found to recommend operation plans which on average were better at avoiding load shedding and wind power curtailment events.  The DA wind integration cost and the stochastic savings were found to be sensitive to the input data and assumptions. The savings were higher in a normal year and a wet year than in a dry year. When the cost of load shedding was reduced, the integration cost decreased substantially and the stochastic savings increased significantly. Furthermore, limiting the AS market depth caused a decrease in both the DA wind integration cost and the stochastic savings. The results from this limited sensitivity analysis indicate that care should be taken when developing assessment strategies for wind integration cost and stochastic savings. The load shedding penalty cost requires a special attention as it was found to significantly affect the stochastic savings. The results from the limited cases investigated in this thesis show that stochastic optimization has a great potential to lower the DA wind integration cost. 131  However, investigation of more water/wind year combinations with more comprehensive sensitivity analysis is required before making final conclusions. The importance of conducting an out-of-sample assessment has also been demonstrated. As compared to the results obtained from the assessment that used the expected optimal values, the out-of-sample assessment results were found to have higher DA wind integration cost and lower stochastic savings. This suggests that an out-of-sample assessment should be carried out when possible to more accurately examine the benefit of a modeling approach.   6.4 The Effect of the Scenarios Generation Method in the Stochastic Savings In the normal water year and wet water year, the SO model with scenarios generated using the MC-based algorithm performed better than the SO model with scenarios generated using the ARMA-based model. This was rather unexpected since the ARMA-based scenarios appeared to be better at capturing the statistical properties of the WFE. The MC-based scenarios led to higher stochastic savings in these water years mainly due to their ability to capture extreme WFEs.  This ability is the unintended consequence of the MC-based model’s overestimation of ACF.  However, the SO model with scenarios from the ARMA-based model performed slightly better than that of the MC-based model in the dry year. The savings from both MC and ARMA based models were also comparable when the load shedding cost was reduced. This indicates that further research with many years of data is needed to confidently identify the advantages of one model over the other in terms of increasing the stochastic savings. From the analysis conducted, it was clear that the property that is most critical in stochastic optimization is the ability to capture the extreme WFEs. This finding also reinforces the need for a better scenario reduction algorithm that preserves the extreme WFEs and better represents the tails of the WFE distribution.  6.5 Future Work The following is a list of suggested extensions to the research carried out in this thesis. 1. A better scenario reduction algorithm that does not use the Euclidean norm as a similarity metric can be investigated. A scenario reduction algorithm that uses the approaches suggested in Bruninx (2016) and Morales et al. (2009) can be investigated to find out if significant 132  improvement can be achieved. The sensitivity of the stochastic saving to the number of scenarios retained after the scenario reduction can also be investigated.  2. A study that treats the load forecast and market price forecast as additional uncertain variables can be conducted. This is especially relevant given that more uncertainty is seen in the market price forecast than in the wind generation forecast. In such study, a scenario generation algorithm that takes into account the correlation between wind generation forecast, electricity market prices, and AS prices should be developed. The effects the WFE and LFE have on RT market prices should also be examined and incorporated in such study. 3. The advantages or disadvantages of using a WFE scenario generation algorithm that is based on wind speed forecasts for individual wind farms instead of a single aggregate wind generation forecast for all the wind farms can be investigated. This scenario generation algorithm should consider the spatial and temporal autocorrelation of the wind speed forecast errors at the various wind farms.  4. In this study, only 15% wind penetration level was considered, but the cost of wind integration and expected benefits of using SO models are known to depend on wind penetration level. Therefore, the sensitivity of the stochastic savings to different wind penetration levels can be assessed. 5. Given the variability of the stochastic savings to water year and wind year, longer study periods with various combinations of water and wind years can be considered to better estimate the value of SO models in reducing the DA wind integration cost. The correlation between the weather variables precipitation and the wind should also be incorporated in such multi-year studies. The comparison between the MC-based model and the ARMA-based model can also be extended for various water and wind year combinations. Such investigations may lead to different conclusions.  6. If long record of WFE data is available, the seasonality of the WFE can be thoroughly investigated and incorporated into a scenario generation algorithm.  133  Bibliography Angarita, J.L., Usaola, J. & Martínez-Crespo, J., 2009. Combined hydro-wind generation bids in a pool-based electricity market. Electric Power Systems Research, 79(7), pp.1038–1046. Archer, C.L. et al., 2017. The challenge of integrating offshore wind power in the U.S. electric grid. Part I: Wind forecast error. Renewable Energy, 103, pp.346–360.  AWEA Data Services, 2017. U.S. Wind Industry First Quarter 2017 Market Report, Available at: http://awea.files.cms-plus.com/FileDownloads/pdfs/1Q2017 AWEA Market Report Public Version.pdf. [Accessed May 18, 2017]. Bathurst, G.N., Weatherill, J. & Strbac, G., 2002. Trading wind generation in short term energy markets. IEEE Transactions on Power Systems, 17(3), pp.782–789. BC Hydro, 2017. 2016/17 Annual Service Plan Report, Vancouver, Canada. Available at: https://www.bchydro.com/content/dam/BCHydro/customer-portal/documents/corporate/accountability-reports/financial-reports/annual-reports/bchydro-2016-17-annual-service-plan-report.pdf. [Accessed June 8, 2018]. BC Hydro, 2010. Clean Power Call Request for Proposals: Report on the RFP Process, Vancouver, Canada. Available at: https://www.bchydro.com/content/dam/hydro/medialib/internet/documents/planning_regulatory/acquiring_power/2010q3/cpc_rfp_process_report.pdf. [Accessed May 7, 2017]. BC Hydro, 2013. Integrated Resource Plan, Vancouver, Canada. Available at: https://www.bchydro.com/energy-in-bc/planning-for-our-future/irp/current-plan/document-centre/reports/november-2013-irp.html. [Accessed May 7, 2017]. Bird, L. et al., 2005. Policies and market factors driving wind power development in the United States. Energy Policy, 33(11), pp.1397–1407. Bludszuweit, H., Dominguez-Navarro, J.A. & Llombart, A., 2008. Statistical Analysis of Wind 134  Power Forecast Error. IEEE Transactions on Power Systems, 23(3), pp.983–991.  Boone, A., 2005. Simulation of short-term wind speed forecast errors using a multi-variate ARMA (1, 1) time-series model. Royal Institute of Technology.  Box, G.E.P., Jenkins, G.M. & Reinsel, G.C., 2008. Time Series Analysis 4th ed., Hoboken, NJ: John Wiley & Sons, Inc. Brokish, K. & Kirtley, J., 2009. Pitfalls of Modeling Wind Power Using Markov Chains. Power Systems Conference and Exposition, pp.1–6. Bruninx, K., 2016. Improved modeling of unit commitment decisions under uncertainty. KU Leuven. Canadian Wind Energy Association, 2016. Wind Energy in Canada, Available at: http://canwea.ca/wp-content/uploads/2017/04/canadas-wind-energy-story.pdf. [Accessed May 18, 2017]. Charles River Associates, 2009. SPP WITF Wind Integration Study, Little Rock, AR. Available at: http://www.uwig.org/CRA_SPP_WITF_Wind_Integration_Study_Final_Report.pdf. [Accessed August 25, 2017]. Chen, P. et al., 2010. ARIMA-based time series model of stochastic wind power generation. IEEE Transactions on Power Systems, 25(2), pp.667–676. Chong, E.K. & Zak, S.H., 2013. An Introduction to Optimization Vol 76., John Wiley & Sons. Conejo, A.J., Carrión, M. & Morales, J.M., 2010. Decision Making Under Uncertainty in Electricity Markets, Boston, MA: Springer US.  Dai, T. & Qiao, W., 2015. Optimal Bidding Strategy of a Strategic Wind Power Producer in the Short-Term Market. IEEE Transactions on Sustainable Energy, 6(3), pp.707–719. DNV-GEC, 2009. BC Hydro Wind Data Study, Available at: https://www.bchydro.com/content/dam/hydro/medialib/internet/documents/environment/w135  inddata/pdf/wind_data_study_report_may1_2009.pdf. [Accessed May 5, 2017]. Dowds, J. et al., 2015. A review of large-scale wind integration studies. Renewable and Sustainable Energy Reviews, 49, pp.768–794. Dupačová, J., Gröwe-Kuska, N. & Römisch, W., 2003. Scenario reduction in stochastic programming: An approach using probability metrics. Mathematical Programming Series A., 511, pp.493–511. Dvorkin, Y., Member, S. & Wang, Y., 2014. Comparison of Scenario Reduction Techniques for the Stochastic Unit Commitment. In PES General Meeting| Conference & Exposition, 2014 IEEE. Seattle, Washington: IEEE, p. 5. EnerNex Corporation, 2011. Eastern Wind Integration and Transmission Study, Golden, Colorado. Available at: http://www.nrel.gov/wind/systemsintegration/pdfs/2010/ewits_final_report.pdf. [Accessed August 25, 2017]. Evans, J.I., 2009. Benefits of wind power curtailment in a hydro-dominated electric generation system. University of British Columbia. Follestad, T., Wolfgang, O. & Belsnes, M.M., 2011. An Approach for Assessing the Effect of Scenario Tree Approximations in Stochastic Hydropower Scheduling. In 17th Power Systems Computation Conference. pp. 271–277. Garcia-Gonzalez, J. et al., 2008. Stochastic Joint Optimization of Wind Generation and Pumped-Storage Units in an Electricity Market. IEEE Transactions on Power Systems, 23(2), pp.460–468. GE Energy, 2016. Pan-Canadian Wind Integration Study (PCWIS), Ottawa, Ontario. Available at: http://canwea.ca/wind-energy/wind-integration-study/. [Accessed November 15, 2017]. GE Energy, 2010. Western Wind and Solar Integration Study:, Golden, Colorado. Available at: http://www.osti.gov/bridge. [Accessed August 25, 2017]. 136  Giebel, G. & Sørensen, Poul Ejnar Holttinen, H., 2007. TradeWind Deliverable 2.2: Forecast error of aggregated wind power. Risø-I, 2567. Goodarzi, E., Ziaei, M. & Hosseinipour, E.Z., 2014. Introduction to Optimization Analysis in Hydrosystem Engineering, Springer International Publishing. Greenhall, A., 2013. Wind Scenarios for Stochastic Energy Scheduling. University of Washington. Gröwe-Kuska, N., Heitsch, H. & Römisch, W., 2003. Scenario reduction and scenario tree construction for power management problems. 2003 IEEE Bologna PowerTech - Conference Proceedings, 3, pp.152–158. Guzman, H.A.R., 2010. Value of Pumped-Storage Hydro for Wind Power Integration in the British Columbia Hydroelectric System. The University of British Columbia. GWEC, 2017. Global Wind Report 2016, Brussels, Belgium. Available at: http://files.gwec.net/files/GWR2016.pdf. [Accessed October 19, 2017]. Heitsch, H. & Römisch, W., 2003. Scenario Reduction Algorithms in Stochastic Programming. Computational Optimization and Applications, 24(2/3), pp.187–206.  Hibiki, N., 2006. Multi-period stochastic optimization models for dynamic asset allocation. Journal of Banking and Finance, 30(2), pp.365–390. Hocaoglu, F.O., Gerek, Ö.N. & Kurban, M., 2008. The Effect of Markov Chain State Size for Synthetic Wind Speed Generation. Probabilistic Methods Applied to Power Systems, 2008. PMAPS’08. Proceedings of the 10th International Conference, pp.2–5. Hodge, B.-M. et al., 2012. A Comparison of Wind Power and Load Forecasting Error Distributions. In 2012 World Renewable Energy Forum. p. 8.  Hodge, B. et al., 2012. Wind Power Forecasting Error Distributions. In The 11th Annual International Workshop on Large-Scale Integration of Wind Power into Power Systems as well as on Transmission Networks for Offshore Wind Power Plants Conference. Lisbon, Portugal 137  Hodge, B.M. & Milligan, M., 2011. Wind power forecasting error distributions over multiple timescales. IEEE Power and Energy Society General Meeting, pp.1–8. Holttinen, H. et al., 2009. Design and operation of power systems with large amounts of wind power, Vuorimiehentie, Finland: Julkaisija, Utgivare. Høyland, K. & Wallace, S.W., 2001. Generating Scenario Trees for Multistage Decision Problems. Management Science, 47(2), pp.295–307.  Hyndman, R.J., 2017. forecast: Forecasting functions for time series and linear models. Available at: http://github.com/robjhyndman/forecast. [Accessed Feb 14, 2018]. Idaho Power, 2013. Wind Integration Study Report, Available at: https://www.idahopower.com/AboutUs/PlanningForFuture/WindStudy/default.cfm. [Accessed May 2, 2017]. IEA, 2015. World Energy Outlook 2015, Paris, France. Available at: https://www.iea.org/publications/freepublications/publication/WEO2015.pdf.  [Accessed November 15, 2017]. Jung, J. & Broadwater, R.P., 2014. Current status and future advances for wind speed and power forecasting. Renewable and Sustainable Energy Reviews, 31, pp.762–777. Keutchayan, J., Gendreau, M. & Saucier, A., 2016. Quality Evaluation of Scenario-Tree Generation Methods for Solving High- Dimensional Stochastic Programs Quality Evaluation of Scenario-Tree Generation Methods for Solving High- Dimensional Stochastic Programs. CIRRELT, (September). Labadie, J.W., 2004. Optimal Operation of Multireservoir Systems: State-of-the-Art Review. Journal of Water Resources Planning and Management, 130(2), pp.93–111.  Latorre, J.M., Cerisola, S. & Ramos, A., 2007. Clustering algorithms for scenario tree generation: Application to natural hydro inflows. European Journal of Operational Research, 181(3), 138  pp.1339–1353. Li, J., Lan, F. & Wei, H., 2016. A Scenario Optimal Reduction Method for Wind Power Time Series. IEEE Transactions on Power Systems, 31(2), pp.1657–1658. Löhndorf, N., 2016. An empirical analysis of scenario generation methods for stochastic optimization. European Journal of Operational Research, 255(1), pp.121–132. Lowery, C. & O’Malley, M., 2012. Impact of Wind Forecast Error Statistics Upon Unit Commitment. IEEE Transactions on Sustainable Energy, 3(4), pp.760–768.  Lowery, C. & O’Malley, M., 2013. Wind Power Scenario Tree Tool: Development and Methodology. In Reliability and Risk Evaluation of Wind Integrated Power Systems. India: Springer India, pp. 13–27.  Ma, X. et al., 2013. Scenario Generation of Wind Power Based on Statistical Uncertainty and Variability. IEEE Transactions on Sustainable Energy, 4(4), pp.894–904. Matevosyan, J. & Söder, L., 2006. Minimization of imbalance cost trading wind power on the short-term power market. IEEE Transactions on Power Systems, 21(3), pp.1396–1404. MEMPR, 2007. The BC Energy Plan:A Vision for Clean Energy Leadership. Ministry of Energy, Mines and Petroleum Resources, p.25. Available at: http://www.energybc.ca/cache/usage/usage3/www.energyplan.gov.bc.ca/PDF/BC_Energy_Plan.pdf. [Accessed May 18, 2017]. Miettinen, J.J. & Holttinen, H., 2017. Characteristics of day-ahead wind power forecast errors in Nordic countries and benefits of aggregation. Wind Energy, 20(6), pp.959–972.  Morales, J.M. et al., 2009. Scenario Reduction for Futures Market Trading in Electricity Markets. IEEE Transactions on Power Systems, 24(2), pp.878–888. Morales, J.M., Conejo, A.J. & Perez-Ruiz, J., 2010. Short-Term Trading for a Wind Power Producer. IEEE Transactions on Power Systems, 25(1), pp.554–564. 139  Morales, J.M., Mínguez, R. & Conejo, A.J., 2010. A methodology to generate statistically dependent wind speed scenarios. Applied Energy, 87(3), pp.843–855.  National Energy Board, 2016. Canada’s Renewable Power Landscape - Energy Market Analysis 2016. , pp.1–96. Available at: https://www.neb-one.gc.ca/nrg/sttstc/lctrct/rprt/2016cndrnwblpwr/2016cndrnwblpwr-eng.pdf. [Accessed November 16, 2017] Nfaoui, H., Essiarab, H. & Sayigh, A.A.M., 2004. A stochastic Markov chain model for simulating wind speed time series at Tangiers, Morocco. Renewable Energy, 29(8), pp.1407–1418. Pacific Corp, 2015. Pacific Corp 2015 IRP: Appendix H –Wind Integration Study, Portland, Oregon. PandŽić, H. et al., 2013. Offering model for a virtual power plant based on stochastic programming. Applied Energy, 105, pp.282–292. Papaefthymiou, G. & Klockl, B., 2008. MCMC for Wind Power Simulation. IEEE transactions on energy conversion, 23(1), pp.234–240. Papavasiliou, A. & Oren, S.S., 2013. Multiarea Stochastic Unit Commitment for High Wind Penetration in a Transmission Constrained Network. Operations Research, 61(3), pp.578–592. Papavasiliou, A., Oren, S.S. & Neill, R.P.O., 2011. Reserve Requirements for Wind Power Integration : A Stochastic Programming Framework. IEEE Transactions on Power Systems, 26(4), pp.2197–2206. Pinson, P. et al., 2009. From probabilistic forecasts to statistical scenarios of short-term wind power production. Wind Energy, 12(1), pp.51–62. R Core Team, 2017. R: A language and environment for statistical computing. Rahimiyan, M., Morales, J.M. & Conejo, A.J., 2011. Evaluating alternative offering strategies for wind producers in a pool. Applied Energy, 88(12), pp.4918–4926. 140  Ruiz, P.A., Philbrick, C.R. & Sauer, P.W., 2009. Wind power day-ahead uncertainty management through stochastic unit commitment policies. In 2009 IEEE/PES Power Systems Conference and Exposition. IEEE, pp. 1–9. Sahin, A.D. & Sen, Z., 2001. First-order Markov chain approach to wind speed modelling. Journal of Wind Engineering and Industrial Aerodynamics, 89(3–4), pp.263–269. Schwabe, P., Lensink, S. & Hand, M., 2011. IEA Wind Task 26, Available at: http://www.nrel.gov/docs/fy11osti/48155.pdf. [Accessed October 19, 2017] Shamshad, A. et al., 2005. First and second order Markov chain models for synthetic generation of wind speed time series. Energy, 30(5), pp.693–708. Sharma, K.C., Jain, P. & Bhakar, R., 2013. Wind Power Scenario Generation and Reduction in Stochastic Programming Framework. Electric Power Components and Systems, 41(3), pp.271–285. Shawwash, Z.K., 2000. A decision support system for real-time hydropower scheduling in a competitive power environment market. University of British Columbia. Sheppard, M., 2012. Fit all valid parametric probability distributions to data. ALLFITDIST Matlab code. Available at: https://www.mathworks.com/matlabcentral/fileexchange/34943-fit-all-valid-parametric-probability-distributions-to-data?focused=5228686&tab=function. [Accessed May 15, 2017] Soder, L., 2004. Simulation of wind speed forecast errors for operation planning of multiarea power systems. 8th International Conference on Probabilistic Methods Applied to Power Systems, pp.723–728. Sumaili, J. et al., 2011. Clustering-based wind power scenario reduction. Proceedings of the 17th Power Systems Computation Conference. Surespanwind, 2017. Canada’s Support of Wind Energy: Grants and Other Incentives. Available 141  at: http://www.surespanwind.com/canadas-support-wind-energy-grants-incentives/ [Accessed October 13, 2017]. Tagliaferri, F. et al., 2016. Wind modelling with nested Markov chains. Journal of Wind Engineering and Industrial Aerodynamics, 157, pp.118–124. Tang, J., Brouste, A. & Tsui, K.L., 2015. Some improvements of wind speed Markov chain modeling. Renewable Energy, 81, pp.52–56. Available at:  Tewari, S., Geyer, C.J. & Mohan, N., 2011. A statistical model for wind power forecast error and its application to the estimation of penalties in liberalized markets. IEEE Transactions on Power Systems, 26(4), pp.2031–2039. Tuohy, A. et al., 2008. Benefits of Stochastic Scheduling for Power Systems with Significant Installed Wind Power. Probabilistic Methods Applied to Power Systems, 2008. PMAPS ’08. Proceedings of the 10th International Conference on, pp.1–7. U.S. Department of Energy, 2013. Federal Incentives for Wind Power. , pp.13–14. Available at: https://www.nrel.gov/docs/fy13osti/57933.pdf. [Accessed November 15, 2017] Vagropoulos, S.I. et al., 2016. ANN-based scenario generation methodology for stochastic variables of electric power systems. Electric Power Systems Research, 134, pp.9–18. Xu, B. et al., 2015. Scenario tree reduction in stochastic programming with recourse for hydropower operations. Water Resources Research, 51(8), pp.6359–6380. Yeh, W.W., 1985. Reservoir Management and Operations Models: A State-of-the-Art Review. Water Resources Research, 21(12), pp.1797–1818. Available at:  Zima-Bočkarjova, M. et al., 2010. Sharing of profit from coordinated operation planning and bidding of hydro and wind power. IEEE Transactions on Power Systems, 25(3), pp.1663–1673.  

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                        
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            src="{[{embed.src}]}"
                            data-item="{[{embed.item}]}"
                            data-collection="{[{embed.collection}]}"
                            data-metadata="{[{embed.showMetadata}]}"
                            data-width="{[{embed.width}]}"
                            async >
                            </script>
                            </div>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.24.1-0369051/manifest

Comment

Related Items