Advanced Scenario Creation Strategies for Stochastic Economic Dispatch with Renewables
Abstract
Realtime dispatch practices for operating the electric grid in an economic and reliable manner are evolving to accommodate higher levels of renewable energy generation. In particular, stochastic optimization is receiving increased attention as a technique for handling the inherent uncertainty in wind and solar generation. The typical twostage stochastic optimization formulation relies on a sample average approximation with scenarios representing errors in forecasting renewable energy ramp events. Standard Monte Carlo sampling approaches can result in prohibitively highdimensional systems for optimization, as well as a poor representation of extreme events that challenge grid reliability. We propose two alternative scenario creation strategies, importance sampling and Bayesian quadrature, that can reduce the estimator’s variance. Their performance is assessed on a week’s worth of 5 minute stochastic economic dispatch decisions for realistic wind and electrical system data. Both strategies yield more economic solutions and improved reliability compared to Monte Carlo sampling, with Bayesian quadrature being less computationally intensive than importance sampling and more economic when considering at least 20 scenarios.
Nomenclature

[]

Dispatch of thermal generators

Cost of thermal generator

Loss function

Dispatch of wind generator

Cost of wind plant

Forecasted wind generation

Error in wind forecast

Max output of generator

Min output of generator

Max ramp up of generator

Max ramp down of generator

Set of thermal generators

Set of timesteps

Set of wind plants

Set of wind forecast error scenarios

Set of load buses

Value of loads

Lossofload

Excess capacity

Cost of Lossofload

Cost of excess capacity

Spilled wind

Cost of spilled wind
I Introduction
Modern practices for economic dispatch of the electrical grid are rooted in deterministic problem formulations. The increasing deployment of stochastic renewable energy generators, such as wind and solar, requires dispatch techniques that are consistent with the reliability and economic operating requirements of the electric grid. Stochastic optimization techniques have received increased attention from the power systems community, and are a promising way to address the randomness in renewable energy generation. This paper explores two novel approaches to scenariobased stochastic optimization that yield lower costs and reduced reserve requirements, facilitating the grid integration and further deployment of renewable energy generators.
As access to high performance computing resources grows, grid operators are able to consider increasingly sophisticated optimization problems in their dispatch decisions. These complexities can include considering AC optimal power flow (ACOPF), larger grids with thousands of buses, multiple stages or time periods, and stochastic optimization techniques. These techniques can improve grid reliability in the face of uncertain generation, and help remove obstacles to further renewable energy deployment.
An emerging approach to dispatching generators in the face of uncertain renewable generation is to formulate and solve a twostage stochastic programming problem with recourse. In twostage recourse problems a first stage decision is made, followed by a recourse stage where adjustments can be made after a random variable is realized. In the context of economic dispatch for the grid, this involves making initial generator dispatch decisions that minimize the combined costs of the first stage and the expectation of the recourse costs. The expected recourse costs are calculated with respect to the distribution of the uncertain renewable generation, and are typically estimated with a sample average approximation using Monte Carlo techniques [1, 2]. However, Monte Carlo techniques exhibit slow convergence with respect to the number of samples and draw samples that are not necessarily a good test of grid robustness in the event of large ramps in generation, creating a need for improved sampling techniques.
In this article we report on two techniques for creating scenarios of errors in forecasting wind generation that improve the convergence rates and accuracy when calculating the expected recourse costs, ultimately achieving more economic dispatch solutions. The created scenarios also contain ramp events that are more challenging for the grid, resulting in dispatch decisions that produce more reliable grid operation. The cost savings and reliability improvements are demonstrated in a weeklong simulation of 5 minute dispatch decisions on realistic wind and grid data.
In Section II we provide a background on stochastic programming techniques for twostage recourse problems, the standard sample average approximation, and power grid applications in economic dispatch. In Section III we introduce new scenario creation strategies using importance sampling and Bayesian quadrature. In Section IV we demonstrate the implementation of these techniques in a stochastic economic dispatch framework and examine their performance. Finally, we conclude in Section V with a discussion of the results and extensions to more complex multistage problems and more complex network topologies.
Ii TwoStage Stochastic Programming
Iia General Problem Formulation
The general twostage stochastic programming problem [3, 4] is defined as:
(1) 
where are the first stage decision variables, is the secondstage cost function defined by,
s. t.  (2) 
and are the secondstage variables, and are random variables. General constraint equations are represented by , and , and will be made specific to grid applications in the next section. We restrict ourselves to using the linear stochastic programming problem, i.e.
(3) 
where the recourse function is defined as the solution to,
s. t.  (4)  
In the context of economic dispatch, the first stage variables represent the dispatch levels of all generators, are linear costs, and the expectation is taken with respect to the random variable , which represents errors in forecasting the available wind generation.
The second stage recourse variable, , represents adjustments to generation using reserves, such as up/down regulation, flex reserve or automatic generation control (AGC), that are necessary based on the amount of wind generation that materializes.
The optimization problem embedded inside the expectation makes it difficult to solve, so a sample average approximation (SAA) is commonly used to represent the expectation. In the SAA approach, we introduce a set of scenarios, , and approximate the expectation as . Commuting the summation and minimization, we can simplify to the extensive form of a twostage stochastic programming problem:
(5) 
A major challenge in the SAA approach is selecting the scenarios such that the expectation is accurately and efficiently calculated, while also ensuring reliable operation of the grid during different levels of wind generation. In Monte Carlo approaches, the samples are drawn from their nominal distribution and given equal weights. This approach yields convergence, requiring prohibitively many samples to achieve an acceptable error level. Other methods for selecting scenarios have been recently proposed, including scenarios motivated by grid reserve requirements [5], chance constraints [6, 7], polynomial chaos expansions [8], or epispline basis functions [9].
IiB TwoStage Stochastic Economic Dispatch
In this is study we restrict ourselves to studying the stochastic economic dispatch problem, i.e. balancing the (stochastic) demand for power across a grid with the (stochastic) power generation from all generators on the network. This study considers aggregated wind and load data, and therefore our model is formulated without network constraints. However, the scenario creation methods presented in this study should generalize to more complicated constraints on the network, e.g. DC and ACOPF. For a review of the economic dispatch problem, as well as the DCOPF and ACOPF problems, we recommend [10, 11, 12]. For an introduction to stochastic grid operations we recommend [13].
In the context of grid operations, the stochastic economic dispatch problem takes the following form:
s. t.  (6)  
where for in the set of generators, , are the first stage decisions representing thermal dispatch, and is the recourse function for a particular realization of the random variable , which represents the error in forecasting wind generation at each wind plant in the set of wind plants . The recourse function is defined as the solution to the following optimization problem:
s. t.  (7)  
where and for each bus in the set of buses . The expectation is taken with respect to the random variable . Generator ramp constraints are imposed between sequential timesteps, and not allowed to participate in recourse. The loss of load or excess capacity represented by is a slack variable that facilitates finding feasible solutions and represents a challenge to grid reliability and incurs substantial costs.
Iii Strategies for Scenario Generation
Iiia Importance Sampling
Solving the twostage stochastic economic dispatch problem requires a method of selecting scenarios. A logical starting point would be to draw from the nominal distribution of forecast errors in the wind power. However, this approach will not take into account costs communicated by the loss function (IIB), potentially resulting in scenarios that do not include information about expensive events. As an example, consider a single wind farm generating power for one time period. Assuming that the nominal distribution of forecast errors is symmetric about zero, simply sampling from this distribution will produce scenarios that represent positive forecast errors, i.e. more wind than forecast, about half the time. If we use a small number of scenarios, e.g. to , as is common in stochastic optimization problems, the scenario set may not include scenarios that represent negative forecast errors. This would be a costly mistake since both lossofload and reserves are more expensive than dispatching extra generation during the first stage decision process.
To address this issue with random sampling we employ importance sampling, see e.g. [16]. The reason for using importance sampling is that by including information from the cost function , we construct a new probability density function that assign more mass to values of corresponding to errors in overpredicting the available wind generation. Additionally, the construction of the importance distribution also provides a formula for properly weighting these samples to ensure the weighted average actually converges to . Finally, importance sampling can also yield a lower variance estimate of the expected recourse cost.
Recall that our expected recourse cost is
(8) 
where is a probability density function on and for all . Introducing a new pdf on (the importance distribution) we can write
(9)  
(10)  
(11) 
where denotes expectation with distributed as . The ratio is an adjusted weight on sample evaluations of that reflects sampling from instead of .
The importance sampling estimate of is then
(12) 
where are realizations of the random variable with distribution . The optimal importance sampling distribution, defined as , is a zerovariance estimator, and computes the correct expectation from a single sample. Such an estimator suffers from the curse of circularity (the quantity being estimated, , is required to properly normalize the importance distribution), but provides guidance on constructing an approximate importance sampling distribution.
Our approach to constructing an importance distribution is to solve a deterministic economic dispatch problem to approximate and , and then use them to build an approximation of . Our first simplification involves approximating the loss function as , i.e. we assume that the wind forecast is perfect. This step reduces the two stage stochastic optimization problem to a simple linear program:
s. t.  (13)  
which we solve for . Given the solution to the simplified dispatch problem we compute the expectation via the trapezoidal rule. We remark that this process is fast since each evaluation of the loss function is computed by solving a linear program. Given the approximated expected cost, , the importance sampling distribution is defined as
(14) 
and is then sampled instead of in the solution of the original stochastic economic dispatch problem.
IiiB Gaussian Process Representation
One shortcoming in the importance sampling approach presented in the preceding subsection is that in reality, is a firststage variable in a twostage stochastic optimization problem, and its value is computed at the same time as the solution to the rest of the problem. By ignoring this coupling, our proxy may be off by a significant amount, leading to an importance distribution with an incorrect shape and location.
We propose an alternative to using the importance sampling approach: representing the loss function as a Gaussian process (GP) and computing its expectation via a Bayesian quadrature technique. One benefit of this approach is that handles uncertainty in the second stage loss function introduced by not knowing the first stage dispatch. In other words, we can treat as a random function when is unknown. Additionally, computing expectations of the GP with Bayesian quadrature does not require knowledge of the loss function to select the quadrature points. Instead, the Bayesian quadrature algorithm spreads out the sample points in the domain based on the nominal distribution and the Gaussian process covariance kernel function. The Bayesian quadrature algorithm also learns quadrature weights, thus ensuring that the weighted sum of the loss function sampled at the quadrature points returns the correct value for the expectation loss. A comprehensive review of Gaussian processes is provided by Rasmussen and Williams [17], while a description of the Bayesian quadrature sampling strategy can be found in [18].
By representing the loss function as a GP, we make the assumption that the joint distribution of any finite number of loss function evaluations is a multivariate Gaussian. The GP is fully specified by a mean function and a kernel function . When represented as a GP, the loss function is
(15) 
and returns a Gaussian distribution when evaluated at a particular value . The kernel function describes the covariance between any two points in the GP and encodes structural properties of the underlying random function like the smoothness, periodicity, and stationarity. The kernel function is a nonlinear measure of similarity between inputs, with the dot product kernel leading to standard linear regression. In this study we use a common nonlinear kernel function, the Gaussian or squared exponential function:
(16) 
where describes the magnitude and describes the smoothness of the function or correlation length between two points.
The GP can be conditioned on a set of loss function training points to improve the posterior estimate of the loss function at new test points. The joint distribution of training outputs and a test output is
(17) 
Conditioning on the training data yields the following estimate for :
(18)  
(19)  
(20) 
where we have introduced the shorthand and . Note that the mean prediction is simply a linear combination of the training points and requires inverting a matrix that scales with the number of training points, N. Additionally, the posterior covariance is given by the prior covariance less the information gained from the training points.
IiiC Bayesian Quadrature
Recall that in stochastic economic dispatch we are faced with calculating the expectation of our GP loss function . In Bayesian quadrature [18, 19, 20] the integral is modeled as a Gaussian random variable because integration is a linear transformation that preserves Gaussianity and the integrand, , is a Gaussian process. The density function must also be a Gaussian, however this can be achieved by reweighting via importance sampling.
Bayesian quadrature allows us to express the expectation and variance of the loss analytically, and then devise a sampling strategy to optimally reduce this variance. Following the derivation in [18], the loss averaged over possible functions is simply the expectation of the mean GP loss function:
(21)  
(22)  
(23)  
(24) 
Similarly, the expected variance of the estimate with respect to possible functions is
(25)  
(26)  
(27)  
(28) 
A useful sampling strategy is to choose the set of samples that minimizes the expected variance:
(29) 
This can be compactly represented as
(30) 
where we have defined the following quantities
(31)  
(32) 
Because the integrated prior covariance has no dependence on the sample points, the optimization can be simplified to
(33) 
which shows that the optimal sampling strategy only depends on the choice of kernel covariance function and the nominal distribution , but not on the actual loss function evaluations . To develop some intuition for Eq. 33, we can see that maximization involves a balance between selecting points that reflect the nominal distribution while avoiding redundancies from selecting similar points.
That the sampling strategy is independent of the loss function evaluations is a key differentiator between importance sampling and Bayesian quadrature approaches. This is an advantage for Bayesian quadrature in cases where the loss function is unknown or stochastic, e.g. due to uncertainties in electrical loads. Furthermore, as a quadrature rather than sampling strategy, the Bayesian quadrature points can be reused at successive timesteps if and are unchanged. The Bayesian quadrature estimate of the expected loss function takes the form of the following quadrature rule
(34) 
Iv Performance of Different Sampling Strategies
Test data
To study the performance of our economic dispatch methods we examine a flattened version of the Reliability Test System, modified by the Grid Modernization Lab Consortium (RTSGMLC) [21]. The RTSGMLC system is an update to the original 1979 [22] and 1996 [23] version of the RTS that contains 73 buses, 96 conventional thermal generators and 4 wind plants. However, for our experiments we flatten the network and only consider power balance, i.e. matching the aggregated power output from the 96 conventional thermal generators and 4 wind plants to the aggregated demand without taking into consideration the network
The wind plant data used in our experiments is drawn from the Wind Integration National Dataset (WIND) Toolkit [24]. To appropriately model wind farms using WIND Toolkit data, wind plant power outputs were aggregated by first assigning geographic locations to nodes in the RTSGMLC system and then summing the output power from the wind plants nearest to the node locations until the prescribed annualized capacity for the node in the RTSGMLC system was reached.
The experiments in this section take place during a week in July. Load data from the RTSGMLC system was combined with simulated wind data from 2013 from the WIND Toolkit. To determine , we fit a Student’s tdistribution to errors in a 5 minute persistence forecast. A tdistribution was chosen because of the nonnormal, heavy tailed behavior of changes in wind generation. The observation heavytailed distributions provide a good fit to wind data has been noted previously, see, e.g. [25]. The distributions of wind forecast errors could be further refined by conditioning on current power, time of day, season, etc, however this study used a fixed tdistribution that best matched errors aggregated over the full month of July.
The actual loads and wind production data used in our experiments are shown in Figure 1. There are two important aspects of the data: first, the peak renewable generation is just over one quarter of the peak demand, corresponding to a moderate renewable energy penetration scenario. Second, there is a wind ramp event that occurs around midnight July 25th. The wind ramp event will pose an economic and reliability challenge to test the different sampling strategies in our experiments.
We examined the performance of each sampling strategy with 5, 10, 20, and 50 scenarios drawn at each 5 minute dispatch decision interval. The resulting costs of the first and second stage decisions for the full week are shown in Figure 3 for all three strategies and the four sets of scenarios. Additionally, a summary of the costs for the full week is shown in Table I. To make the differences in sampling strategies more clear, our model does not explicitly implement reserves. Instead any shortfall in generation results in an expensive lossofload penalty incurred during the second stage. This results in the bulk of the costs attributed to first stage dispatch decisions, and rare, but comparatively large second stage costs. In practice, these lossofload events would largely be addressed by ISO policies on reserve requirements.
Representative samples produced by each sampling technique, along with the nominal distribution, importance distribution, and loss function, are shown in Figure 2. The outcomes of these sampling strategies are shown in Figure 3 which displays a weeklong timeseries of first and second stage dispatch cost. The costs broken down by stage are additionally summarized in Table II and Table III.
The Monte Carlo method that draws from the nominal distribution has the most lossofload events. Alternatively, drawing from the importance sampling distribution leads to the smallest number of lossofload events. The reason for this behavior is clear: the shape of the loss function, large penalties for lossofload and relatively small penalties for overload or wind spilling, causes the importance distribution to favor sampling scenarios of wind overprediction. Thus, we are more likely to draw scenarios that represent lower than forecast values for wind from the importance distribution. These scenarios allow the stochastic optimization algorithm to hedge against expensive loss events. Drawing from the nominal distribution with Monte Carlo sampling, however, leads to scenarios that represent up and down ramps in wind power generation with equal probability. Therefore, for many of the 5 minute timeperiods, the stochastic optimization algorithm was not equipped to hedge against lossofload events with Monte Carlo samples.
The Bayesian quadrature approach results in a consistent decrease in maximum second stage costs as the number of scenarios is increased. This consistent performance improvement is a strength of a quadraturebased approach. In contrast, importance sampling does not show a monotonic decrease in peak second stage costs as more scenarios are considered. This is likely due to errors arising from approximating the true importance distribution.
Sampling Strategy Costs ($)  

# of Scenarios  MC  IS  BQ 
5  \num1.453e7  \num1.383e7  \num1.394e7 
10  \num1.409e7  \num1.382e7  \num1.394e7 
20  \num1.391e7  \num1.383e7  \num1.378e7 
50  \num1.386e7  \num1.382e7  \num1.378e7 
Sampling Strategy Costs ($)  

# of Scenarios  MC  IS  BQ 
5  \num1.355e7  \num1.363e7  \num1.356e7 
10  \num1.358e7  \num1.365e7  \num1.356e7 
20  \num1.360e7  \num1.366e7  \num1.360e7 
50  \num1.365e7  \num1.366e7  \num1.360e7 
Sampling Strategy Costs ($)  

# of Scenarios  MC  IS  BQ 
5  \num9.839e5  \num2.010e5  \num3.838e5 
10  \num5.101e5  \num1.688e5  \num3.838e5 
20  \num3.073e5  \num1.732e5  \num1.858e5 
50  \num2.099e5  \num1.648e5  \num1.834e5 
With only 5 or 10 scenarios to consider, the importance sampling approach yields the lowest total costs. We attribute this performance to the additional information available from the deterministic solution. Recall that an estimate of the loss function is required to determine the importance distribution. We approximate the loss function by solving a deterministic dispatch problem by assuming a perfect wind forecast. We then use the trapezoidal rule to calculate the normalizing factor that makes the importance distribution a proper pdf. Although this is a flawed estimate of the loss function, it is still valuable when only a few scenarios are available for optimization.
For 20 or more scenarios, the Bayesian quadrature approach yields a lower total cost. At this level of sampling, the Gaussian process approximation of the loss function is more accurate than the deterministic approximation made in importance sampling. We also note that the Bayesian quadrature approach is much less computationally expensive than importance sampling because it avoids solving an additional deterministic problem and integrating to normalize the importance distribution. Interestingly, the Bayesian quadrature solution yields a lower total cost than importance sampling despite triggering lossofload penalties more often. This is a reminder that the economically optimal dispatch is not necessarily the most reliable.
V Conclusion
This paper introduced two new techniques, importance sampling and Bayesian quadrature, to create scenarios of wind forecast error events for use in solving stochastic economic dispatch via a twostage recourse problem. The proposed techniques were compared against a Monte Carlo sampling strategy in a weeklong study of 5 minute economic dispatch using realistic generation, load, and wind forecast data from the RTSGMLC test case and NREL’s WIND Toolkit dataset. Both importance sampling and Bayesian quadrature achieve a reduction in estimator variance as compared to the baseline Monte Carlo sampling, which ultimately yield lower total costs for the twostage economic dispatch decisions. Furthermore, the scenarios selected by the improved strategies are drawn further out into the tails of the distribution of possible wind ramp events than with Monte Carlo sampling. The resulting dispatch decisions improve grid reliability and robustness by finding feasible solutions that reduce loss of load, while still maintaining economically optimal dispatch decisions.
We also found that the computational cost of determining sample locations is lower for Bayesian quadrature. Importance sampling suffers from the curse of circularity in that the normalization factor to make the importance distribution a proper probability density function requires knowing the expectation one is trying to estimate in the first place. To overcome this, we solve a deterministic economic dispatch problem to estimate the loss function and then numerically integrate it to approximate the normalization constant. This additional solve increases the computation cost of importance sampling relative to Bayesian quadrature. Furthermore, in the Bayesian quadrature approach the sample points can be computed once and reused at successive timesteps because it is a true quadrature approach rather than a stochastic sampling approach. The Bayesian quadrature samples only depend on the nominal wind pdf and the Gaussian process kernel, allowing them to be computed offline for wind distributions that are specific to a time of day, season, or current power output. Our analysis simply used a Student’s tdistribution to represent the distribution of wind ramps over the course of a year, however more conditionspecific distributions would likely improve performance even further.
The most economic dispatch decisions are found when using importance sampling with ten or fewer scenarios, and with Bayesian quadrature for more than ten scenarios. We attribute this to the extra information about the loss function that is gained through its deterministic approximation while computing the importance distribution. With more than ten scenarios, the Bayesian quadrature approach can develop a better loss function surrogate and find a lower cost dispatch strategy. While both proposed strategies improve the reliability from the baseline, the most economic strategy is not necessarily the most reliable. The solutions found via Bayesian quadrature for 20 and 50 scenarios produce the lowest total cost as shown in Table I, but as shown in Table III, they involve more recourse or second stage costs than the less economic strategies found via importance sampling, reflecting increased number of expected loss of load events.
Applying these strategies to a realistic test case of 5 minute economic dispatch with time varying loads, stochastic wind generation, and realistic costs yields improvements in total cost and in reliability represented by the second stage recourse costs. Depending on the number of scenarios used, total costs are reduced by  relative to the Monte Carlo baseline. A large component of the improvements in overall economic optimality is actually due to improvements in second stage costs that represent consumption of reserves or loss of load. The economically optimal dispatch strategy reduces these costs by  , demonstrating significant benefits to grid reliability and reserve requirements. Employing these advanced scenario generation strategies can significantly reduce costs and improve the reliability of operating the grid with high levels of stochastic renewable energy production.
Acknowledgment
This work was supported by the U.S. Department of Energy under Contract No. DEAC3608GO28308 with Alliance for Sustainable Energy, LLC, the Manager and Operator of the National Renewable Energy Laboratory. Funding provided by the Department of Energy’s Grid Modernization Lab Consortium and Exascale Computing Project.
References
 E. M. Constantinescu, V. M. Zavala, M. Rocklin, S. Lee, and M. Anitescu, “A computational framework for uncertainty quantification and stochastic optimization in unit commitment with wind power generation,” IEEE Transactions on Power Systems, vol. 26, no. 1, 2011.
 A. Soroudi and T. Amraee, “Decision making under uncertainty in energy systems: State of the art,” Renewable and Sustainable Energy Reviews, vol. 28, pp. 376–384, 2013.
 A. Shapiro, D. Dentcheva, and A. Ruszczyński, Lectures on stochastic programming: modeling and theory. SIAM, 2009.
 J. R. Birge and F. Louveaux, Introduction to stochastic programming. Springer Science & Business Media, 2011.
 A. Papavasiliou, S. S. Oren, and R. P. O’Neill, “Reserve requirements for wind power integration: A scenariobased stochastic programming framework,” IEEE Transactions on Power Systems, vol. 26, no. 4, 2011.
 J. Cheng, R. L.Y. Chen, H. N. Najm, A. Pinar, C. Safta, and J.P. Watson, “Chanceconstrained economic dispatch with renewable energy and storage,” Computational Optimization and Applications, vol. 70, no. 2, pp. 479–502, 2018.
 M. Lubin, Y. Dvorkin, and S. Backhaus, “A robust approach to chance constrained optimal power flow with renewable generation,” IEEE Transactions on Power Systems, vol. 31, no. 5, pp. 3840–3849, 2016.
 C. Safta, R. L.Y. Chen, H. N. Najm, A. Pinar, and J.P. Watson, “Efficient uncertainty quantification in stochastic economic dispatch,” IEEE Transactions on Power Systems, vol. 32, no. 4, 2017.
 A. Staid, J.P. Watson, R. J.B. Wets, and D. L. Woodruff, “Generating shortterm probabilistic wind power scenarios via nonparametric forecast error density estimators,” Wind Energy, vol. 20, no. 12, 2017.
 A. J. Wood and B. F. Wollenberg, Power generation, operation, and control. John Wiley & Sons, 2012.
 J. Zhu, Optimization of power system operation. John Wiley & Sons, 2015, vol. 47.
 J. A. Taylor, Convex optimization of power systems. Cambridge University Press, 2015.
 A. J. Conejo, M. Carrión, and J. M. Morales, Decision making under uncertainty in electricity markets. Springer, 2010, vol. 1.
 G. Giebel, L. Landberg, G. Kariniotakis, and R. Brownsword, “Stateoftheart methods and software tools for shortterm prediction of wind energy production,” in EWEC 2003 (European Wind Energy Conference and exhibition), 2003.
 A. M. Foley, P. G. Leahy, A. Marvuglia, and E. J. McKeogh, “Current methods and advances in forecasting of wind power generation,” Renewable Energy, vol. 37, no. 1, pp. 1–8, 2012.
 J. A. Bucklew, Introduction to Rare Event Simulation, ser. Springer Series in Statistics. New York, NY: Springer New York, 2004.
 C. E. Rasmussen and C. K. I. Williams, Gaussian processes for machine learning. Cambridge, MA: MIT Press, 2006.
 Z. Ghahramani and C. E. Rasmussen, “Bayesian Monte Carlo,” in Advances in Neural Information Processing Systems 15, S. Becker, S. Thrun, and K. Obermayer, Eds. MIT Press, 2003, pp. 505–512.
 A. O’Hagan, “BayesHermite quadrature,” Journal of Statistical Planning and Inference, vol. 29, no. 3, pp. 245–260, Nov. 1991.
 P. Hennig, M. A. Osborne, and M. Girolami, “Probabilistic numerics and uncertainty in computations,” Proc. R. Soc. A, vol. 471, no. 2179, p. 20150142, Jul. 2015.
 Grid Modernization Lab Consortium, “Reliability test system  Grid Modernization Lab Consortium,” 2017, [Online; accessed 11December2017]. [Online]. Available: https://github.com/GridMod/RTSGMLC
 P. Chairman, M. Bhavaraju, B. Biggerstaff et al., “IEEE reliability test system: a report prepared by the reliability test system task force of the application of probability methods subcommittee,” IEEE Trans Power Appar Syst, vol. 98, no. 6, pp. 2047–2054, 1979.
 C. Grigg, P. Wong, P. Albrecht, R. Allan, M. Bhavaraju, R. Billinton, Q. Chen, C. Fong, S. Haddad, S. Kuruganty et al., “The IEEE reliability test system1996. a report prepared by the reliability test system task force of the application of probability methods subcommittee,” IEEE Transactions on power systems, vol. 14, no. 3, pp. 1010–1020, 1999.
 C. Draxl, A. Clifton, B.M. Hodge, and J. McCaa, “The wind integration national dataset (wind) toolkit,” Applied Energy, vol. 151, 2015.
 B.M. Hodge and M. Milligan, “Wind power forecasting error distributions over multiple timescales,” in Power and Energy Society General Meeting, 2011 IEEE. IEEE, 2011, pp. 1–8.