Risk-adjustable stochastic schedule based on Sobol augmented Latin hypercube sampling considering correlation of wind power uncertainties

The risks associated with wind power forecast (WF) deviations are of paramount importance to many power system participants (PSPs). However, traditional sampling approaches are computationally prohibitive to model these deviations. Additionally, setting a risk level for satisfying different PSPs receives little attention. This paper constructs a risk-adjustable stochastic day-ahead scheduling (RSDS) model to balance the risk requirements of PSPs, and proposes a Sobol-augmented Latin Hypercube Sampling (SaLHS) approach to improve sampling efﬁciency for scenario generation process in RSDS. At ﬁrst, SaLHS and D-vine copula are combined to generate WF error scenarios for RSDS considering correlations of wind farms. Speciﬁcally, SaLHS improves the uniformity and removes the correlation of random samples. Then, a Glue-VaR-based generation adequacy index (GVGAI) is proposed to measure operational risk. By adjusting the parameters of GVGAI, a desirable risk level can be obtained considering requirements of different PSPs. Furthermore, a multi-objective RSDS model is constructed considering operational cost and GVGAI. At last, an entropy-Weighted Aggregated Sum Product Assessment method is proposed to ﬁnd the best compromise solution for RSDS model based on the Pareto front obtained by an 𝜖 -constraint method. A modiﬁed IEEE-RTS system is used to validate the effectiveness of proposed method via numerical simulations.


INTRODUCTION
Day-ahead unit commitment (UC) is one of the most important applications in power system scheduling and operation. The solution methods for UC problems are categorized into deterministic and stochastic approaches [1], [2]. However, as the wind penetration rate increases, deterministic approaches are less desirable as it neglects the significant variability of wind energy. Stochastic UC (SUC) has been introduced as a promising tool to deal with uncertainties. The SUC considers scenariobased uncertainties to find the optimal solution. Therefore, a large number of stochastic scenarios should be generated to represent the randomness of wind power outputs. have correlations on each other. To generate wind power forecast error scenarios for meteorological correlated wind farms, generating multidimensional samples can be a prerequisite. Simple random sampling (SRS) is the most commonly employed sampling technique, but it faces immense computational challenges. Improvements have been made by stratified techniques such as Latin Hypercube Sampling (LHS), whereas LHS cannot provide good uniformity properties in a n-dimensional unit hypercube. Xie et al. [3] applied Quasi-Monte Carlo (QMC) methods in the sampling procedure of probabilistic optimal power flow because of their better multidimensional uniformity property. However, it has been found that QMC shows spurious correlations when applied to high-dimensional problems. When dealing with multivariate problems, sampling accuracy is affected not only by sample values, but also correlations among samples. Wang et al. [4] coupled LHS with Hammersley sequence to achieve uniformity in one and multi-dimensions. However, this method has not been evaluated for large dimensions. Nishant and Urmila [5] proposed a new sampling technique based on LHS and Sobol sequence (SS) to circumvent the issues of spurious correlations, and maintained the multidimensional uniformity. However, the permutation method used in SS cannot minimize the undesired correlations between samples. This motivates us to find a more reliable and efficient sampling method that maintains the multidimensional uniformity as well as minimizes the undesired correlations among samples.
System operators typically have to deal with two conflicting demands: management of operational costs and risks. A thorough SUC should consider evaluating benefits and risks of wind power integration and thus supporting trade-off decisions. Finding a trade-off between these two demands is a challenging task for managers to face on a daily basis. Besides, they need to quantify risks.
For SUC problems, the objective function is generally formulated as the minimization of expected cost. Some researchers adopted a stochastic security-constrained UC (SCUC) method to manage operational risks [6]- [8]. Risk measures for SCUC generally adopt probabilistic indices, for example, loss of load probability (LOLP) or expected load not served (ELNS), which cannot tell the risk level for a specific decision. The meanvariance (MV) model not only pays attention to the economic cost but also attaches great importance to the economic risk. Li et al. [2] constructed a risk constrained MV for the stochastic economic dispatch problem. This method is based on the variance of fuel cost, which can be further improved by incorporating risk measures. Some researchers adopted VaR and CVaR as risk measures in UC [9]. However, what constitutes a suitable risk measure requires further investigation.
From the system operators' perspective, controlling the operational risk is fundamental to protecting the benefit of customers and utilities, which may have conflicting objectives. Besides, different participants of the power industry, such as independent system operators (ISO) and system regulators, could have different attitudes towards risk. However, to the best of authors' knowledge, little attention has been paid on risk measurement considering different participants. To this end, this paper introduces a GlueVaR based generation adequacy index (GVGAI) to measure the risk. GlueVaR (GVaR) [10] is a flexible risk measure, which provides a risk assessment that lies between VaR and CVaR under different confidence levels.
In the light of previous discussion, a multi-objective (MO) stochastic day-ahead scheduling (SDS) model considering operational cost and risk can be formulated based on SUC. To solve the MO model, most literature adopted weighted sum method (WSM). However, this approach fails to coordinate conflicting objectives. Another method, Pareto optimization aims to search a Pareto solution set, and one solution is extracted as final decision [11]. Inspired by decision-making science, an entropy-Weighted Aggregated Sum Product Assessment (WAS-PAS) method is proposed to acquire the best compromise solution. Based on the above analysis, there remain research gaps which can be summarized by the following aspects: 1. Scenario generating methods for meteorologically correlated wind farms. A considerable part of works in this field have studied sampling techniques of wind power forecast error scenarios. However, there remain challenges because current sampling methods still face immense computational burdens. The sampling efficiency can be improved by improving multidimensional uniformity and minimizing the undesired correlations among samples. 2. Suitable risk measure for different power system participants. Traditional risk measures may limit the efficiency of the generation schedules, and may even lead to infeasible solutions in some cases. In addition, these measures fail to consider the risk preference of decision makers. Therefore, it is necessary to develop a flexible approach to improve the overall efficiency of scheduling, which could not only enable flexibility to adjust the benefits and risks, but also consider the SO's risk preference by adjusting the risk level.
The main contributions of this paper are listed as follows: 1. A stochastic scenarios generation method using a Sobolaugmented Latin Hypercube Sampling (SaLHS) and D-vine Copula is proposed considering correlation of wind speed forecast error. 2. A mean-risk (M-R) risk-adjustable SDS (RSDS) is proposed based on GVaR, which provides a flexible tool to control operational risks for different participants. 3. An entropy-WASPAS method is proposed to obtain the best compromise solution based on -constraint method to solve the RSDS model.
In order to illustrate the proposed procedure, a flowchart of the procedure is illustrated in Figure 1.
The rest of this paper is organized as follows: Section 2 describes scenario generation method based on SaLHS and D-vine Copula. Section 3 presents theoretical foundation of GVaR and proposes GVGAI. Section 4 constructs the RSDS model and explains the entropy-WASPAS to solve the RSDS model. Section 5 presents case studies of the proposed approach. Conclusions of this work are presented in Section 6.

2.1
Scenario generation based on D-vine copula

Copula and D-vine copula
In the case of multiple wind farms, wind speed forecast errors of different wind farms are dependent. The most straightforward way to address correlations of multiple wind farms' forecast errors is to obtain the joint distribution. Copula theory provides a convenient way to model multivariate joint distribution. As stated by the Sklar theorem, any multivariate joint distribution can be written in terms of N univariate marginal distribution functions and a copula function [12]. Elliptical and Archimedean copulas are common options for copulas. Although both classes are implemented successfully in bivariate applications, their multivariate extensions suffer from significant deficiencies.
To cope with their drawbacks, vine copula was proposed to construct high-dimensional dependence structures by combining different bivariate copulas. Two special types of vine copula, namely C-vine and D-vine, have been used extensively in literatures [2]. C-vine is skilled in modeling correlations between the dominated variable and other variables, while D-vine specializes in capturing correlations among arbitrary multi-variables [2]. As wind farms for this paper are geographically similar, D-vine is selected to construct vine structures.
Let D ⊂ {1, … , N } and m, n ∈ {1, … , N } ⧵ D, c mn|D (⋅, ⋅|x D ) is the copula density associated with the conditional distribution of (X m , X n ) given X D = x D . F m|D (⋅|x D ) denoting the conditional distribution of variable X m given X D = x D . For D-vine copula, a random vector X = {X m }, m = 1, … , N with marginal distribution functions F 1 , … , F N , the joint density f of the random vector X can be written as: To intuitively illustrate the D-vine structure, decomposition structure of D-vine copula with five variables are shown in Figure 2. The nodes in Tree1 represent the marginal densities of wind speeds for five wind farms. The edges in all trees are (conditional) bivariate copula densities. Model selection and sampling from D-vine copulas The procedure for fitting a joint distribution function using the D-vine copula can be briefly summarized in four steps.
Step 1. Model marginal distributions. Kernel density estimation method with Gaussian kernel function is applied to estimate marginal distribution of wind forecast errors.
Step 2. Select an order of variables for D-vine copula structures. The structure is determined by maximizing the sum of absolute values of pairwise Kendall correlation coefficients [13].
Step 3. Choose a bivariate copula for each pair-copula. Select the sequential estimation approach with Akaike Information Criterion as criterion [14]. Commonly adopted bivariate copulas include Gaussian, t, Frank, Gumbel, and Clayton copulas [12].
Step 4. Estimate all copula parameters sequentially by maximum likelihood estimation of each bivariate copula.
As the analytical formula obtained by the D-vine copula is complex, it is difficult to achieve the distributions of aggregated wind power forecast errors by multiple integrals calculation. The distributions from the D-vine copula can also be obtained by inverse transformation sampling. For the gener- (2)

SaLHS technique
For scenario generation based on sampling from the D-vine copula, the first step is to generate i.i.d random vector which follows uniform distribution. However, it is almost impossible for the commonly applied random permutation method to remove the undesired correlations among samples of different random variables [15]. SRS is a classical method to generate random samples. However, this method has disadvantages of poor repeatability and low efficiency [16]. As the sample size required for Monte Carlo Simulation (MCS) is directly proportional to the sample variance, variance reduction techniques (VRT) are designed to reduce the estimate variance without increasing the sample size [4]. LHS is one of the most popular VRT, whose efficiency is proved compared to SRS.
LHS process can be shown as follows: Consider the range of input variable divided into N intervals of the equal length 1∕N . One point is selected at random from each interval forming a sequence of N points in 1-dimension unit hypercube The two sequences can be paired to populate a bi-dimensional space, and so on until an n-dimensional sequence is formed. This algorithm can be shown as follows. Let {R k }, k = 1, .., K be independent random permutations of {1, … , N }. Then, the input sample is generated based on the inverse transform method and given by where U k i is independent and randomly sampled points on [0,1] interval, F k i is the inverse transformation of input variable k, R k (i ) is the ith element of R k .
As shown above, conventional methods pair input variables randomly, which brings undeserved correlation to generated samples.
In addition, stratification scheme of LHS cannot provide good uniformity properties in a n-dimensional unit hypercube H n . Sobol sequence has been proven to have better multidimensional uniformity. Nishant et al. [5] has proposed LHS-SOBOL method to pair samples generated from LHS using the ranks of N × k Sobol matrix M s . However, the sample correlation matrix of input variables generated by LHS-SOBOL is not exactly equal to identity matrix I, which reduces the efficiency and accuracy of sampling.
To cope with above problems, a pairing process is proposed in this paper to make the sample correlation matrix equal to I, and also maintain the d-dimensional uniformity of samples by taking advantage of the properties of Sobol matrix. The paring process of SaLHS is shown as follows.
Let R stand for the sample correlation matrix associated with Sobol matrix M s . The generation of Sobol matrix can be referred to [5]. R is used to find a matrix S so that where C is the is the desired sample correlation matrix. Let P be a matrix such that PP T = C. The Cholesky factorization can be used to find a lower triangular matrix Q such that Equation (4) can be rewritten as Therefore, S can be given by S = PQ −1 . According to [17], the matrix M * s = RS T has a correlation matrix exactly equal to C. The LHS samples can be paired by M * s rather than M s . Based on the pairing process, Sobol matrix and LHS, this paper proposes SaLHS. The procedure for SaLHS can be summarized as follows: where N sample stands for number of sample size, N w f stands for the number of wind farms. Each column comply with uniform distribution of [0,1].
Step 3. Generate permutation matrix M * s by paring process.
Step 4. Pair the Z LHS with M * s and generate sampling matrix U SaLHS .
SaLHS proposed in this paper generates samples with multi-dimensional uniformity and removes sample correlations. Numerical studies shown in Section 5 present its ability in improving sampling efficiency.

Scenario generation from SaLHS-based D-vine copulas
The LHS samples paired by M * s can be used as the input variables for D-vine construction. Based on SaLHS and D-vine copula, scenarios for wind power forecast errors considering correlations of wind farms can be generated. The steps are as follows: Step 1. Obtain wind speed forecast error N WSE × N w f matrix E WS , and then fit a joint distribution function using the model selecting method introduced in Section 2.

GLUE VaR-BASED GENERATION ADEQUACY INDEX
The most widely used risk measurements for SUC are VaR and CVaR. However, VaR fails to consider tail loss, and lacks subadditivity and convexity. On the contrary, CVaR is able to quantify the risk potential beyond VaR. In addition, CVaR is a coherent risk aversion measure [18]. However, so far, little risk measures consider the risk preference of different participants.
This paper proposed a GVGAI. GVaR can be defined based on a distortion function [10] and acts as a bridge between VaR and CVaR. GVaR can be expressed as linear combinations of standard risk measures. Each combination of risk scenarios reflects a specific risk attitude of decision maker. Therefore, GVaR is a flexible and simple risk measure.
be a function such that g(0) = 0, g(1) = 1, and g is nondecreasing. Then g is called a distortion function.

Definition 2. Let g be a distortion function. Consider a random variable X and its survival function S X
Any GVaR risk measure can be described by means of its distortion function. Given two confidence levels and , the distortion function for GVaR is where , ∈ [0, 1] and ≤ , h 1 ∈ [0, 1] and h 2 ∈ [h 1 , 1]. Parameters h 1 and h 2 are called the heights of distortion function.
If the following notation is used, Then, GVaR can be expressed as a linear combination of three risk measures: CVaR at confidence level and and VaR at confidence level .
In addition, when dealing with fat tail losses (i.e. lowfrequency and large-loss events), risk managers are especially interested in the tail region. GVaR has been proven to be subadditive in the tail region, which is a desirable property for risk measures. The property of GVaR inspires us to measure the risk of generating capacity shortfall over a certain threshold.
The power that load demand excesses available capacity is defined as the shortage of available capacity (SAC). That is, where P SAC t is the SAC for time t . SAC is a shortfall between required operational capacity and real operating capacity of the system. Based on the definition of GVaR and SAC, GVGAI is defined as follows: GVGAI can be seen as the weighted sum of VaR and CVaR for different participants under threshold and , respectively. In addition, as GVGAI is a coherent risk measure in the tail region, the convexity of GVGAI is kept in the tail region.
Weights and the parameters of GVGAI can be adjusted according to the risk preference of different power system participants. Therefore, weights and the confidence level can be seen as controlling parameters in GVGAI.
For the weights in GVGAI, as the weights are determined by the distortion function, one possible way to select weights is to determine the parameters of distortion function, that is, h 1 and h 2 . Once the confidence levels are determined, h 1 and h 2 can be used to determine the weights. The determination of h 1 and h 2 indicates how the decision makers distort the survival function. The determination of h 1 and h 2 could amplify or diminish the value of survival function, therefore increasing or decreasing the risk value. h 1 and h 2 control the shape of distortion function. Higher h 1 and h 2 means the higher gradient in the first part and second part of distortion function, respectively. The other way is to directly determine the weights 1 , 2 and 3 . Higher 1 indicates the decision maker puts more emphasis on risk level VaR . Similarly, higher 2 and 3 indicates the decision maker puts more emphasis on CVaR and CVaR , respectively.
For the confidence level in GVGAI, 95% and 99% are selected because these two confidence level are commonly adopted and compared in literatures [21]. The decision makers can select proper confidence level according to their risk preference.

Objective function formulation
For SUC, the objective function is generally formulated as minimization of expected operational cost. However, with increasing penetration of wind power, the uncertainty of wind farms' power outputs has a great impact on the economical and reliable operation. Considering the risk brought by wind power uncertainty to generation capacity adequacy of power systems, we propose a M-R MO RSDS. The risk is measured by the GVGAI. By adjusting the weights, GVGAI is able to reflect the risk attitudes of decision makers. The objectives of RSDS model is formulated as follows: where C u i,t stands for the start-up/shutdown cost of generator i at time t , s is the probability of scenario s, C p i,t ,s is the fuel cost of generators, RC t ,s is the reliability cost, which consists of loss of load and wind curtailment cost, and are confidence levels which can be specified by decision makers.

Constraints formulation
Given the previous modeling assumptions, the constraints for RSDS model are presented as follows: Constraint (17) is the demand-supply balance constraint. Constraints (18) and (19)

Solution algorithm based on entropy-WASPAS
MO problems are difficult to solve because there is a set of acceptable solutions within a range of values, which is called a Pareto front. Traditionally, the WSM is adopted to convert MO problem into a single-objective problem. However, the solution of WSM strongly depends on the weight coefficients. The weight coefficients of WSM are easily influenced by the subjectivity. In addition, WSM is not necessarily equivalent to the original MO problem.
-constraint method is efficient in obtaining Pareto front [19]. The first step of -constraint method is to establish the payoff table. An example of payoff table is shown in Table 1. Take objective function f i as the objective to solve RSDS model. At this time, the optimal value of f i is f i,i and the other objective functions under this solution are f i, j .
Next, consider one of the objective, such as f 1 , as the main function. Range of f 2 values is divided into q ranges, varying the value of k 2 . The -constraint method is formulated as: where k 2 is the kth range of f 2 , r 2 and q 2 are the range of f 2 , and the number of the equal parts, respectively. Therefore, the Pareto solution ps i p can be obtained.
The selection of the best compromise solution is done by evaluating the relative significance for each solution. WASPAS technique, suggested initially by Zavadskas et al. [20], is considered as a more accurate approach than WSM. However, the weight coefficients of WASPAS are selected by subjective empowerment methods. To avoid subjectivity, this paper introduces entropy weighting method into WASPAS and formulate a entropy-WASPAS method to evaluate the relative significance for each solution so as to select the best compromise solution. The entropy weight method uses the information entropy to calculate the entropy weight of each index according to the variation degree of the objective function. The steps of entropy-WASPAS technique are illustrated as follows: Step 1. Calculate the entropy E i of the objective function f j : Step 2. Calculate the weight w j of the objective function f j : Step 3. Normalize objective functions of the minimum operational cost and operation risk by Step 4. Determine the total relative importance of the ps i according to WSM: Step 5. Determine the total relative importance of the ps i according to WPM: Step 6. To determine the relative significance, apply a joint weighted aggregation criterion: i . The weight generally equals to 0.5. Step 7. Select the optimal compromise solution by the rank of Q i .

Verification of SaLHS-based D-vine copula scenario generation model
The hourly wind speed data from 8 geographically closely related wind farms were collected. The data spans 6 years from 1st January 2007 to 31st December 2012. Denote marginal distribution of wind speed forecast errors by u1∼u8. The frequency histograms for bivariate joint distribution of u1∼u8 are shown in Figure 3, which are ordered by the method illustrated in Section 2.1.2. As shown in Figure 3, the histograms show asymmetric tail distribution. In addition, there is no general method to describe the correlation of all variables. Therefore, it is crucial to consider the tail information so as to capture the dependence structure.
To improve the sampling efficiency for D-vine copula, SaLHS is applied in the sampling procedure.
SaLHS reduces the computational burden, and keeps the accuracy of sampling results. The superiority of SaLHS results Sampling results of 6 different sampling techniques, that is, MCS, LHS, Halton Sequence, Sobol Sequence, LHS-SOBOL [5], and SaLHS, are compared. Figure 4 compares 1dimensional uniformity of different sampling techniques. Better 1-dimensional uniformity is demonstrated by the proximity to the 45 • line with a uniform distance among the neighboring sample points. As shown in Figure 4, SaLHS shows better one-dimensional uniformity. Figure 5 shows generated 2-dimensional samples. It is clear that Halton and Sobol have better uniformity than MCS and LHS for multi-dimensional approximation. However, with the dimension increasing, these methods show spurious correlations. SaLHS is shown to break these correlations. This can be quantified by Variance Inflation Factor (VIF), which is calculated to identify the pairing correlations [4].
The VIFs of each sampling technique are compared in Figure 6. The VIFs of dimensions 1, 2, 74, and 75 of each sampling method are intuitively shown in Table 2. As shown in Figure 6 and Table 2, the VIF of SaLHS is the lowest, which shows its superiority in overcoming undesired correlations.
In order to compare sampling methods, accuracy of the sampling results is estimated by the relative errors of the expectation = | a − s a |, where a and s stand for the expectation of actual and sampling results, respectively. Three commonly applied techniques are compared based on̄V , which is the mean from 8 wind farms. Results are shown in Figure 7. As shown in Figure 7, SaLHS exhibits obvious superiority with respect tōV , whereas MCS performs worst among sampling techniques. Although LHS performs better than SaLHS when sample size is less than 200, the sample size could  be over 1500 for LHS to get a stable result. On the contrary, when the sample size grows over 200, SaLHS performs obviously better than LHS. For SaLHS,̄V becomes stable when the sample size is equal or over 500. The simulation time of the three sampling methods is shown to further compare the computation cost. The terminal criterion is set as̄V ≤ 0.005. The simulation time of the three sampling techniques is shown in Table 3. As shown in Table 3, the SaLHS has the lowest computational cost, while the computation cost of MCS is the highest among the three sampling techniques.
In conclusion, comparative results verify SaLHS improves multi-dimensional uniformity and removes sample correlations. Therefore, SaLHS can obtain accurate sampling results with a smaller sample size, which helps reduce the computation burden. The comparison of simulation time proves the effectiveness of SaLHS over MCS and LHS. In the following study, the sample size is selected as 500 to generate wind power output scenarios.

Test systems
This paper selects the modified RTS system as the simulation system. The scheduling period is considered the first day of the year with 24 equal time slots, whose maximum and minimum loads are 2284. 73   The confidence level is assumed to be = 0.95, = 0.99. In addition, CPLEX is called by Yalmip to solve the problem based on Matlab coding performed on the aforementioned personal laptop.

Simulation and results
Based on the formulation of dependence among multiple wind farms, 500 wind power output samples for each wind farm at each hour are generated. To improve the efficiency of computing, the fast forward selection method based on Kantorovich distance is applied in RSDS model to reduce 500 wind output scenarios into 10 scenarios for each hour. For the convenience of visualization, the reduced scenarios are illustrated in Figure 8. Based on the reduced scenarios, RSDS model can be constructed. To solve the RSDS model, the payoff table is firstly obtained, as shown in Table 4. The Pareto front of the model obtained by -constraint method is shown in Figure 9. To select the best compromise solution based on Pareto optimal set, the method illustrated in Section 4.3 is used. The optimal compromise solution is marked red in Figure 9. The optimal compromise solution is f 1 = 603485.66, f 2 = −928.52, and the corresponding relative significance Q is 0.5963.

Comparative studies
This paper proposes a GVGAI, which provides a flexible approach that enables decision makers to consider risk preference by adjusting the risk level and weights. To verify the effectiveness of the GVGAI in day-ahead scheduling, we have pro- Pareto front of SUC model vided comparative studies on operational cost, traditional risk measure CVaR at different confidence levels, and GVGAI. Table 5 compares the performance of RSDS with 3 cases in terms of the operational cost and risk.
Case 1: The objective is modeled as the minimization of cost. Case 2: The objective function is modeled as the minimization of cost and CVaR at confidence level of 95%. Case 3: The objective function is modeled as the minimization of cost and CVaR at confidence level of 99%.
Observed from Table 5, Case 1 has the lowest operational cost and highest GVGAI because of its ignorance of operational risk. In addition, the optimal compromise solution is also affected by the risk preference of decision makers. Compared with Case 3, Case 2 is more conservative in terms of the operational risk, which leads to a higher operational cost.
The adjustable weights of GVGAI ( 1 , 2 , and 3 ) represent the preference of system operators. Figure 10 illustrates the variation in operational cost and risk with regards to the different values of weights, separately.
As observed in Figure 10, the operational cost decreases with w 1 and w 2 increasing. That is, with the increase of w 3 , system is considered to have a higher operation risk level. Scheduling strategies are taken to reduce the cost associated with risk. Therefore, the operational cost decreases.
In conclusion, GVGAI is more flexible than traditional risk measures such as VaR and CVaR. This results from that weights  and confidence level of GVGAI can be tuned to consider risk preference in the decision-making process. Comparative studies provide an illustrative example of how the scheduling result changes with the variation of GVGAI parameters. Results can provide support for decision makers of power systems.

Numerical simulation with IEEE 118-bus system
In order to illustrate the scalability and applicability of the proposed method with a large number of buses and generators, a numerical study is implemented on a modified IEEE 118-bus test system. Eight wind farms with the same capacity of 55 MW are connected to buses 1, 6, 32, 36, 62, 70, 73, and 104. Table 6 compares the operational cost, risk, and computational time of IEEE 118-bus system with three sampling methods. As shown in Table 6, the scheduling result of the proposed method is as accurate as traditional sampling methods. The superiority of SaLHS lies in its high efficiency in sampling. As observed in Table 6, the model solution time is not sensitive to the sampling techniques, but the total computational time of the proposed method is the lowest. Therefore, the proposed method reaches the same accuracy as traditional sampling techniques with a less computational burden. The proposed method can be effectively scaled up to larger and more complex power systems.

CONCLUSION
This paper focuses on SDS problem considering correlation of wind farms. The contribution of this paper is three-fold. First, to improve scenario generation efficiency, a SaLHS sampling method is proposed. SaLHS is introduced in D-vine copula to generate scenarios considering correlations of wind farms. Second, to improve decision-making flexibility, GVGAI is proposed. This makes operational risk adjustable according to risk preference of power system participants. Third, an entropy-WASPAS method is proposed to avoid subjectivity in compromise solution selection. Numerical results show that SaLHS has superiority in improving sampling efficiency. RSDS is a flexible tool for system operators to balance the operational cost and risk. Results compared the influence of decision makers' risk preference on the operational cost and risk. However, there still exist some limitations that require further discussion. The first one is the proposed method requires generating a lot of stochastic scenarios. The scenario generation and reduction are commonly adopted in stochastic optimization methods. This paper focuses on the scenario generation methods to improve sampling efficiency. Scenario reduction method is still necessary and important in reducing the computational burden. Therefore, the scenario reduction method can be further improved for computation efficiency. The second is the proposed method is limited to generation adequacy indices.
How to extend the proposed method in the reliability analysis and operational optimization in bulk power system is another interesting topic in future researches.
The proposed method gives a day-ahead scheduling strategy for system operators. Future work is planned for extending the application of RSDS model considering correlations of generators and transmission line failures. P i,t ,s Power output of conventional units u i,t , y i,t , z i,t on/off, startup, shutdown state of unit i during t P lc b l ,t ,s Load curtailment on b l for s during t P w i w ,t ,s Wind power output of i w for s during t P wc i w ,t ,s Wind power curtailment of i w for s during t P min i , P max i Minimum/maximum generation limit of i P L ib,tb Power flow limit of lines from ib to tb ib Phase angle of bus ib ORCID Jianwei Gao https://orcid.org/0000-0003-4203-0374