Capacity Credit Evaluation of Correlated Wind Resources Using Vine Copula and Improved Importance Sampling

This paper concentrates on the capacity credit (CC) evaluation of wind energy, where a new method for constructing the joint distribution of wind speed and load is proposed. The method is based on the skew-normal mixture model (SNMM) and D-vine copulas, which is used to model the marginal distribution and the correlation structure, respectively. Then a cross entropy based importance sampling (CE-IS) is improved to enhance the efficiency of the power system reliability assessment, which is a crucial part of the CC evaluation. After that, the proposed methods are adopted to combine with the secant method to develop a complete algorithm to calculate the CC of wind energy. Numerical tests are designed and carried out based on the IEEE-RTS 79 system and wind speed data obtained from four wind farms in Northwest China. In order to show the superiority of SNMM and D-vine copula, the goodness-of-fit is quantified by different statistics. Besides, the improved CE-IS method is validated by comparison with Monte Carlo sampling (MCS) and traditional CE-IS in the efficiency of reliability assessment. Finally, the proved methods are combined with the secant method to calculate the CC of four wind farms, which can provide information for wind farm planning.


Introduction
Clean energy resources such as renewable energy sources and flexible demand have been introduced into power systems in order to establish an environment-friendly society [1][2][3][4].Wind energy, a kind of clean energy, has experienced rapid development in recent years to aid lowering of greenhouse gas emission [5,6].However, wind farms cannot provide a steady power capacity, restricting the current power system from transitioning to one with high penetration of wind energy [7].The reason lies in the uncertain and intermittent nature of this resource.Therefore, in order to identify proper power planning schemes and reliably satisfy load demands, we need to estimate the contribution of uncertain wind energy to power systems.
Conventionally, the contribution of wind energy is quantified by CC, which is also called capacity value in some literature [8].CC can be defined as the additional amount of load the system can support after the integration of wind energy (or other power resource to be evaluated), with the reliability level unchanged [9].The calculation of CC is usually performed by the secant method, which is a relatively mature and steady algorithm [10,11].Hence, the precision of the evaluation of CC depends mainly on two factors, i.e., the accuracy of the system reliability assessment and the wind speed model of uncertainty.The wind speed model is used to provide large numbers of stochastic scenarios for the assessment.
Traditional ways of modeling the marginal probability density function (PDF) of wind speed can be categorized into sequential or non-sequential methods.
In non-sequential models, the most commonly used paradigm is the Weibull distribution, which can be given by two parameters.The utilization of Weibull distribution was also reported in [12,13], which aimed at estimating wind energy resources, wind power forecasting, and so on.However, individual distribution, such as the Weibull and Rayleigh-Rice distributions, can hardly be suitable for capturing the characteristics of wind speed under various topographical and climatic conditions.Therefore, the focus of recent researchers has shifted to applying mixture distributions to fit the distribution of wind speed, among which the Gaussian mixture model (GMM) is typical.Ke, Chung and Sun in [14] proposed a customized GMM to deal with the discontinuities of the PDF, based on which a novel probabilistic optimal power flow model was established.GMM was also combined with neural network to forecast short-term wind power generation in [15].
The sequential models of wind speed are usually constructed by time series models.Morales, Mínguez and Conejo in [16] employed an auto-regressive moving average (ARMA) process to model wind speed characteristics, which involves transformation of the historical data into normalized Gaussian time series before fitting the model.ARMA models have also been used for wind speed modeling in reliability studies of electrical power systems [17].Besides, auto-regressive integrated moving average (ARIMA) models are also reported to perform well in forecasting wind speed precisely [18].The time series models have also been combined with empirical mode decomposition (EMD) to deal with the high nonlinearity and instability of wind speed series in some researches [19,20].
Besides methods to model the marginal PDFs of wind speed, efforts have been made to construct the spatial dependence between wind speed and load at different sites to form an integral joint PDF.
A prevailing methodology to cover the dependence is the copula theory.For example, Gaussian copula function was applied in [21] to build a correlation structure using wind speed data available at the database of the National Oceanic and Atmospheric Administration.Xie, Li and Li in [22] employed multivariate Archimedean copulas to model wind speed dependence, which performs better than Gaussian copula in modeling asymmetrical dependence structure [23].However, it was concluded in [24] that the performance of multivariate Archimedean copulas deteriorates drastically if the number of wind and load sites increases.The deficiency cannot be overcome by simply using a mixture of multiple Archimedean copula functions.
Hence, researches on the pair copula method, which specializes in high dimension dependence construction, have been developing rapidly.The pair copula method is capable of representing a complex multivariate correlation structure as the product of multiple bivariate conditional copula functions, each of which can be chosen from any copula function family individually.The method greatly improves the flexibility and accuracy of the multivariate dependence modeling [25].In recent years, there have been researchers introducing the pair copula into the optimal planning and reliability assessment of power systems, whose graphical structure is called C-vine structure [26,27].C-vine copulas are suitable for cases where there is a dominant variable among all the correlated random variables (RV).For instance, if the wind speed at one site is relatively highly related with the wind speed and load at all other sites, C-vine will exhibit satisfactory efficiency and precision.
Using the joint PDF of load and wind speed, we can evaluate the CC of wind energy based on the calculation of reliability indices.In the context of high reliability of power systems, various methods have been developed to accelerate the process of reliability assessment, among which the CE-IS method is a most popular one [28,29].By combining CE-IS and the secant method, we can obtain a complete algorithm for CC evaluation.
However, there are drawbacks in this algorithm.Firstly, GMM fails to consider the higher-order statistics of the marginal PDFs, which may result in an excessive number of components.Secondly, C-vine copulas need the existence of a dominant RV.Otherwise, its performance declines greatly.Thirdly, the traditional CE-IS method assigns one RV to every conventional unit or wind turbine generator (WTG) to represent its working state, while too many RVs hinder the performance of CE-IS.
To solve the first and second problems, we develop a method of constructing a joint PDF by combining SNMM with D-vine copulas.SNMM owns superiority over GMM by considering higher-order statistics, which enables it to capture features of the PDFs with fewer components [30].D-vine copulas are introduced to link the marginal PDFs fitted by SNMM to form an integral joint PDF.The D-vine structure needs no dominant RV, and can thus be applied to more general situations [31].
As for the drawback of CE-IS, we seek to improve it by offering a dimension reduction (DR) idea.The newly developed method is called DRCE-IS.The DR idea is implemented by using binomial and trinomial distributions to group two-state and three-state generators, which reduces discrete RVs.DRCE-IS is then designed to accommodate the two distributions to fulfil the speed-up of the power system reliability assessment.
On the basis of the SNMM, D-vine copulas and DRCE-IS, we can obtain reliability indices rapidly and precisely, offering necessary information for the secant method to further evaluate the CC of each wind farm.
The rest of the paper is organized as follows.Section 2 details the process of SNMM for marginal PDF construction.Section 3 describes the basic ideas of pair copulas and D-vine structure.In Section 4, the DRCE-IS method is presented with SNMM and D-vine copulas integrated; the pre-simulation and main simulation stages are detailed.Section 5 discusses the concept and procedures of the secant method.Numerical tests are performed and discussed in Section 6.The conclusions are summarized in Section 7.

Model Marginal Distribution of Load and Wind Speed Using SNMM
For an arbitrary PDF, there is no unique methodology to model it accurately enough, which propels the study on GMM to provide a rather satisfactory approximation.On this basis, SNMM was developed as an improvement of GMM by covering the skewness of the PDF, observed in the expression of SNMM PDF later.
Therefore, SNMM inherits almost all the merits of GMM.For example, any distribution can be represented by as a convex combination of several skew-normal (SN) distributions with respective means and variances.If infinite components are adopted, SNMM can inerrably depict distribution of any kind.Besides, due to the improvement, SNMM is capable of representing the asymmetry, heavy tail and multimodality of PDFs, which is superior to GMM [32].
Because a PDF must be nonnegative and the integral of a PDF over the sample space of RVs must evaluate to unity, the mixture weights must be nonnegative and the sum of all the weights must be equal to one.Hence, for the univariate case, the parametric form of the PDF of SNMM f (x lw ) is: where x lw represents the RV of load or wind speed.lw is the index for the sites of load and wind speed.
ω m is the weight of the m-th component of the SNMM.The sum over ω m is 1. n SN is the number of components in SNMM.f SN,m (•) is the PDF of a single SN distribution, whose expression is shown as: where ϕ is the PDF of standard Gaussian distribution.Φ is the corresponding cumulative distribution function (CDF) of ϕ. µ m , σ m and λ m are the parameters reflecting the location, scale and skewness of the m-th SN distribution, respectively.The estimation of all parameters is conducted by the expectation maximization (EM) algorithm, which is a very mature and steady method applied in various mixture models [33].
According to [34], the CDF of SNMM F(x lw ) can be expressed by the addition of the CDFs of the SN components, shown as: where the CDF of a single SN distribution F SN,m (x lw ) is: where T(•) is the Owen's T function.
Based on (3) and ( 4), we derived the theoretical expression of the CDF of SNMM.Therefore, when we draw samples from SNMM, we can directly use the inverse CDF (iCDF) of SNMM, although the iCDF may not be an explicit expression [35], i.e., x lw = F −1 (U), where U is uniformly distributed in [0, 1].

Link the Marginal PDFs of Wind Speed and Load by D-Vine Copulas
A large numbers of scenarios is needed to approximate the reliability indices of power systems.Historical data are finite, so many extreme scenarios, such as the case of high load and low wind speed, may not be observed, which by contrast contributes more to the reliability indices.
Consequently, we should try to unveil the underlying multivariate distribution as a continuous function from the historical data.We can then generate sufficient scenarios to compute the reliability indices with satisfying accuracy.In this paper, we resort to the copula theory.

Sklar's Theorem
According to Sklar's theorem, for any n-dimensional CDF F(x 1 , . . ., x n ), there exists a function which satisfies the equation [36]: where C(•) is called copula function and F i (x i ) is the marginal CDF of x i , i = 1, . . ., n.
The copula density function c(•) can be defined as: The copula density function c(•) owns a particular property, shown as: where f (•) and f i (•) are the PDFs corresponding to F(•) and F i (•), respectively.According to (7), one can model the marginal PDFs and choose a proper copula to form the joint PDF.The form of the copula function can be determined by the test of goodness fit [27].

Pair Copula Theory
In terms of selecting a copula function, a single multivariate copula may not be a good choice.Hence, it is necessary to find a more accurate and flexible method.
In this context, pair copula theory is proposed to decompose a multivariate copula into multiple conditional bivariate copulas [37].The concept is introduced as follows.
Firstly, the joint PDF f (x 1 , . . ., x n ) can be decomposed as: we then try to rewrite the conditional PDFs in (8) with conditional bivariate copula functions, which is shown as: where x v-o is the subset formed by excluding x o from the set x v , which is a subset of {x 1 , . . ., x n }. c io|v-o is the conditional bivariate copula of x i and x o under the condition of x v-o .The ( 9) is a recursive expression.Consider a 3-variable example to present the recursion as: where c 1,2|3 is the conditional copula function of x 1 and x 2 under the condition of x 3 .
We have already modeled the marginal PDF by SNMM, so the difficulty in (9) lies in how to obtain F(x i |x v-o ).F(x i |x v-o ) can be derived as [38]: where x v-o-e is the subset of x v with both x o and x b excluded.The proposal of ( 11) renders ( 9) practical, so we can use (9) to replace all the conditional PDF in (8).After the substitution, f (x 1 , . . ., x n ) can be expressed by the product of n marginal PDFs and n(n − 1)/2 conditional copula functions.

Graphical Structure of Decomposition
As shown in (9), x o is picked up from x v arbitrarily.When we decompose (9) recursively, picking up x o in different orders generates various conditional copulas, which determines the graphical structure of decomposition for (8).
Generally, there are two graphical structures frequently used-C-vine and D-vine.C-vine is more suitable when there is a dominant RV, while D-vine performs better when the RVs have arbitrary correlation structures.
There is no dominant RV in the correlation modeling of the wind speed and load, so we adopt D-vine copulas to capture the dependence features.We then decompose Equation (8) as: where D(x 1 , . . ., x n ) stands for the D-vine copula function.
Graphically, the D-vine structure adopted in Equation ( 12) can be illustrated by Figure 1, which is a 4-dimension case.
Graphical structure of four-dimensional D-vine copulas.

Parameter Estimation
We can choose a suitable copula type for each pair of F(xk|xk+1, …, xk+j−1) and F(xk+j|xk+1, …, xk+j−1) according to the property of their dependence, such as the tail behavior [39].
After the function form of each pair is determined, the parameter estimation process is carried out, in which we use the maximum log-likelihood estimation (MLE).
The MLE process proceeds row by row down the vine shown in Figure 1, because F(xk|xk+1, …, xk+j−1) and F(xk+j|xk+1, …, xk+j−1) are derived from the conditional copula CDF in the row above them by (11).The expression of MLE is: where Y is the number of observations used for estimation.Θ is the set of all θk,k+j, which is the parameter of ck,k+j|k+1,…,k+j−1.
After obtaining Θ by MLE, the expression of the joint PDF can be determined according to (12).

Generate Correlated Samples by D-Vine Copula
Consider x1, …, xn as the variables representing wind speed and load.In fact, when we use CE-IS to assess system reliability, only in the first iteration of the pre-simulation stage we may use Dvine copulas to draw samples.In the subsequent iterations, we draw samples from the suboptimal distribution derived by CE-IS.If there are enough historical observations, we can use the simple random sampling of the observations in the first iteration as well.
However, to ensure the integrity of the methodology, the procedures of sampling from D-vine copulas is briefed below, based on Quasi-Monte Carlo method [40].
Step 1: Generate n independent samples distributed uniformly on [0, 1], which are denoted by z1, …, zn. Step , which has been saved while the bivariate conditional copulas are constructed by (11).F1(x1) is the marginal CDF of x1 modeled by SNMM.
When we evaluate the CC of wind energy, we use the data of wind power rather than wind speed.Hence, the wind speed samples are further transformed into the power of the WTG by: ( ) ( ) ( ) where vci, vr and vco are the cut-in, rated and cut-out wind speed of the WTG.xi can stand for any element in {x1, …, xn} that represents wind speed.PR is the rating of the WTG.

Parameter Estimation
We can choose a suitable copula type for each pair of F(x k |x k+1 , . . ., x k+j−1 ) and F(x k+j |x k+1 , . . ., x k+j−1 ) according to the property of their dependence, such as the tail behavior [39].
After the function form of each pair is determined, the parameter estimation process is carried out, in which we use the maximum log-likelihood estimation (MLE).
The MLE process proceeds row by row down the vine shown in Figure 1, because F(x k |x k+1 , . . ., x k+j−1 ) and F(x k+j |x k+1 , . . ., x k+j−1 ) are derived from the conditional copula CDF in the row above them by (11).The expression of MLE is: where Y is the number of observations used for estimation.Θ is the set of all θ k,k+j , which is the parameter of c k,k+j|k+1, . . .,k+j−1 .
After obtaining Θ by MLE, the expression of the joint PDF can be determined according to (12).

Generate Correlated Samples by D-Vine Copula
Consider x 1 , . . ., x n as the variables representing wind speed and load.In fact, when we use CE-IS to assess system reliability, only in the first iteration of the pre-simulation stage we may use D-vine copulas to draw samples.In the subsequent iterations, we draw samples from the suboptimal distribution derived by CE-IS.If there are enough historical observations, we can use the simple random sampling of the observations in the first iteration as well.
However, to ensure the integrity of the methodology, the procedures of sampling from D-vine copulas is briefed below, based on Quasi-Monte Carlo method [40].
Step 2: Set . Convert z i into x i by the inverse function of F(x i |x 1 , . . ., x i−1 ), which has been saved while the bivariate conditional copulas are constructed by (11).F 1 (x 1 ) is the marginal CDF of x 1 modeled by SNMM.
When we evaluate the CC of wind energy, we use the data of wind power rather than wind speed.Hence, the wind speed samples are further transformed into the power of the WTG by: where v ci , v r and v co are the cut-in, rated and cut-out wind speed of the WTG.x i can stand for any element in {x 1 , . . ., x n } that represents wind speed.P R is the rating of the WTG.

DRCE-IS Combined with D-Vine Copula and SNMM to Assess the Power System Reliability
The essence of CE-IS is to search for the optimal distribution for every RV, in order to reduce the variance of the reliability indices estimated by MCS.Based on the traditional CE-IS method, we propose the DRCE-IS method, which specifically aims at finding a better way to deal with the large numbers of 2-and 3-state units in the reliability assessment.

Overview of CE-IS
Denote all the RVs in reliability assessment by a vector x ss , i.e., where x L , x CU and x WT are the RV vectors representing the states of transmission lines, conventional units and WTGs.x LW is the RV vector representing the values of the load and wind speed at different sites.
Then the reliability index γ can be determined by: where I(x ss ) is a function to indicate the system state, e.g., in the normal or failure states.f (x ss ;θ) is the original joint PDF of x ss with θ as the parameter vector.Equation ( 16) can be further transformed by the optimal joint PDF g(x ss ) used in importance sampling (IS), which is called IS PDF in this paper.The transformed ( 16) is shown as: where the rightmost term is an expectation in which g(x ss ) and I(x ss )Q(x ss ) act as the PDF and a function of the RVs.The theoretical expression of the optimal IS PDF is [28]: where g(x ss ) is obviously entangled with γ, which means that the explicit expression of g(x ss ) is unavailable.
To solve the problem, the concept of cross entropy is introduced into the process of IS to form the traditional CE-IS method.The basic idea is to select a function form in advance, regardless of the actual form of the optimal PDF.Denote the selected function form as f (•;ω), in which ω is the parameter vector unknown.By minimizing the Kullback-Leibler (KL) divergence between f (•;ω) and g(•), we can derive a satisfactory suboptimal PDF as the IS PDF.
The KL divergence K(g, f ) is expressed as: where the first integral term is independent of ω, so it can be regarded as a constant C when we search for the optimal ω.Then, we plug (18) into (19) and get the objective function used for optimizing ω, shown as: The integral term in (20) can be further converted as: where E ω (•) stands for the expectation of a function with f (•;ω) as the PDF.n s is the sample size for estimating the expectation.Q(x ss ) is called the likelihood ratio function (LRF), which can be divided into several parts, shown as: There is no constraint imposed on (21), so the best ω can be achieved when the derivative of the last term in ( 21) is equal to zero.

Pre-Simulation Stage of DRCE-IS
DRCE-IS consists of a pre-simulation stage and a main simulation stage.We detail the pre-simulation stage first.The purpose of the pre-simulation stage is to derive IS PDFs for continuous RVs and IS probability mass functions (PMF) for discrete RVs.
We use an iterative process to eliminate the randomness of the samples in (21).In each iteration, we draw n s new samples from f (•;ω) to solve (21) again to renew f (•;ω) for the next iteration, until the terminating condition is satisfied.
IS PMFs, particularly, will hold the same function form as the original PMFs of the discrete RVs, while IS PDFs will be assigned a new function form, totally different from the original ones.

IS PMFs of Transmission Lines
The PMFs of the transmission lines are shown as: where ω l is the optimal failure probabilities of transmission lines.x l is the RV representing the states of the l-th transmission line, whose potential values are 0 and 1, indicating failure state and operation state, respectively.n L is the number of transmission lines.When we update ω l , only the terms containing ω l in (21) affect the IS PMF.Therefore, we retain them in the objective function, while other terms are replaced by a constant C.Then, the objective function is shown as: In the t-th iteration, draw n s new samples and make the derivative of (24) equal to zero, i.e., 1 The solution of ( 25) is the optimal failure probabilities of transmission lines, which is shown as: Appl.Sci.2019, 9, 199 9 of 21 Equation ( 26) is valid for any l = 1, . . ., n L .The set {x l |l = 1, . . ., n L } forms the vector x L in (22).According to (21), LRF is defined as the ratio of the original PDF (or PMF) to the IS one.Hence, the LRF of transmission lines Q(x L ) is given as: where θ l is the original failure probability of the transmission line l.

IS PMFs of Homogeneous 2-and 3-State Units
The wind farms in this paper are all centralized and integrated with the transmission network, which means that the WTGs on the same wind farm are definitely connected to the same bus in the power system.Hence, the WTGs of the same type on the same wind farm make equal contributions to the power system.These WTGs should have the same IS PMF.Similarly, the 2-state conventional units, which are of the same type and connected to the same bus, should have the same IS PMF.These units/WTGs are called "homogeneous units" in this paper.
If we search for the IS PMFs of the homogeneous units individually, the IS PMFs may significantly differ from each other due to the randomness of the samples, which are also far away from the theoretical optimal PMF.Besides, the state of each unit is regarded as a RV, which results in a high-dimension problem and a bulky solution space, increasing the difficulty and time in searching for optimal PMFs.
To solve the problem, we aggregate the homogeneous units into a single group and represent the group by one RV.For the group of 2-state conventional units, the RV follows binomial distribution, while for the group of 3-state WTGs, the RV obeys trinomial distribution.
By this means, the dimension of the solution space is reduced, which improves the quality of the IS PMFs.
Denote the number of 2-state conventional units in the c-th group as n c .The failure probability is marked as ω f,c .Denote x f,c as the number of units in the failure state in the c-th group.The PMF of the c-th group p(x c ) is shown as: where When we update ω f,c , we only retain the terms containing ω f,c in the objective function, i.e., max Similarly with the processing of transmission lines, we can obtain the solution of ( 29), which is the updating formulae for ω f,c , shown as: Equation ( 30) is valid for any c = 1, . . ., n cg .n cg is the number of conventional unit groups.
The set {x c |c = 1, . . ., n cg } forms the vector x CU in (22).Then the LRF of conventional units Q(x CU ) can be expressed as: where θ f,c is the original failure probability of the c-th group of conventional units.
Denote the number of 3-state WTGs in the w-th group as n w .The failure and derated probabilities are marked as ω f,w and ω d,w , respectively.The PMF of the w-th group p(x w ) is shown as: where x w = [x f,w , x d,w , x o,w ]. x f,w , x d,w and x o,w are the number of WTGs in the failure, derated and operation states in the w-th group, respectively, which meet: Similarly with ( 24), when we update ω f,w and ω d,w , the objective function is shown as: Similarly with (26), the solutions of ( 34) are the updating formulae for ω f,w and ω d,w , which are shown as: Equations ( 35) and ( 36) is valid for any w = 1, . . ., n wg .n wg is the number of WTG groups.
The set {x w |w = 1, . . ., n wg } forms x WT in (22).Then the LRF of WTGs Q(x WT ) is expressed as: where θ f,w , θ d,w and θ o,w are the original state probabilities of the w-th WTG group.We only deal with 2-and 3-state units in this paper.However, the proposed method is not limited to 2-and 3-state units.The method can be expanded to the case of multi-state units by replacing p(x c ) and p(x w ) with the PMF of the multinomial distribution, so it is widely adaptive.

IS PDF of the Wind Speed and Load Based on D-Vine Copulas and SNMM
Load and wind speed are represented by continuous RVs, so we need to acquire the IS joint PDF, rather than the PMF.According to the basic idea of CE-IS, there is no restriction to f (•;ω) in (19).Actually whatever the form of f (•;ω) is, we can obtain a practical IS PDF by optimizing ω.
For simplicity, we select the multivariate Gaussian PDF as the IS joint PDF, in which all variables are supposed to be independent.The expectation µ LW and standard deviation σ LW are the parameters to be optimized, just like ω in (19).Therefore, when we update µ LW and σ LW , the objective function can be represented by: max where diag(•) is a transformation by which a vector is converted into a diagonal matrix.
Similarly, with the processing of the PMFs, ( 38) is solved, whose solutions are the updating formulae for µ LW and σ LW , shown as: where (•) •2 stands for the Hadamard square.
The LRF of wind speed and load is defined based on the D-vine copulas and SNMM, because we need to combine them to obtain the original joint PDF.The LRF is expressed as: where n lw is the number of RVs of wind speed and load.f lw (•) is the marginal PDF modeled by SNMM.

Procedures of Pre-Simulation Stage
The procedures of the pre-simulation stage are detailed as below.
Step 1: Input the original PMFs and joint PDF of RVs in the reliability assessment.Set t = 0. Specify the sample size n s in each iteration.
Step 3: Quantify the performance of each sample in the term of reliability.The frequently used quantifying method is the minimum load shedding (MLS) model, whose objective function is the quantifying function S(x ss,s ).
The MLS model is expressed as: min P sh ,P g ,δ S(x ss,s ) = where P g , P sh and δ are the control variable of the model.P g is the vector of generation power at all buses including wind power, whose upper and lower limits are P gmax and P gmin .P sh is the vector of load shedding amount at all buses.δ is the vector of the voltage phase angle of all buses.P d is the vector of the load demand at all buses.b is a diagonal matrix, whose diagonal elements are the reciprocal of the reactance of all transmission lines.A is the branch-node association matrix.P linemax is the vector of maximum transmission capacity of all lines.
Step 5: In the t-th iteration, set a threshold α t for {S 1 , S 2 , . . ., S ns } by: According to the value of α t , we can define the system state I(x ss,s ) as: Step 6: By ( 26), ( 30), ( 35), ( 36), ( 39) and ( 40), we can update the parameters of the IS PMFs and IS PDF with the samples generated in Step 2 and the system states obtained in Step 5.
Step 7: If α t = 0, end the pre-simulation stage and switch to the main simulation stage.Otherwise, set t = t + 1 and return to Step 2.

Main simulation Pre-simulation
Accurate enough?Y N Compute S(x ss,s ), s=1,…,n s Compute S(x ss,q ) and I(x ss,q ) by ( 38) and (40).Store the results

Figure 2.
Flowchart of the DRCE-IS method.

Main Simulation Stage of DRCE-IS
The main simulation stage is implemented by adding a scaling operation to the MCS, which is to use the LRF Q(xss,s) to multiply the I(xss,s).The procedures are detailed as follows.
Step 8: Set the counter q = 0.The PMFs and PDF in the last iteration of the pre-simulation stage is accepted here.
Step 9: q = q + 1. Generate a new sample xss,q with the IS PMFs and PDF.
Step 12: Derive the reliability indices, the loss of load probability (LOLP) and expected power not supplied (EPNS) according to the stored data by: Step 13: Calculate the relative error of LOLP and EPNS.Compare them with the relative error limit.If LOLP and EPNS are accurate enough, end the algorithm and output them.Otherwise, return to Step 9.
The main simulation process is also shown in Figure 2.

Main Simulation Stage of DRCE-IS
The main simulation stage is implemented by adding a scaling operation to the MCS, which is to use the LRF Q(x ss,s ) to multiply the I(x ss,s ).The procedures are detailed as follows.
Step 8: Set the counter q = 0.The PMFs and PDF in the last iteration of the pre-simulation stage is accepted here.
Step 9: q = q + 1. Generate a new sample x ss,q with the IS PMFs and PDF.
Step 12: Derive the reliability indices, the loss of load probability (LOLP) and expected power not supplied (EPNS) according to the stored data by: Step 13: Calculate the relative error of LOLP and EPNS.Compare them with the relative error limit.If LOLP and EPNS are accurate enough, end the algorithm and output them.Otherwise, return to Step 9.
The main simulation process is also shown in Figure 2.

Capacity Credit Assessment of Wind Energy by the Secant Method
The CC of wind energy is frequently defined as the equivalent load carrying capacity (ELCC) [41].The ELCC is the load increment that is needed to maintain the reliability level unchanged after wind energy integration.The load increment is defined as the CC of the wind energy.The idea of ELCC is expressed as: where I(•) is LOLE, EENS or other reliability indices.G and L are the total generation capacity and load demand of the whole system, respectively.∆G stands for the wind energy to be evaluated.∆L is the load increment that the power system can support due to the integration of the wind energy.
Therefore, we need to find the maximum value of ∆L, which is usually implemented by the secant method [11].The schematic diagram of the secant method is shown in Figure 3. ΔG stands for the wind energy to be evaluated.ΔL is the load increment that the power system can support due to the integration of the wind energy.
Therefore, we need to find the maximum value of ΔL, which is usually implemented by the secant method [11].The schematic diagram of the secant method is shown in Figure 3.
Input initial data and set r = 0 In the secant method, we need to compute the reliability index of the original system, which excludes the wind energy to be evaluated.The value of this index is denoted as R0.Then we add the wind energy to form the new system and start an iterative process to compute the CC value.
From the schematic diagram, one can perceive that the points Xr and Xr+1 in the secant method are approaching the terminal point, whose ordinate is R0.Actually, once the position of the terminal point is located, its abscissa minus L is the CC of the wind energy.
The procedures of the secant method are detailed below.
Step 1: Set r = 0. Determine two points X0 and B by computing the reliability index with the load demands of the new system being L and L + ΔG, respectively.
Step 2: Connect Xr and B with a straight line, in which we search for a point whose ordinate is R0.The abscissa of this point is defined as Lr+1.
Step 3: Find the new point Xr+1 by computing the reliability index with the load of the new system being Lr+1.
Step 4: If the ordinate of Xr+1 is not close enough to R0, set r = r + 1 and return to Step 2. Otherwise, end the iteration.Output the value of the last Lr minus L, which is just ΔL in (47).ΔL is the CC of the wind energy.
The flowchart of the whole procedures is also presented in Figure 3.In the secant method, we need to compute the reliability index of the original system, which excludes the wind energy to be evaluated.The value of this index is denoted as R 0 .Then we add the wind energy to form the new system and start an iterative process to compute the CC value.

Numerical Tests and Discussion
From the schematic diagram, one can perceive that the points X r and X r+1 in the secant method are approaching the terminal point, whose ordinate is R 0 .Actually, once the position of the terminal point is located, its abscissa minus L is the CC of the wind energy.
The procedures of the secant method are detailed below.
Step 1: Set r = 0. Determine two points X 0 and B by computing the reliability index with the load demands of the new system being L and L + ∆G, respectively.
Step 2: Connect X r and B with a straight line, in which we search for a point whose ordinate is R 0 .The abscissa of this point is defined as L r+1 .
Step 3: Find the new point X r+1 by computing the reliability index with the load of the new system being L r+1 .
Step 4: If the ordinate of X r+1 is not close enough to R 0 , set r = r + 1 and return to Step 2. Otherwise, end the iteration.Output the value of the last L r minus L, which is just ∆L in (47).∆L is the CC of the wind energy.
The flowchart of the whole procedures is also presented in Figure 3.

Simulation Setting
All tests are coded with MATLAB 2016a and run on an Intel Core i7-5500U personal computer with 16 G memory.We use LOLP and EPNS as the reliability indices.
In the test cases, the MCS method and traditional CE-IS method are compared with the DRCE-IS method in a modified IEEE-RTS 79 system.The data of the original IEEE-RTS 79 system can be acquired in [42].The loads at all buses are expected to be perfectly positively related.
Four wind farms are integrated into the modified system.We arbitrarily connect them to Bus 1, 2, 13 and 21, respectively.Each of them includes 50 3MW-WTGs of the same type, whose parameters are shown in Table 1 [43].

Parameter
Value Parameter Value WTGs may work under operation, derated or failure states.λ f , µ f , λ d , µ d , λ df and µ df are the transition rates between any two of the three states, as shown in Figure 4. Theoretically, the WTG can transition from any of the three states to another, as shown by Figure 4.However, according to the "partial failure mode" developed in Page 26-27 of [44], WTGs are always restored to the operation state rather than the derated state in the actual repair process, so it is difficult to obtain λ df and µ df by statistics.Hence, the transfer between the derated and failure states is always omitted.We also employ the assumption and set λ df and µ df as zero here.WTGs may work under operation, derated or failure states.λf, μf, λd, μd, λdf and μdf are the transition rates between any two of the three states, as shown in Figure 4. Theoretically, the WTG can transition from any of the three states to another, as shown by Figure 4.However, according to the "partial failure mode" developed in Page 26-27 of [44], WTGs are always restored to the operation state rather than the derated state in the actual repair process, so it is difficult to obtain λdf and μdf by statistics.Hence, the transfer between the derated and failure states is always omitted.We also employ the assumption and set λdf and μdf as zero here.

Derated state
Failure state Under the assumption, the probabilities of failure and derated states are derived by: where θf is the failure probability and θd is the derated probability, which are consistent with the notation in (37).
According to the data in Table 1, we can calculate θd and θf, which are 0.1050 and 0.1074, respectively.
According to the practical experience in the wind farms in Northwest China [43], the output of WTGs under derated state is set as 60-70% of its real-time maximum available output.Hence, we set Under the assumption, the probabilities of failure and derated states are derived by: where θ f is the failure probability and θ d is the derated probability, which are consistent with the notation in (37).
According to the data in Table 1, we can calculate θ d and θ f , which are 0.1050 and 0.1074, respectively.
According to the practical experience in the wind farms in Northwest China [43], the output of WTGs under derated state is set as 60-70% of its real-time maximum available output.Hence, we set the derated ratio as 60%, shown in Table 1.
Wind speed data at Bus 1, 2, 13 and 21 are extracted from the historical data of four areas in Northwest China from 1 January to 31 December in 2013.The areas are Ganhekou, Changma, Liuyuan and Mazongshan.The wind speed PDFs are shown as curves in Figure 5.The data can be referred to by DOI: 10.13140/RG.2.2.29281.97125.

Validation of SNMM
Before using SNMM and D-vine copulas to provide scenarios for reliability assessment, we need to validate them.Firstly, we use the wind speed data at Mazongshan as an example to illustrate the performance of SNMM.
The fitted PDF by using GMM and SNMM with eight components is depicted in Figure 6, respectively.A distinct difference between the two PDFs is marked by a red dashed circle, which shows that SNMM performs better than GMM when they have the same number of components.
Wind Speed / m•s -1 To present the superiority quantitatively, we compute the Chi-square values of the models fitted by SNMM, GMM [45], Weibull model and kernel density estimation (KDE) [46], as shown in Table 2.

Validation of SNMM
Before using SNMM and D-vine copulas to provide scenarios for reliability assessment, we need to validate them.Firstly, we use the wind speed data at Mazongshan as an example to illustrate the performance of SNMM.
The fitted PDF by using GMM and SNMM with eight components is depicted in Figure 6, respectively.A distinct difference between the two PDFs is marked by a red dashed circle, which shows that SNMM performs better than GMM when they have the same number of components.
to validate them.Firstly, we use the wind speed data at Mazongshan as an example to illustrate the performance of SNMM.
The fitted PDF by using GMM and SNMM with eight components is depicted in Figure 6, respectively.A distinct difference between the two PDFs is marked by a red dashed circle, which shows that SNMM performs better than GMM when they have the same number of components.To present the superiority quantitatively, we compute the Chi-square values of the models fitted by SNMM, GMM [45], Weibull model and kernel density estimation (KDE) [46], as shown in Table 2.

Wind Speed
Apparently, SNMM holds the minimum Chi-square values among the three for all wind speed.Hence, SNMM is validated and can be further used to model the marginal PDFs of wind speed and load in reliability assessment.To present the superiority quantitatively, we compute the Chi-square values of the models fitted by SNMM, GMM [45], Weibull model and kernel density estimation (KDE) [46], as shown in Table 2.
Apparently, SNMM holds the minimum Chi-square values among the three for all wind speed.Hence, SNMM is validated and can be further used to model the marginal PDFs of wind speed and load in reliability assessment.Since the load at all buses is perfectly positively correlated, there are five correlated RVs in the test system, i.e., the wind speed at four sites and the total load of the system.
We apply D-vine copulas to construct the correlation model based on historical data.Each bivariate copula can be chosen from Normal, Student t, Clayton, Gumbel, Frank and Joe copulas, whose parameters are shown in Table 3.To evaluate the goodness-of-fit of D-vine copulas, we use the statistic recommended by [47], which is based on Cramér-Von Mises (CvM) criterion and empirical copula process (ECP).We entitle it ECP-CvM, which indicates the difference between the fitted model and the empirical model.
The ECP-CvMs of D-vine copulas, C-vine copula [27], Normal copula and t copula are calculated and shown in Table 4.According to the results in Table 4, D-vine copulas exhibit superiority over Normal and t copulas due to a much smaller ECP-CvM value.It performs much better than any individual multivariate copula model, so we can adopt it for reliability assessment of power systems.

Validation of DRCE-IS
In this section, we validate the DRCE-IS method by comparing the values of LOLP and EPNS obtained by MCS, CE-IS [28] and DRCE-IS.Then the CC for wind farms will be evaluated by the DRCE-IS and the secant method.
Relevant parameters of DRCE-IS are set as follows: (1) n s in Step 1 of the pre-simulation stage is 50,000.
(2) ρ in the Step 4 of the pre-simulation stage is 0.05.
(3) The relative error limit in the main simulation stage is set as 0.05 for both LOLP and EPNS.
Based on the parameters above, we can assess the reliability of power systems with the four integrated wind farms, whose values are displayed in Table 5.Since the LOLP and EPNS values obtained by the three methods are remarkably close to each other under the same relative error limit 0.05, we conclude that DRCE-IS is correct.Besides, the sample size of DRCE-IS is smaller than that of CE-IS, because of the proposed DR technique.Hence, the DRCE-IS is validated.
The convergence process of the three methods is depicted in Figure 7.The ordinate is the relative error of the LOLP computed by samples.With the increase of size, the errors of the three all decrease.One can perceive that the relative error of DRCE-IS is the first to decrease to 0.05, which is also consistent with the results in Table 5.
sample size of DRCE-IS is smaller than that of CE-IS, because of the proposed DR technique.Hence, the DRCE-IS is validated.
The convergence process of the three methods is depicted in Figure 7.The ordinate is the relative error of the LOLP computed by samples.With the increase of sample size, the errors of the three all decrease.One can perceive that the relative error of DRCE-IS is the first to decrease to 0.05, which is also consistent with the results in Table 5.

Relative Error
Sample Size/10

Evaluate the Capacity Credit of Wind Energy by DRCE-IS and the Secant Method
Finally, we can calculate the CC of the four wind farms with the secant method and all the methods validated above.The parameters of DRCE-IS is exactly the same as those in Section 6.2.3.
According to Figure 3, we need to compute the reliability indices in three distinct situations to start the iteration of the secant method, which correspond to Rd, Ru and R0 in Figure 3.
We employ LOLP to compute Rd, Ru and R0, whose values are listed with the CC of each wind farm in Table 6.
Besides, we also compute the ratio of CC of each wind farm to its installed capacity, which is called capacity credit factor (CCF).Since the installed capacity of each wind farm is 3 × 50 = 150 MW, we can easily obtain the CCF values, shown in Table 6 as well.
The results in Table 6 unveil the fact that the CCF of wind energy is rather small (below 30%), which reflects a low efficiency in utilizing wind energy.This is definitely due to the characteristics of uncertainty and fluctuation of wind.Therefore, how to overcome the disadvantages to improve the CCF of wind energy by available measures, such as by energy storage systems and demand response dispatch, will be researched in our future work.Finally, we can calculate the CC of the four wind farms with the secant method and all the methods validated above.The parameters of DRCE-IS is exactly the same as those in Section 6.2.3.
According to Figure 3, we need to compute the reliability indices in three distinct situations to start the iteration of the secant method, which correspond to R d , R u and R 0 in Figure 3.
We employ LOLP to compute R d , R u and R 0 , whose values are listed with the CC of each wind farm in Table 6.
Besides, we also compute the ratio of CC of each wind farm to its installed capacity, which is called capacity credit factor (CCF).Since the installed capacity of each wind farm is 3 × 50 = 150 MW, we can easily obtain the CCF values, shown in Table 6 as well.
The results in Table 6 unveil the fact that the CCF of wind energy is rather small (below 30%), which reflects a low efficiency in utilizing wind energy.This is definitely due to the characteristics of uncertainty and fluctuation of wind.Therefore, how to overcome the disadvantages to improve the CCF of wind energy by available measures, such as by energy storage systems and demand response dispatch, will be researched in our future work.

Conclusions
This paper focuses on the capacity credit evaluation of wind energy, whose key points are modeling the uncertainty of the correlated wind speed and load and assessing the power system reliability.The main achievements and conclusions are listed below: (1) A novel method of fitting the joint PDF is proposed for wind speed and load.Particularly, SNMM is used to fit the marginal PDFs, while the correlation model is established by D-vine copulas.(2) The goodness-of-fit of SNMM, GMM, Weibull model and KDE is quantified by the Chi-square values, among which that of SNMM is the smallest.Similarly, the goodness-of-fit of D-vine, C-vine, Normal and t copulas is quantified by the ECP-CvM values, among which that of D-vine is the smallest.Therefore, SNMM and D-vine copulas are validated.

FindFigure 3 .
Figure 3. Illustration of the secant method.(a) schematic diagram of the secant method; (b) flowchart of the secant method.

Figure 3 .
Figure 3. Illustration of the secant method.(a) schematic diagram of the secant method; (b) flowchart of the secant method.

Figure 5 .
Figure 5. Wind speed PDFs of four areas.

Figure 6 .
Figure 6.Fitted PDFs of the wind speed at Mazongshan.

Figure 5 .
Figure 5. Wind speed PDFs of four areas.

Figure 6 .
Figure 6.Fitted PDFs of the wind speed at Mazongshan.

Figure 6 .
Figure 6.Fitted PDFs of the wind speed at Mazongshan.

Figure 7 .
Figure 7.The convergence of the three methods.

Figure 7 .
Figure 7.The convergence of the three methods.6.3.Evaluate the Capacity Credit of Wind Energy by DRCE-IS and the Secant Method

1 T P sh
T P g − 1 T P d + 1 T P sh = 0 P gmin ≤ P g ≤ P gmax 0 ≤ P sh ≤ P d −P linemax ≤ bAδ ≤ P linemax A T bAδ = P g + P sh − P d ,

Table 2 .
Chi-square values of three models.

Table 2 .
Chi-square values of three models.

Table 3 .
Parameters of the constructed D-vine copulas.

Table 4 .
Comparison of ECP-CvMs of three copula models.

Table 5 .
Reliability indices of the modified IEEE-RTS 79 system with four dependent wind farms.

Table 6 .
Three LOLP values and CC of four wind farms.

Table 6 .
Three LOLP values and CC of four wind farms.