A critical step in stochastic optimization models of power system analysis is to select a set of appropriate scenarios and significant numbers of scenario generation methods exist in the literature. This paper develops a clustering based scenario generation method, which aims to improve the performance of existing scenario generation techniques by grouping a set of correlated wind sites into clusters according to their cross-correlations. Copula based models are utilized to model spatiotemporal correlations and the Gibbs sampling is then used to generate scenarios for day-ahead markets. Our results show that the generated scenarios based on clustered wind sites outperform existing approaches in terms of reliability and sharpness and can reduce the total computational time for scenario generation and reduction significantly. The clustering-based framework can therefore provide a better support for real-world market simulations with high wind penetration.

## I. INTRODUCTION

Stochastic programming represents a potentially promising technique to address uncertainties associated with wind power^{1–3} and has been studied widely in the recent decades^{4–6} in a variety of fields.^{7–10} However, one challenge in its practice is to select a set of appropriately weighted scenarios to represent the space of uncertainty. Such techniques usually involve fitting forecasted wind power or forecast errors to specific distributions,^{11} where scenarios are then generated by sampling the derived distributions.^{12} Such distributions can be either parametric distributions, such as Gaussian^{5} and Generalized Gaussian distributions,^{13} the Beta distribution,^{14,15} the Weibull distribution,^{16} and the logit-normal distribution,^{12} or non-parametric distributions.^{17} Ma *et al.*^{18} characterized forecast errors using empirical distributions and the inverse transformation method is applied to obtain a set of scenarios. Cui *et al.*^{13} developed a generalized Gaussian mixture model to fit forecast errors of aggregated wind power from hundreds of wind farms and the fitted distribution is used to sample scenarios for probabilistic wind ramp forecasting.

While sampling methods can be applied to these distributions to generate scenarios for a single wind farm or aggregated power output of multiple wind farms,^{5,11,13,15,17,18} an increasing number of studies have placed emphasis on spatiotemporal correlation in scenario generation.^{19–21} Typically, such correlation is modeled with multivariate joint distributions. For example, Pinson *et al.*^{12} employed a multivariate Gaussian distribution to describe correlations between wind power forecasts made at different lead times, and this method has been widely adopted in numerous studies recently. Other commonly used methods include Cholesky decomposition,^{22} moment matching techniques,^{23} and the generalized dynamic factor model.^{24} In addition, machine learning-based methods are also employed to generate scenarios in the latest studies. Such methods include artificial neural network (ANN),^{25} quantile regression convolutional neural network (QRCNN),^{26} radial basis function neural network (RBFNN),^{18} generative adversarial networks (GAN),^{27,28} and support vector machine (SVM).^{29} Particularly, techniques that are used to generate probabilistic forecasts can be used to produce discrete scenarios by sampling the predictive distributions to avoid reliance on prior distributions.^{30} Nevertheless, modeling high dimensional multivariate non-Gaussian distributions can be challenging, and a common approach is to use copula.^{31} By applying the marginal cumulative distribution functions to stochastic variables, the original variables are transformed from the original space into a common uniform domain, in which correlations among the original variable can be further characterized using copulas. The copula-based methods have been used extensively in the literature recently.^{12,32–35} For example, Zhang *et al.*^{19} modeled spatiotemporal correlations of clustered wind farms using a copula based model and implemented a scenario generation method. A similar method is developed by Tang *et al.*^{20} and the Gibbs sampling method is adopted to generate scenarios for a stochastic economic dispatch (ED) model.

The above methods represent a considerable number of scenario generation approaches. However, real-world power systems usually involve hundreds or even thousands of wind farms, while sampling random variables of such high dimensions is usually computationally prohibitive, especially when spatial-temporal correlation is explicitly considered. Recently, data mining techniques have been applied to power system studies to process large volumes of data in a diverse set of literature studies.^{36,37} Typical applications include cable layout design, evaluation of transfer capabilities, reliability of networks, and also wind power forecasting.^{38} In the study by Vallee *et al.*,^{37} a clustering method is developed to group 94 wind areas into 14 clusters, which are then used to generate samples of wind speed to evaluate system reliability through a non-sequential Monte Carlo simulation.

In this paper, we present a scenario generation framework based on cluster analysis, where wind sites are grouped into clusters by correlation. The framework is aimed to improve the quality of generated scenarios at less time expense by maximizing intra-cluster correlation and minimizing inter-cluster correlation. By clustering wind sites into a number of smaller clusters, existing methods can be employed to generate scenarios for each cluster. The generated scenarios of wind power are then reduced and used as inputs to a stochastic day-ahead unit commitment (DAUC) model, which uses a hedging strategy to clear the day-ahead market to determine commitment statuses of thermal units. Given the commitment statuses of thermal generators, a real-time ED model is then solved to evaluate the system cost-effectiveness and reliability. A synthetic, 2000-bus network built on the footprint of the Electric Reliability Council of Texas (ERCOT)^{39} is employed to demonstrate the applicability of the proposed method in real-world power systems. Our approach represents an enhanced framework based on existing methods and can be applied to real-world power systems with high wind penetration.

The main contributions of our work include (i) enhancing existing scenario generation methods by substantially improving the quality of generated scenarios with reduced time consumption, (ii) demonstrating the proposed framework with a real-world sized power system and scrutinizing the trade-off between the number of clusters and a variety of evaluation metrics, and (iii) using the generated scenarios as inputs to a stochastic power market simulation of a synthetic network built on the footprint of the Texas power system and assessing its impact in market operation and scheduling. The proposed framework can be extended by integrating with other scenario generation methods and applied to real-world power systems with high wind penetration.

The remainder of this paper is organized as follows: Sec. II details the developed clustering-based framework. Section III formulates the power market simulation model. Section IV summarizes the data and assumptions. Sections V and VI present the generated scenarios and the power market simulation, respectively. Section VII concludes this paper.

## II. CLUSTERING BASED SCENARIO GENERATION

### A. Cluster analysis

While numerous studies suggest that including correlations across wind farms can provide better support for existing scenario generation methods,^{19–21} these studies typically include all wind farms in a single group. However, wind farms in a real-world power network usually scatter over a wide area since geographical diversification can effectively decrease fluctuations in aggregate output and improve overall economic viability.^{40,41} In addition, correlations normally decay with distances, which results in a weak correlation between two wind farms if they are remotely located. Therefore, it is questionable whether the generated scenarios can benefit from including weakly correlated wind farms and a more appropriate approach is to group wind farms into clusters based on their cross-correlations.^{42}

Cluster analysis groups data objects based only on information provided by the data itself, where data objects in one cluster often share similar characteristics and are distinct from other clusters. The similarity or distance between data objects is usually expressed in terms of proximity measures, such as $\u2113p$ norms, or the cohesion of the cluster. In this study, we use Pearson's correlation coefficient $\rho (xi,xj)$ as the proximity measure since the wind farms are clustered based on correlation, i.e., wind farms in the same cluster should be more closely correlated with each other, while correlations between two clusters are less strong. The formula for $\rho (xi,xj)$ is given by the equation below:

A variety of methods exist in cluster analysis.^{43} In this study, we use the *k*-means method. The *k*-means method is a widely used prototype-based, partitional method that selects the mean of all data objects in a cluster as the prototype. Let *C* denotes a set of data objects, which can be grouped into *K* clusters and let $x\u2208Ck$ denotes the data objects in the *k*th cluster, the objective function of the *k*-means method becomes

where $dist()$ represents the selected proximity measure. Note that since higher *ρ* represents stronger correlation, hence higher similarity, the objective function minimizes $1\u2212\rho $ instead:

To solve the optimization problem, one can exhaustively enumerate all possible ways of grouping objects into *K* clusters. However, this is often computationally prohibitive and heuristic iteration-based algorithms are typically used. Such heuristic algorithms often require initial conditions (i.e., centroids) as input and iteratively converge to a solution where the objective value no longer decreases. Although the heuristic algorithms can converge within an acceptable amount of time, they are sensitive to initial conditions and often converge to a local minimum rather than a global minimum. Therefore, multiple runs with randomly selected initial conditions are typically required. In this paper, the k-means function in MATLAB is employed and repeated ten times using different initial cluster centroids.

### B. Scenario generation

After all wind sites are grouped into *K* clusters, existing scenario generation methods can be applied to each cluster. Given a cluster with *N* wind sites, the spatial correlations of wind power production across all sites in the cluster can be modeled using a multivariate joint probability distribution, whose cumulative density function (CDF) and probability density function (PDF) are given below:

where $Xi,ta$ and $Xi,tf$ denote actual and forecasted wind power at site *i* during time *t*, respectively. Note that the forecasts can be deterministic forecasts with any lead time and in our analysis, day-ahead forecasts are used. Therefore, scenario generation given wind power forecast becomes sampling the distribution of actual wind power $(X1,ta,\u2026,XN,ta)$ conditioned on the forecasted wind power $(x1,tf,\u2026,xN,tf)$,

Here, we follow Tang *et al.*^{20} by modeling the multivariate distribution with marginal distributions of individual variables and a copula to link their interdependence. The spatial correlation is modeled with the copula, while the temporal correlation is modeled with a multivariate Gaussian distribution. Next, Gibbs sampling is adopted to simulate sampling of the multivariate random variable by sequentially sampling each component.

#### 1. Gibbs sampling

The Gibbs sampling method is a Markov chain Monte Carlo algorithm^{44,45} and the idea is to iteratively sample only one variable or a block of variables at a time from its distribution conditioned on the remaining variables. Suppose we have a two-dimensional random variable (*X*, *Y*), the Gibbs sampler generates a sequence (known as a Gibbs sequence) by sampling each component conditioned on the other component:

The initial value $y(0)$ is specified and the rest of the sequence is obtained by iteratively sampling,

The process of generating a Gibbs sequence is Gibbs sampling. It shows that under certain conditions, the distributions of $X(k)$ and $Y(k)$ converge to *f*(*X*) and *f*(*Y*) (i.e., their marginal distributions), respectively, as $k\u2192\u221e$. Without knowing the joint distribution, we can generate a set of samples that follow their marginal distributions. It can also be extended to high dimensional space as well. Therefore, Gibbs sampling converts sampling a multivariate distribution into sequentially sampling a set of univariate distributions, i.e., sampling the multivariate distribution in Eq. (5a) becomes sampling $Fai|a1,\u2026,ai\u22121ai+1,\u2026,\u2009aNf1,\u2026,fN$ iteratively. A detailed procedure is given in Algorithm 1. Note that, in this analysis, we use the forecasted wind power $xi,tf$ to initialize $xi,ta$.

procedure Gibbs sampler |

Initialize $xi,ta,0\u2190xi,tf,\u2200i=1,\u2026,N$ |

for scenario $\xi \u2208\Xi $ do Sample: |

$X1,ta|x2,ta,\xi \u22121,\u2026,xN,ta,\xi \u22121,x1,tf,\u2026,xN,tf$ |

$X2,ta|x1,ta,\xi ,x3,ta,\xi \u22121,\u2026,xN,ta,\xi \u22121,x1,tf,\u2026,xN,tf$ |

… |

$XN,ta|x1,ta,\xi ,\u2026,xN\u22121,ta,\xi ,x1,tf,\u2026,xN,tf$ |

end for |

end procedure |

procedure Gibbs sampler |

Initialize $xi,ta,0\u2190xi,tf,\u2200i=1,\u2026,N$ |

for scenario $\xi \u2208\Xi $ do Sample: |

$X1,ta|x2,ta,\xi \u22121,\u2026,xN,ta,\xi \u22121,x1,tf,\u2026,xN,tf$ |

$X2,ta|x1,ta,\xi ,x3,ta,\xi \u22121,\u2026,xN,ta,\xi \u22121,x1,tf,\u2026,xN,tf$ |

… |

$XN,ta|x1,ta,\xi ,\u2026,xN\u22121,ta,\xi ,x1,tf,\u2026,xN,tf$ |

end for |

end procedure |

#### 2. Spatial correlation

To obtain the univariate distributions required in Algorithm 1, we start from the conditioned probability distributions from Eq. (4a). By applying Bayes' theorem, the CDF of $Xi,ta$ conditioned on other variables in Algorithm 1 is

By applying the copula function, the equation becomes

which only includes a copula function $C2N$ and marginal distribution functions associated with each variable $Fai$. Note that detailed derivation can be found in Appendix.

Once the univariate CDF is determined, Algorithm 1 is implemented sequentially to produce samples. In this study, the logit-normal distribution is used to fit the marginal distributions of wind power,^{11} and the copula is selected from among Gaussian, student *t*, and the Archimedean family based on the Bayesian information criterion (BIC).^{13}

#### 3. Temporal correlation

Samples produced from the above method represent simultaneous power production from multiple spatially correlated wind farms. However, as suggested by Pinson *et al.*,^{12} power production from a wind farm at different times are also temporally correlated, which can be modeled as follows: Given wind farm *i* and let $ui,ta=Fai(xi,ta)\u2208(0,1)$ denote the CDFs associated with its actual power production at *t* = 1 to *T*, suppose $Yi=(Yi,1,Yi,2,\u2026,Yi,T)$ is a *T*-dim variable whose entries are given by $Yi,t=\Phi \u22121(ui,ta)$, then, $Yi$ can be modeled using a multivariate Gaussian distribution

where $\Phi \u22121$ is the inverse of the standard Gaussian CDF function, $0$ is a *T*-dim vector of zeros, and $\Sigma i$ is a covariance matrix whose entries are given by

Therefore, Eq. (8) implies that the temporal correlation decays exponentially with time and is dictated by the parameter *v*. In this study, the value of *v* is selected by fitting $\Sigma i$ to historical data using the method of least squares.

### C. Scenario reduction

After $|\Xi |$ scenarios are generated for each of the *K* clusters, the generated scenarios can be directly used as inputs to the stochastic programming models. However, the exhaustive combination of all scenarios from *K* clusters results in $|\Xi |K$ scenario, which can be computationally prohibitive as *K* increases. Therefore, the number of generated scenarios must be reduced to an adequate number before the model can be solved in a time efficient manner. Numerous scenario reduction techniques exist in the literature,^{46–49} where similar scenarios are aggregated based on particular distance measures. However, most of these reduction techniques do not involve clusters. In addition, the process of scenario generation itself can also be time-consuming, particularly given the great number of scenarios if the scenarios in all clusters are exhaustively combined. Therefore, we adopt the method of forward reduction from Dupacova *et al.*^{49} and apply it in a heuristic way to reduce the number of scenarios from $|\Xi |K$ to 10.

For convenience, we define a *T*-dim vector $xka,\xi $ to represents scenario *ξ* in the *k*th cluster,

where *N _{k}* is the number of wind sites in the

*k*th cluster. Therefore, each scenario in a cluster is represented by the sum of wind power from all wind sites in that cluster. Let's denote the original set of sums as $Sk={xka,\xi ,k=1,\u2026,K,\xi \u2208\Xi}$ and use

**S**as the reduced set. Using the above definition, the detailed implementation is given in Algorithm 2.

procedure Scenario reduction |

Initialize $S\u2190$ empty set $\u2205$ |

for $k=1,\u2026,K$ do |

Add scenarios from cluster k: $S\u2190S\xd7Sk$ |

Reduce S into 10 scenarios using forward reduction |

end for |

Return the final scenarios $S={xa,j,j=1,\u2026,10}$ |

end procedure |

procedure Scenario reduction |

Initialize $S\u2190$ empty set $\u2205$ |

for $k=1,\u2026,K$ do |

Add scenarios from cluster k: $S\u2190S\xd7Sk$ |

Reduce S into 10 scenarios using forward reduction |

end for |

Return the final scenarios $S={xa,j,j=1,\u2026,10}$ |

end procedure |

Note that during the *k*th iteration, a new cluster $Sk$ is added and **S** is updated by taking the Cartesian product of **S** and $Sk$, i.e., the element in **S** is updated by taking all ordered pairs $(xa,j,xka,\xi )$, where $xa,j\u2208S,j=1,\u2026,10$ and $xka,\xi \u2208Sk,\xi \u2208\Xi $. Therefore, the dimension of set **S** increases by *T* after each iteration and after all iterations, $dim\u2009S=K\xb7T$. Upon completion of Algorithm 2, each single element in set **S** includes the sums of wind power output from all clusters. Based on the indices of the selected scenarios in **S**, site-specific wind power output can then be picked up and used in the power market simulation.

### D. Performance evaluation

To evaluate the performance of the generated scenarios, two prediction interval (PI) based metrics are introduced: reliability and sharpness. These metrics are used widely in probabilistic forecast to assess the quality of the PIs.^{50} Probabilistic forecasts can take a variety of forms,^{51,52} such as PDFs, quantiles, and scenarios. While the evaluation metrics are often defined based on PDFs, other forms can be converted into PDFs as needed. Since the generated scenarios can be viewed as realizations of probabilistic wind power forecasts, the PIs are estimated based on the generated scenarios and the PI-based metrics are then adopted to evaluate the quality. In addition, the continuous ranked probability score (CRPS) is also employed to measure the quality of our scenarios.

#### 1. Prediction interval

Given the forecasted mean *x* and variance $\sigma 2$ of a random variable *X* probabilistically forecasted random variable *x* and its variance $\sigma 2$, a PI with $100(1\u2212\alpha )%$ confidence level can be given by $[L\alpha ,U\alpha ]$,

where $z1\u2212\alpha /2$ is the $100(1\u2212\alpha /2)$ percentile of the standard Gaussian distribution. Note that the random variable *X* here can be the wind power production at site *i* in time *t*, i.e., $Xi,ta$, while *x* and $\sigma 2$ in Eqs. (10) and (11) are the estimated first and second moments from the generated scenarios. Since the PI of normalized wind power is bounded in $[0,1]$, anything out of the interval will be floored at 0 and capped at 1, i.e., the resulting distribution is assumed to be censored Gaussian. Note that a previous study^{50} suggests that even if the actual distribution is non-Gaussian, the assumption of a censored Gaussian can still be applied with satisfactory performance.

#### 2. Reliability

Given confidence level $100(1\u2212\alpha )%$, PIs can be derived from the generated scenarios for each random variable. The future measured wind power is expected to lie within the PI. The given confidence level is also known as the PI nominal confidence (PINC), and it is expected that the probability of the obtained PIs covering the measured wind power $xi,ta$ will asymptotically reach the nominal level of confidence. Therefore, the reliability of the PIs is measured by PI coverage probability (PICP),

where the indicator of PICP, $Ii,t\alpha $, is defined by

Ideally, the PICP should equal the PINC, i.e., the average coverage error (ACE)^{50} should be as close to zero as possible. In our analysis, we use the average absolute coverage error over the entire forecast time horizon *T* and across all *N* wind sites as the indicator of reliability,

where $PINC=100(1\u2212\alpha )$.

#### 3. Sharpness

The definition of PICP implies that a high PICP can be realized by widening PIs. Therefore, another metric, sharpness, is introduced to measure the size of PIs. Sharpness is evaluated in the form of interval score. Given the width of a PI: $\omega i,t\alpha =Ui,t\alpha \u2212Li,t\alpha $, the interval score $Sci,t\alpha $ is

where $Sci,t\alpha $ rewards a narrow PI and gives penalty if the PI fails to cover the realized value $xi,ta$.

Similar to the definition of ACE, we use an average interval score (AIS) to evaluate the sharpness of the generated scenarios over time horizon *T* and across all *N* wind sites,

#### 4. Overall metric

Clearly, there is a trade-off between reliability and sharpness, since a wider PI may result in a higher coverage rate at the cost of a lower interval score, and *vice versa.* We, therefore, introduce a synthetic evaluation metric (SEM) as the linear combination of ACE and AIS to give a comprehensive evaluation of the PI,

Note that a lower SEM indicates a better overall performance since a negative AIS is used here. *λ*_{1} and *λ*_{2} are weighting coefficients, and we assume $\lambda 1=\lambda 2=0.5$ to put equal emphasis on both reliability and sharpness. The weighting coefficients can be varied by decision makers to put more weights on either aspect.

#### 5. Score

There are numerous numerical scores proposed to measure the quality of probabilistic forecasts.^{53} We use the CRPS score, a proper score that is designed to measure both the reliability and sharpness of a probablistic forecast.^{54} Suppose the predicted CDF of the wind power production in site *i* at time *t* is $F\u0302i,t(x)$ and the observed wind power is $xi,ta$, the definition of CRPS can be given by the following equation by measuring the difference between the predicted and observed CDFs,

where $1x\u2265xi,ta:{x|x\u2208\mathbb{R}}\u2192{0,1}$ is the indicator function of event $x\u2265xi,ta$. We further take the average over all wind sites and all time intervals to measure the overall performance,

Note that the predicted CDF is required before calculating the CRPS. In the realm of weather prediction, the predicted CDF can be obtained by three methods: parametric method, non-parametric method, and ensemble forecast.^{51} By giving a set of discrete scenarios, our study falls in the category of ensemble forecast. Therefore, to estimate the predicted CDF, the empirical CDF is derived from the generated scenarios by using the following equation:

where $xi,ta,\xi $ denotes the generated wind power scenarios for wind site *i* at time *t*.

## III. POWER MARKET SIMULATION

A critical motivation of this study is to apply the enhanced approach of scenario generation to real-world power systems with high renewable penetration. Power markets in real world are typically operated on multiple timescales ranging from days ahead to real time. Typically, as the markets evolve forward, forecasts with shorter lead times are used as they become available. In most real-world power markets, the day-ahead market is often cleared one day before using day-ahead forecasts and the schedules are updated in the real-time market several hours ahead of the actual operation by using more up-to-date forecasts. The imbalance between the forecasts and the actual values is then offset by deploying operating reserves, which is known as the load frequency control (LFC) process. The purpose of the power market simulation is to show the market benefits from the improvement in forecasts due to clustering. However, we believe that including forecasts with multiple lead times would be unnecessarily complex to illustrate the benefits of the methods utilized in this work. Therefore, although forecasts with different lead times are used in real-world markets, we only focus on day-ahead forecasts and assume that the real-time market is cleared based on actual values to neglect the LFC process. In addition, we assume that both markets are operated on a 24-h horizon at an hourly resolution. To simulate the two markets, we develop a day-ahead unit commitment (DAUC) model and a real-time economic dispatch (RTED) model. The DAUC model takes in day-ahead wind power forecasts and determines the commitment statuses of all thermal units. Next, the RTED model is solved based on actual wind power production to determine the final financial settlement. In addition, the DAUC model is solved both deterministically and stochastically. The stochastic run solves for a hedging strategy that explicitly considers future uncertainties associated with wind power, which is represented by scenarios. By contrast, the deterministic run solves each wind scenario independently and all decision variables vary by scenarios. Once the DAUC model is solved, commitment statuses of all thermal units are fixed and the RTED model is solved by using the actual wind output. We then compute the financial settlement only based on the solutions from the RTED model.

### A. The deterministic DAUC model

The DAUC model is formulated as a mixed integer program by following Ref. 55 and enhanced by constraints associated with regulating and responsive reserves. In addition, active power flow on transmission lines is calculated using shift factors and bounded by a set of pre-defined line limits. The objective of the model is to minimize total costs, which include actual production costs, fixed costs, and load and curtailment costs. The full DAUC model formulation is given below. Note that decision variables in each scenario (*ξ*) are independent of other scenarios in the deterministic runs.

System load balance

Maximum and minimum power generation

Note that if wind unit $g\u2208GW$ is installed at site *i*, the maximum available wind power $Pg,tF,\xi $ is given by multiplying the normalized wind power $xi,ta,\xi $ with the rated capacity of unit *g*.

Maximum available power of thermal generators accounting for ramp rates

Power generation during ramp-down or shut-down

Active power flow, $\u2200\u2113\u2208L$,

Up-time of thermal generators, $\u2200g\u2208GT$,

Down-time of thermal generators, $\u2200g\u2208GT$,

Reserve availability, $\u2200g\u2208GT$,

Spinning reserve balance, $\u2200t\u2208T$,

Regulating up/down reserve balance, $\u2200t\u2208T$,

Cost of thermal generators, $\u2200g\u2208GT,t\u2208T$,

Curtailment cost, $\u2200t\u2208T$,

Constraint (20) ensures system-wide load balance. Different from Ref. 55, we include slack variable $\delta b,t\xi $ for each load bus $b\u2208BL$ to reduce the chances of infeasibilities. Equations (21a) and (21b) allow thermal units to operate within their output limits, while the output of renewable units is limited by power forecast in Eq. (21c). In addition, Eqs. (22a) and (22b) limit the maximum output from thermal generators during start-up and shut-down. The output of thermal units is further limited by Eq. (23) during ramp-down or shut-down. In Eq. (24), the active power flow on transmission line $\u2113$ is calculated using shift factors $SF\u2113,b$ and is limited by line limit $\Lambda \u2113$. The minimum online/offline time of thermal generators are limited by Eqs. (25) and (26). The amount of available reserves ($rg,tS,\xi $ for responsive reserve, $rg,tRU,\xi $ and $rg,tRD,\xi $ for up/down regulating reserves) from thermal units is provided in Eq. (27) and requirements of reserves are met in Eqs. (28) and (29). Similar to the load balance constraint, slack variables are also introduced to account for reserve shortages.

Costs of thermal generators consist of production costs and fixed costs. Typically, production costs of thermal generators are represented by quadratic functions; here for simplicity, they are modeled as linear functions of power production $pg,t\xi $, as shown in Eq. (30a). The fixed costs, as given by Eqs. (30b) and (30c), include start-up and shut-down costs, which are determined solely by commitment statuses. In addition to the costs of thermal units, the objective function also includes curtailment costs, which are used to avoid load shedding and reserve shortages, as given in Eq. (31). Finally, the objective function in the deterministic run is given in Eq. (32) by summing up all costs:

### B. The stochastic DAUC model

To account for uncertainties associated with wind power, a two-stage stochastic DAUC model is developed based on the deterministic DAUC model. We assume that the commitment statuses ($vg,t\xi $) are the only first-stage decision variables and must remain equal across all scenarios, while all other variables are the second-stage decision variables and are thus scenario-dependent. Therefore, a non-anticipativity constraint is introduced in the following equation:

The objective function of the stochastic run is given by Eq. (34), where the probability-weighted costs are minimized. The fixed costs are equal across all scenarios, since they are only dependent on the first stage variables. By contrast, the production costs and curtailment costs are functions of the second stage variables, which are scenario dependent; therefore, Eq. (34) minimizes their probability-weighted expectation. Note that probabilities of scenarios ($Pr\xi $) are determined in the scenario reduction process

### C. The RTED model

After the DAUC model is solved, the RTED model is solved by fixing the commitment status of all thermal units and replacing the forecasted wind power with actual values. In addition, the minimum online/offline time constraints in Eqs. (25) and (26) are also removed. Therefore, the RTED model is formulated as a continuous linear program, due to the removal of all binary variables. The full model consists of Eqs. (20)–(24) and (27)–(31), and the objective function is represented by Eq. (32). Note that while each deterministic scenario results in an independent set of commitment status, hence an independent RTED run, the stochastic run returns only one set of commitment status for all scenarios and is followed by one single RTED run.

## IV. DATA AND ASSUMPTIONS

To simulate real-world power systems, we use a well-developed, publicly available synthetic network developed by Ref. 39, which is built on the footprint of the Electric Reliability Council of Texas (ERCOT). The 24-h load profiles from the IEEE 118 bus system are mapped into load profiles of the Texas system. Note that the load profile of one load bus in the 118 bus system may be mapped into multiple load buses in the Texas system, due to the latter's significantly larger size. In addition, the load profiles are scaled up or down to match the original load magnitudes of the Texas system.

Generating units in the Texas system can be grouped into six categories based on their fuel types, as shown in in Table I. Most techno-economic parameters of thermal units, such as the maximum and minimum output and ramp rates, are directly drawn from the original dataset. The production cost functions are mapped from the units with similar fuel types and rated capacities in the IEEE 118 bus system. The minimum online and offline time of a thermal units is determined primarily by its prime driver, while the start-up and shut-down costs are also affected by fuel costs. Since the original dataset does not give the types of prime drivers, the prime drivers are determined by their capacities following Ref. 56. Fuel costs are drawn from US EIA's Monthly Energy Review.^{57} Note that we assume that all nuclear units are constantly online. In addition, wind and solar units have zero production and start-up/shut-down costs.

. | Synthetic ERCOT system . | Real-world Texas system . | |||
---|---|---|---|---|---|

Fuel type . | No. of units . | Capacity (MW) . | Share . | Capacity (MW) . | Share . |

Coal | 39 | 14 502 | 15% | 23 589 | 19% |

NG | 367 | 63 810 | 66% | 69 386 | 56% |

Nuclear | 4 | 5139 | 5% | 4960 | 4% |

Wind | 87 | 9587 | 10% | 22 583 | 18% |

Solar PV | 22 | 651 | 1% | 1240 | 1% |

Hydro | 25 | 2603 | 3% | 670 | 1% |

Biomass | 0 | 0 | 0% | 425 | 0% |

Other | 0 | 0 | 0% | 659 | 1% |

Total | 546 | 96 292 | 100% | 123 512 | 100% |

. | Synthetic ERCOT system . | Real-world Texas system . | |||
---|---|---|---|---|---|

Fuel type . | No. of units . | Capacity (MW) . | Share . | Capacity (MW) . | Share . |

Coal | 39 | 14 502 | 15% | 23 589 | 19% |

NG | 367 | 63 810 | 66% | 69 386 | 56% |

Nuclear | 4 | 5139 | 5% | 4960 | 4% |

Wind | 87 | 9587 | 10% | 22 583 | 18% |

Solar PV | 22 | 651 | 1% | 1240 | 1% |

Hydro | 25 | 2603 | 3% | 670 | 1% |

Biomass | 0 | 0 | 0% | 425 | 0% |

Other | 0 | 0 | 0% | 659 | 1% |

Total | 546 | 96 292 | 100% | 123 512 | 100% |

For simplicity, we assume that outputs from all solar and hydro units are zero. Neglecting solar and hydro power may affect the commitment statuses of thermal units; however, we do not think it introduces a systematic error given their small capacity shares. As can be seen from Table I, both hydro power and solar PV only contribute marginally to the total capacity in both the real-world and the synthetic systems. In addition, it allows us to focus on the uncertainties associated with wind power.

We assume that load curtailment cost is $9000/MWh, according to the value of load loss (VOLL) used by ERCOT.^{58} We use $2000/MWh, the low system-wide offer cap (LCAP) provided by the Texas Public Utility Commission,^{59} as the penalty price for responsive reserve shortage and $5500/MWh, the average of LCAP and VOLL, as the penalty price for the regulating up/down reserve shortages. The penalty prices are selected such that responsive reserves are curtailed ahead of regulating reserves, and load will be the last to be curtailed.^{60} Historical data from ERCOT indicate that the regulating up/down reserves account for 1%–2% of system-wide load, while the responsive reserves range from 5% to 11%. Therefore, we assume 10% for responsive reserve margin and 2% for both up and down regulating reserve margins.

The data drawn from NREL's Wind Integration National Dataset (WIND) toolkit^{62} are used as input to the scenario generation model. The WIND toolkit provides synthetic wind turbine output power as well as meteorological data for more than 126 000 sites in the contiguous United States from 2007 to 2012. The data from 2007 to 2011 are used as input to the cluster analysis. The WIND toolkit simulates wind power production every 5 min and applies an NWP model to provide hourly power forecasts at multiple lead times. In our study, the day-ahead forecast is selected to generate scenarios for the day-ahead market. Since the temporal granularity is every 5 min for the actual power data and hourly for all forecasts, we calculate the hourly average of the actual wind power in the training dataset to match the time resolution of the day-ahead forecast. To capture seasonal variations of load profiles in the power market simulation, one day is selected from each season in 2012: February 10, 2012 (winter), May 10, 2012 (spring), August 10, 2012 (summer), and November 10, 2012 (fall). The Gibbs sampler is trained using actual wind power and wind power forecast from January 1, 2007 to the day before each target day. Next, scenarios are generated by the Gibbs sampler based on the day-ahead forecasts for the target day. The 87 wind units in the Texas system are mapped from 45 sites (shown in Fig. 1) in the WIND toolkit based on their geographical coordinates. Therefore, all wind units in one site have identical normalized wind power, although they may have different rated capacities. Figure 2 shows a heat map of the cross-site correlation coefficients across all 45 wind sites in this study, indicating strong correlations.

## V. RESULTS: CLUSTERING-BASED SCENARIO GENERATION

In this section, the clustering based framework developed in Sec. II is applied to the 45 wind sites from the Texas system over the four target days. The wind sites are first divided into a given number (*K*) of clusters and the Gibbs sampler is applied to each cluster independently. Since Gibbs sampling requires a large number of iterations before the samples reach a stationary distribution, 9000 samples are generated for each cluster and the first 8000 samples are discarded (known as the “burn-in” process) and only the last 1000 samples are kept for market simulation. The performance metrics are then derived from the generated scenarios and scrutinized. The goal is to examine the relation between the number of clusters and the quality of generated scenarios. The optimal number of clusters is jointly determined by the performance metrics and time consumption. In this study, the number of clusters *K* ranges from 1 to 20. Note that since when *K* = 1, all wind sites are put into one single group, it is equivalent to the case where cluster analysis is not applied. Therefore, the case when *K* = 1 represents our baseline case. The implementation is coded in MATLAB and executed on a high performance computing platform to exploit its parallel computing capability. A compute node equipped with a 20-core Xeon E5–2698 CPU and 256 GB memory is used, where the model training and scenario generation are both parallelized.

Figures 3 and 4 show the results from an ex-post validation analysis by scrutinizing the autocorrelation coefficients and cross correlation coefficients. The autocorrelation coefficients of historical wind power from site *i* are given by the following equation:

where $\Delta t$ is time lag and $xi,ta$ and $xi,t+\Delta ta$ are taken from the training dataset. Since the generated scenario $xi,ta,\xi $ can be viewed as an estimator of $xi,ta$, the estimated autocorrelation coefficient $\rho \u0302i(\Delta t)$ can be given by the above equation when $xi,ta$ is substituted by $xi,ta,\xi $. The autocorrelation is caused by temporal correlations in the wind power time series. By comparing $\rho i(\Delta t)$ and $\rho \u0302i(\Delta t)$ from all 45 sites, Fig. 3 indicates good agreement between the model input and output. The declining trends indicate weaker temporal correlations as time lag increases, as modeled in Eq. (8).

As presented in Fig. 2, the wind sites are strongly correlated. As spatial correlation is modeled by copula functions, the resulting scenarios should present strong correlations across sites. We examine the linear correlation coefficients using the generated scenarios and compare it with the original coefficients from historical data. The correlation coefficient is given by

Similar to the autocorrelation coefficients, the estimated cross correlation coefficients $\rho \u0302i,j$ are given when $xi,ta$ are substituted by $xi,ta,\xi $. Due to the large number of sites, we calculate the relative differences between $\rho i,j$ and $\rho \u0302i,j$,

Figure 4 visualizes the relative differences as heat maps, where each pixel represents one pair of sites. By comparing it with Fig. 2, it shows approximation errors using the copula functions. However, it also indicates that while greater errors are presented when the spatial correlations are weak, most errors for pairs with strong spatial correlations are not significant, indicating good agreement between the model input and output.

The results from the winter day (February 10, 2012) are shown in Fig. 5. As illustrated in Fig. 5(a), the diagonal indicates the ideal results where PICP equals PINC; therefore, greater distance between the PICP and the diagonal implies lower reliability. It shows that the reliability is affected by both the number of clusters and PINC. In most cases, the distance between PICP and the diagonal is greater when PINC = 50%. This inspires us to examine the average distance between the PICP and the diagonal over the entire PINC range, i.e., ACE, which measures the overall reliability by taking the mean. Figure 5(c) shows the ACE as a function of *K*. As can be seen, the ACE drops from around 15% to less than 1% when *K* increases from 2 to 9, and rebounds to 4% when *K* increases to 20. The considerable drop of ACE before *K* = 9 implies improved reliability due to clustering, which can be further attributed to the exclusion of sites that are weakly correlated. By contrast, the rise of ACE beyond *K* = 9 indicates less reliable scenarios with further refined clusters, probably due to the exclusion of those sites that present strong correlations.

Similar results are presented by interval scores. Figure 5(b) exhibits monotonically increasing interval scores as PINC increases. In addition, the interval scores are also largely affected by *K*. Similar to ACE, AIS gives a comprehensive evaluation to the interval scores of a given *K* by taking the mean over the entire PINC range. As presented in Fig. 5(d), the AIS reaches its peak at *K* = 6 and declines thereafter. Since higher interval scores indicate better performance in terms of sharpness, the results suggest *K* = 6 for the optimal sharpness.

As displayed in Fig. 5, both metrics suggest that a medium number of clusters results in the optimal performance of the generated scenarios. However, the two metrics do not agree on the optimal number of clusters: the highest reliability is present when *K* = 9 while *K* = 6 maximizes sharpness. Therefore, SEM is introduced and used as a comprehensive evaluation metric by factoring in both reliability and sharpness. Figure 6 shows the SEM values and the CRPS scores as a function of *K* in all four representative days. As can be seen, the winter, spring, and fall days exhibit very similar profiles, where the SEM is minimized when a medium number of clusters is selected. For example, the SEM is minimized when $K=$ 10, 12, and 10 in the winter, spring, and fall day, respectively. Although the summer day presents a slightly different profile, where the SEM is minimized at *K* = 15 and only rebounds slightly when *K* increases to 20, the minimum SEM still indicates the optimal overall performance of the generated scenarios when a medium number of clusters is selected. Interestingly, the CRPS scores present different profiles from the SEM: the scores remain low before *K* = 10 and increase when *K* > 10. Therefore, both the SEM and CRPS scores suggest against a high number of clusters, which can reduce the quality of the generated scenarios due to failure to account for correlations across wind farms. In addition, the CRPS scores also suggest benefits from clustering: the CRPS scores are at their lowest when *K* ranges from 2 to 4 in all seasons except for winter.

The generated scenarios can be directly used as inputs to the stochastic DAUC models. However, the number of generated scenarios must be reduced to an adequate number before the model can be solved in a time efficient manner. We follow the procedures in Algorithm 2 to reduce the generated scenarios. Figure 7 shows the time consumption as a function of *K*. The total time consumption is composed of scenario generation time and scenario reduction time. As shown in all four representative days, the total time consumption sees a trade-off between scenario generation and scenario reduction. The time consumed by scenario generation contributes to over 90% of the total time consumption when *K* < 3, while scenario reduction consumes over 65% of the total time when *K* increases to 20. In addition, as *K* increases, the time consumption of scenario generation presents a decreasing trend, while scenario reduction presents an increasing trend. For example, the average time consumption of scenario generation across all four days decreases from 140 min to 25 min when *K* increases from 1 to 20. By contrast, the average time consumption of scenario reduction increases from less than 1 min to more than 50 min over the same range of *K*. It is not surprising to observe the trade-off, since as *K* increases, each cluster is more likely to see a smaller number of wind farms, which reduces the scenario generation time. However, the number of iterations in scenario reduction is directly related to the number of clusters, and an increased *K* leads to longer scenario reduction time.

The trade-off of time consumption between scenario generation and reduction again suggests that a medium number of clusters should be selected. In addition, a medium number of clusters can also significantly reduce the total time consumption. As shown in Fig. 7, the total time consumption tends to be minimized when $5\u2264K\u226410$.

In addition, Table II compares the baseline (*K* = 1) with the clustered cases (*K* > 1). The comparison indicates that the clustering-based framework can significantly improve the generated scenarios in terms of SEM and time consumption. As can be seen, the SEM in the clustered cases (*K* = 2–20) can be reduced by up to 33%, while the total time consumption can be reduced by over 60%. In addition, the CRPS score can also benefit marginally from clustering. Therefore, the clustering-based framework can significantly enhance the quality of the generated scenarios and reduce the total time consumption.

. | . | Winter . | Spring . | Summer . | Fall . |
---|---|---|---|---|---|

SEM | Optimal K | 10 | 12 | 16 | 10 |

Baseline | 32% | 32% | 30% | 23% | |

Clustered | 26% | 27% | 20% | 16% | |

Improvement | 20% | 16% | 33% | 28% | |

CRPS | Optimal K | 2 | 4 | 2 | 3 |

Baseline | 0.249 | 0.242 | 0.277 | 0.191 | |

Clustered | 0.259 | 0.231 | 0.268 | 0.186 | |

Improvement | −4% | 5% | 3% | 3% | |

Total time (min) | Optimal K | 8 | 6 | 8 | 7 |

Baseline | 123 | 125 | 134 | 178 | |

Clustered | 49 | 49 | 51 | 54 | |

Improvement | 61% | 61% | 62% | 70% |

. | . | Winter . | Spring . | Summer . | Fall . |
---|---|---|---|---|---|

SEM | Optimal K | 10 | 12 | 16 | 10 |

Baseline | 32% | 32% | 30% | 23% | |

Clustered | 26% | 27% | 20% | 16% | |

Improvement | 20% | 16% | 33% | 28% | |

CRPS | Optimal K | 2 | 4 | 2 | 3 |

Baseline | 0.249 | 0.242 | 0.277 | 0.191 | |

Clustered | 0.259 | 0.231 | 0.268 | 0.186 | |

Improvement | −4% | 5% | 3% | 3% | |

Total time (min) | Optimal K | 8 | 6 | 8 | 7 |

Baseline | 123 | 125 | 134 | 178 | |

Clustered | 49 | 49 | 51 | 54 | |

Improvement | 61% | 61% | 62% | 70% |

## VI. RESULTS: POWER MARKET SIMULATION

As mentioned in Sec. III, the DAUC model will be solved deterministically and stochastically. Each deterministic run solves only one wind scenario, while the stochastic run solves all scenarios simultaneously. The 1000 wind scenarios generated in Sec. V with the minimum SEM values are reduced to ten scenarios and are represented with “1” to “10.” Note that the stochastic run with scenarios generated from the clustered case is denoted by “SC” and we use “SU” to denote the stochastic run with scenarios generated from the baseline, i.e., the unclustered method. In the clustered cases, the numbers of clusters are selected based on the minimum SEM. In addition, another run (denoted by “A”) takes the actual wind power as inputs to the DAUC model as if we have perfect foresight into the future. In summary, the DAUC model is solved 13 times, including 11 deterministic solves and two stochastic solves. Next, the RTED model is solved by using the actual wind power and committing all thermal units, where the commitment statuses are given by the DAUC model. In this study, Pyomo Stochastic Programming (PySP),^{63} an open-source Python library, is employed to implement the stochastic programming models and Gurobi is used to solve them. On the same compute node as used in Sec. V, it takes 20–25 min to solve a deterministic run, around 8 h to solve a stochastic run, and less than 5 min to solve a RTED model.

Results of the power market simulations are presented in Fig. 8 and Table III. As shown in Fig. 8, scenario A is always the least expensive in all seasons due to perfect foresight of wind power. Meanwhile, the more expensive total costs in other scenarios can be largely attributed to uncertainties associated with the forecast. In addition, great variations of total costs are presented in the deterministic scenarios. For example, in the fall day, the total cost from scenario 2 is 66% more expensive than that of scenario A, while scenario 3 is only 7% more expensive. Overall, the increased costs in the deterministic runs range from 2% to 103% across all seasons and scenarios, implying great uncertainties when the day-ahead market is cleared based on the deterministic wind power forecast. In contrast, the total costs of scenario SC only increase by 4%, 2%, and 4% in the spring, summer, and fall days, respectively. Note that although scenario SC results in a 50% more expensive total cost in the winter day, it is still lower than most of the deterministic runs in the same day, which range from 49% to 103%. The comparison between scenarios SC and SU indicates that the total costs are increased when the clustering-based framework is applied. However, the increases are only marginal in most cases: the total costs increase by 1%–3% in seasons 2, 3, and 4.

Season . | 1 . | 2 . | 3 . | 4 . | 5 . | 6 . | 7 . | 8 . | 9 . | 10 . | A . | SC . | SU . |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|

Winter | 6321 | 4034 | 5824 | 3047 | 3146 | 4095 | 5474 | 6411 | 4471 | 5141 | 0 | 3134 | 0 |

Spring | 4432 | 1586 | 4571 | 2632 | 2062 | 3795 | 4716 | 3921 | 1556 | 4326 | 803 | 1253 | 803 |

Summer | 735 | 2935 | 2949 | 1336 | 3972 | 905 | 523 | 978 | 432 | 191 | 0 | 177 | 0 |

Fall | 2479 | 3583 | 365 | 807 | 3304 | 1498 | 2758 | 2124 | 2891 | 2138 | 0 | 223 | 0 |

Season . | 1 . | 2 . | 3 . | 4 . | 5 . | 6 . | 7 . | 8 . | 9 . | 10 . | A . | SC . | SU . |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|

Winter | 6321 | 4034 | 5824 | 3047 | 3146 | 4095 | 5474 | 6411 | 4471 | 5141 | 0 | 3134 | 0 |

Spring | 4432 | 1586 | 4571 | 2632 | 2062 | 3795 | 4716 | 3921 | 1556 | 4326 | 803 | 1253 | 803 |

Summer | 735 | 2935 | 2949 | 1336 | 3972 | 905 | 523 | 978 | 432 | 191 | 0 | 177 | 0 |

Fall | 2479 | 3583 | 365 | 807 | 3304 | 1498 | 2758 | 2124 | 2891 | 2138 | 0 | 223 | 0 |

A closer examination of the cost breakdown indicates that the differences in total costs are primarily driven by variations in reserve shortages. As shown in Fig. 8, production costs are almost equal across all scenarios in the same day, while costs of reserve shortages vary significantly. Note that regulating reserve requirements are well met in all scenarios and shortages of responsive reserves are shown in Table III. Similar to the results of total costs, scenario A always results in the least amount of shortages and scenario SC results in less shortages than most deterministic scenarios. Therefore, the stochastic program represents a better strategy for market scheduling than the deterministic approach in terms of reliability and economy.

## VII. CONCLUSION AND DISCUSSION

This paper presents an enhanced scenario generation framework by integrating existing approaches with cluster analysis. As demonstrated in the results, the generated scenarios can benefit significantly from the clustering based framework in terms of reliability, sharpness, and time consumption. While only one scenario generation approach is investigated in this study, the framework can also be extended to many other scenario generation techniques. Therefore, this framework can provide proper support to real-world power system analysis with large-scale wind penetration.

Several insights can be drawn from our study. First, as shown in Sec. V, the SEM of the generated scenarios is optimal when a medium number of clusters is used, due to trade-off between intra-cluster and inter-cluster correlations. In addition, the time consumption is also affected by the number of clusters: a small cluster with less members requires less time for scenario generation; however, the increased number of clusters leads to an increased amount of time for scenario generation. Although the CRPS scores suggest a smaller optimal number of clusters, they also indicate that the quality of the generated scenarios clustering can be improved by clustering. The results also imply that it may negatively affect the quality of generated scenarios when all wind farms are grouped into one single cluster regardless of the strength of correlations, which highlights the necessity of the clustering-based framework.

In addition, as indicated in the power market simulation, the hedging strategy given by the stochastic DAUC model presents great advantage over most deterministic runs in terms of both economy and power system reliability. The reserve shortages are effectively reduced without significantly increasing the total costs. Note that even though some deterministic runs may result in lower total costs than the hedging strategy, i.e., scenario 4 in the winter day, it is more likely to result in greater loss by adhering to the deterministic strategy, since all other deterministic runs lead to greater total costs and reserve shortages. By contrast, the stochastic strategy outperforms most deterministic runs in terms of both economy and power system reliability.

Several caveats exist. Figure 8 indicates reserve shortages are present across all scenarios, which can be attributed to multiple factors. First, we assume ancillary services are only provided by thermal units, while other types of qualified resources, such as hydro or even load resources, are excluded. In addition, the economic loss due to reserve shortages may be underestimated as a result of the fixed penalizing prices, since the dynamic scarcity pricing mechanism in ERCOT allows the prices to rise when a shortage of operating reserves is present.^{58} Finally, ERCOT may procure replacement ancillary service capacity from a supplemental ancillary services market (SASM),^{60} which is not co-optimized with the energy market. Therefore, the reserve shortages do not necessarily reflect resource insufficiency. However, although this study may benefit from including the above factors, the main emphasis is placed on demonstration of the clustering-based scenario generation framework, and modeling of the above factors is beyond our scope.

Besides, we should note that although the stochastic run provides a comparatively better strategy for the power market operation, the time consumption is less acceptable to the operators. In this study, it takes around 8 h to solve its extensive form using state-of-the-art commercial solvers, while it takes 1.5 h to execute day-ahead reliability unit commitment (RUC) in ERCOT.^{60} The slow performance is largely intrinsic to mixed integer programs (MIPs), which is proved to be NP-hard.^{64} Therefore, in the future, additional techniques such as in Refs. 65 and 66 should be employed for a better and more stable time performance.

Another caveat, as shown in Fig. 8, is the more expensive total costs when the clustering-based framework is used to generate scenarios for stochastic runs. However, this does not indicate that the clustering-based framework is outperformed by the baseline in terms of solution quality, since only one number of clusters (*K*) are examined in each season due to the extremely long computational time of stochastic runs. Ideally, the solution quality should be examined over a range of cluster numbers. In addition, the optimal *K* is selected when the SEM is minimized and an underlying hypothesis is that a better SEM of probabilistic scenarios is related to a better solution quality of stochastic optimization. The fact that SEM is weakly correlated with the solution quality implies that further examination of other evaluation metrics and comprehensive studies should be conducted to examine the relationship between the evaluation metrics of probabilistic forecasting and the solution quality of stochastic optimization. Moreover, a real-world decision-maker may need a unique and comprehensive evaluation metric to determine the best set of scenarios, which should account for a variety of factors, such as scenario quality (reliability and sharpness), computational time, and solution quality. In fact, there has been little research into the trade-offs between the quality of stochastic optimization solutions and the quality of probabilistic scenarios and such studies are critical given the discrepancies presented in our study.^{67}

Finally, the clustering strategy in Sec. II is based purely on the cross-correlations across wind sites, which only captures the spatial correlation and fails to capture time dependencies. While temporal correlations are accounted for in the scenario generation in the form of autocorrelation, the quality of the generated scenarios could benefit from taking into account more complex dependence structures. Particularly, the spatiotemporal dependencies across wind sites can be complex, and a localized understanding of this structure may provide more targeted information to improve forecasts. For example, previous research suggests that the characterization of spatial correlation across wind sites can be improved by identifying the wind regime that is dominated by a wind direction since sites that are located along the wind direction are more likely to be correlated.^{68} Future studies will explore more clustering methods and criteria, such as distance-based methods, as studies have suggested superior performance when advanced techniques are adopted.^{69,70}

## ACKNOWLEDGMENTS

This work was supported by the National Renewable Energy Laboratory under Subcontract No. XAT-8–82151-01 (under the U.S. Department of Energy Prime Contract No. DE-AC36–08GO28308). This work was authored in part by Alliance for Sustainable Energy, LLC, the manager and operator of the National Renewable Energy Laboratory for the U.S. Department of Energy (DOE) under Contract No. DE-AC36- 08GO28308. Funding provided by U.S. Department of Energy Office of Energy Efficiency and Renewable Energy Wind Water Power Technologies Office. The views expressed in the article do not necessarily represent the views of the DOE or the U.S. Government. The U.S. Government retains and the publisher, by accepting the article for publication, acknowledges that the U.S. Government retains a non-exclusive, paid-up, irrevocable, worldwide license to publish or reproduce the published form of this work, or allow others to do so, for U.S. Government purposes.

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

## NOMENCLATURE

### Sets

### Parameters

*C*_{g}Slope of the cost function of generator $g\u2208GT$, in $/MWh

- $C\Delta L$
System-wide load curtailment cost, $/MWh

- $C\Delta R$
System-wide regulating up/down reserve curtailment cost, $/MWh

- $C\Delta S$
System-wide spinning reserve curtailment cost, $/MWh

- $CgSU,CgSD$
Start-up and shut-down cost of generator $g\u2208GT$

- $Cg0$
Intercept of the cost function of generator $g\u2208GT$, in $

- $Db,t$
Demand at bus $b\u2208BL$ in time $t\u2208T$

*K*Number of clusters during scenario generation

- $pg,0$
Initial power production of generator

*g*- $Pg,tF$
Forecasted power at time

*t*from non-thermal generator $g\u2208G\u2216GT$- $Pgmax,Pgmin$
Upper and lower bounds of power production from thermal generator $g\u2208GT$

- $RgSU,RgSD$
Start-up/shut-down ramp rates of generator

*g*- $RgU,RgD$
Ramp up/down rates of generator

*g*- $RQtRD$
System regulation down reserve requirement at time

*t*- $RQtRU$
System regulation up reserve requirement at time

*t*- $RQtS$
System spinning reserve requirement at time

*t*- $SF\u2113,b$
Power shift factor, power flow increment on branch $\u2113$ due to power injected at bus

*b*- $TgU,TgD$
The minimum online/offline time of generator

*g*- $TgU0,TgD0$
Number of initial time slices generator

*g*must be online/offline- $ui,ta$
Marginal CDF of $xi,ta$

- $ui,tf$
Marginal CDF of $xi,tf$

- $vg,0$
Initial commitment status of generator

*g*- $\Lambda \u2113$
Power limit of branch $\u2113\u2208L$

### Variables

- $pg,t\u2208\mathbb{R}+$
Power production from generator $g\u2208G$ in time $t\u2208T$

- $p\xafg,t\u2208\mathbb{R}+$
Maximum power available from generator $g\u2208G$ in time $t\u2208T$

- $rg,tRD\u2208\mathbb{R}+$
Regulation down reserve provided by generator $g\u2208GT$ in time

*t*- $rg,tRU\u2208\mathbb{R}+$
Regulation up reserve provided by generator $g\u2208GT$ in time

*t*- $rg,tS\u2208\mathbb{R}+$
Spinning reserve provided by generator $g\u2208GT$ in time

*t*- $vg,t$
Commitment status of generator $g\u2208GT$ in time $t\u2208T$, binary

- $xi,ta,xi,tf$
Samples of random variables $Xi,ta$ and $Xi,tf$

- $xi,ta,\xi $
Realization of variables $Xi,ta$ in scenario

*ξ*- $Xi,ta,Xi,tf$
Random variables indicating actual and forecasted wind power at site

*i*in time*t*- $zg,tP\u2208\mathbb{R}+$
Production cost of generator $g\u2208GT$ in time $t\u2208T$

- $zg,tSD\u2208\mathbb{R}+$
Shut-down cost of generator $g\u2208GT$ in time $t\u2208T$

- $zg,tSU\u2208\mathbb{R}+$
Start-up cost of generator $g\u2208GT$ in time $t\u2208T$

- $zt\Delta \u2208\mathbb{R}+$
Load and reserve curtailment cost in time $t\u2208T$

- $\delta b,t\u2208\mathbb{R}+$
Load curtailment at bus

*b*in time*t*- $\delta tRD\u2208\mathbb{R}+$
System-wide regulating down reserve curtailment in time

*t*- $\delta tRU\u2208\mathbb{R}+$
System-wide regulating up reserve curtailment in time

*t*- $\delta tS\u2208\mathbb{R}+$
System-wide spinning reserve curtailment in time

*t*

### APPENDIX: DERIVATION OF EQUATIONS

Copulas are employed to model the multivariate distributions in the denominator and numerator of Eq. (6). According to Sklar's theorem,^{71} any multivariate joint distribution can be written in terms of the marginal distributions of each component and a copula that describes the dependence structure between the components. Therefore, the joint CDF in Eq. (4a) becomes

where the copula $C2N:[0,1]2N\u2192[0,1]$ is a continuous, real-valued function and $uj,ta,uj,tf\u2208[0,1]$ are marginal CDFs associated with $Xj,ta$ and $Xj,tf$,

Similarly, the joint PDF in Eq. (4b) can be expressed as

where $c2N:[0,1]2N\u2192\mathbb{R}+$ is the 2*N*th order partial derivative of $C2N$,

and $faj$ and $ffj$ are marginal PDFs associated with $Xj,ta$ and $Xj,tf$,

Using the above copula representation, the numerator in (6) can be expressed as follows using the copula function,

Similarly, the denominator can be written as the following equation: