Risk-Based Probabilistic Voltage Stability Assessment in Uncertain Power System

The risk-based assessment is a new approach to the voltage stability assessment in power systems. Under several uncertainties, the security risk of static voltage stability with the consideration of wind power can be evaluated. In this paper, we first build a probabilistic forecast model for wind power generation based on real historical data. Furthermore, we propose a new probability voltage stability approach based on Conditional Value-at-Risk (CVaR) and Quasi-Monte Carlo (QMC) simulation. The QMC simulation is used to speed up Monte Carlo (MC) simulation by improving the sampling technique. Our CVaR-based model reveals critical characteristics of static voltage stability. The distribution of the local voltage stability margin, which considers the security risk at a forecast operating time interval, is estimated to evaluate the probability voltage stability. Tested on the modified IEEE New England 39-bus system and the IEEE 118-bus system, results from the proposal are compared against the result of the conventional proposal. The effectiveness and advantages of the proposed method are demonstrated by the test results.


Introduction
With the rapid growth of power demand, modern power grids have been changing incredibly throughout the world.Despite the sustainable features, wind power generation introduces uncertainties to the entire power system.Additionally, load forecasting errors with respect to the change of operating conditions also contribute to the system uncertainties.Therefore, the risks of bearing these uncertainties during operating states must be quantified and fully investigated.The security of the voltage stability problem considered in this paper is motivated by the increasing presence of stochastic elements in the power system as a result of integrating intermittent renewables at both the wholesale and retail levels.The surge of solar power integration in recent years, for example, has fundamentally changed the overall net load characteristics.In some areas, the traditional load profile is being transformed to the so-called "duck curve" profile where a steep down-ramp during the hours when a large amount of solar power is injected into the network is followed by a steeper up-ramp when the solar power drops sharply in the late afternoon and early evening hours.
While the duck curve phenomenon represents an average net load behavior, it is the highly stochastic and spatial-temporal-dependent ramp event that presents difficult operational challenges to system operators.For this reason, power system voltage stability is one of major topics of concern for secure operation of transmission networks.It is defined as the ability of a power system to maintain load voltage magnitudes after disturbances [1].The static method is widely analyzed for voltage stability assessment when large demand and limited network expansion occur.The relationship between transmitted power and receiving end voltage (PV) curve is one of the powerful tools for stability analysis [2].The work in [3] provides an index that accounts for both the future uncertainties on the system and the consequences associated with voltage collapse and violation of limits.Besides, Reference [4] proposes an on-line voltage security method based on wide-area measurements and a decision tree.While the operating point reaches the critical voltage collapse point, the system may lose controllability and face a series of cascading outages problem.Furthermore, due to the load fluctuation and intermittent energy behavior, it is hard to predict the system operational status.The voltage stability analysis should consider how to accurately express the system's risks.
Probabilistic risk-based approaches have been adopted for not only power system planning, but also the security assessment for the last decades [5][6][7].The work in [8] provides an AFTER (A Framework for electrical power sysTem vulnerability identification, dEfense and Restoration) in-depth risk assessment tool to deal with extreme events for static security assessment, such as threats, vulnerability, contingency and impact.Otherwise, Reference [9] introduces a transmission line overload risk index for measuring the security level of wind integrated power systems considering the wind and load-power generation correlation.Besides, a fuzzy risk assessment of small-disturbance stability has been proposed in [10] to handle the rotor angle instability.The literature for voltage risk assessment is rather limited.In order to quantitatively evaluate the voltage stability risk, Reference [11] provides a stability index for determining the bifurcation point of the power system.The voltage instability probability is calculated based on a joint distribution for the entire system [12,13].The work in [14] addresses the voltage stability topic, including unstable states caused by modeling resolvability and loss of controllability.The risk of high distribution generation penetration in a low-voltage secondary distribution network affecting voltage stability is analyzed in [15].The work in [16] proposes a set of risk indices to tell the weak lines and buses of the system, which can bring in static voltage instability, and the author also provides an assessment technique to calculate them.Furthermore, risk can also be used for the application of security-constrained optimal power flow to describe the system health [17,18].
The stochastic nature of the system uncertainties requires the probabilistic approach to cope with them.These techniques could be classified under three main categories: cumulant method, point estimate method and Monte Carlo (MC) simulation.The cumulant method is widely used for solving probabilistic power flow [19][20][21], which is computationally effective, but requires some mathematical assumptions in order to simplify the problem.The cumulant method combined with the Gram-Charlier expansion is used to estimate the output distribution [20,22].Later, the method of cumulants and the von Mises function is developed to deal with discrete distributions [21].Furthermore, point estimate methods provide an approximate statistical description of the output random variables with less computational burden [8,9,23].The Monte Carlo simulation is widely used as the benchmark method for the reason that circuit size, simulation time, etc., often make Monte Carlo analysis time expensive [24].In this paper, the Quasi-Monte Carlo (QMC) simulation is made to accelerate the speed of Monte Carlo simulation by improving the technique of sample generation.
Although researchers have made great efforts, these probabilistic risk-based methods are not appropriate enough for large wind power-penetrated system, because these indices cannot reflect the tail risk of the voltage static loss.Furthermore, the randomness and violation are not considered in their forecast model.Hence, a risk-based probabilistic voltage stability assessment is of great value to be learned.As we know, Conditional Value at Risk (CVaR) is used as risk hedging tool in a power system for its advantages of consistency and subadditivity [25][26][27].Sometimes, it is difficult to calculate CVaR for the reason that it is hard to find an analytical representation of the loss function, which is associate with the distribution of the input random variables [28].Hence, Monte Carlo simulation is employed to the optimization and calculation of the CVaR to solve this problem.However, the risk-based security assessment is a state-of-the-art method, which will be used in the near future.The Monte Carlo simulation cannot handle the online computation, because of the drawback of the great number of simulations required to obtain convergence [29].We provide a new voltage stability index based on QMC and CVaR, which reveals important aspects of static voltage stability and is used for the online evaluation of tail risk-based security levels of wind-penetrated power systems.
The paper is organized as follows.In Section 2, the probabilistic model of wind power generation is built by the Gaussian mixture model, whose parameters are obtained according to the Expectation Maximization (EM) algorithm.In Section 3, a security risk index of voltage stability based on CVaR and the PV curve is described.In Section 4, the QMC simulation is adopted to forecast the real-time voltage distribution.In Section 5, two comprehensive case studies are designed, and the corresponding results are analyzed to highlight the advantages of the proposed method.

Probabilistic Model of Wind Power Generation
Wind power ultra-short-term forecast error is usually assumed to follow a weighted sum of beta distribution.The work in [30] uses the Gaussian mixture distribution to model any non-Gaussian distribution.Furthermore, in order to represent the stochastic power output of wind farms, References [31][32][33] use the Gaussian mixture model to approximate the development for the probability density function (PDF) of wind power generation.The ultra-short-term wind power forecast error distribution, following the Gaussian mixture model at each time interval, had also served in in this paper.Hence, the probability density function of measured wind power can be modeled by a Gaussian mixture model corresponding to different given predicted wind powers.The observed value of the measured wind power is available from the stratified sampling of the wind power scenario tree for a single time interval.According to the observed value, maximum likelihood mixture densities' parameters can be found via the EM algorithm [34].In addition, the Stochastic Approximation Expectation Maximization (SAEM) algorithm [35], compared to the EM algorithm, performs better and converges almost surely to a local maximum of the likelihood of the function.In this paper, we obtained 140,160 sampling data from the four years, each of the 15-min wind power historical data on [36], the website of Belgium's transmission system operator in Europe.The histogram of 15-min forecast error compared with the Gaussian mixture distribution is plotted in Figure 1, which shows that the Gaussian mixture distribution can fit the shape of the forecast error PDF.The PDF and parameters of the Gaussian distribution are calculated by the following equations: (1) where f (x) is the PDF of the total forecast error and w i is the weight of the forecast power bin i.

Security Risk Evaluation of Voltage Stability Issues
IEEE defines voltage collapse as the process by which voltage instability leads to the loss of voltage in a significant part of the power system.Voltage stability is a complex phenomenon that involves many power system components at once.As mentioned in [37], the voltage collapses take place on different time scales ranging from seconds to hours: (1) electromechanical transients (e.g., generators, regulators condensers) and power electronics (e.g., SVC (Static Var Compensator) , HVDC (High-Voltage Direct Current) reaching reactive power limits) in seconds; (2) discrete switching devices (e.g., action of tap-changers transformers) acting at tens of seconds; (3) load recovery processes in several minutes.Besides, increasing in loading, line tripping or generator outages are also contributing to the voltage collapses, and the method proposed in this paper captures partially of such a complex phenomenon in minutes time scales.
Nowadays, it is of high value to evaluate the safety distance from a normal operational status to a process of cascading trip incident as typical centralized wind farms are integrated.However, the operational status will be estimated to be a probability distribution, instead of a determined point, in the near future when considering the wind uncertainties.In this section, a new index of the risk-based voltage stability margin has been proposed for security risk evaluation of voltage stability analysis when considering the small probability of a serious event in an uncertain environment.

Definition of Conditional Value at Risk in Security Risk Evaluation
Risk is the potential of gaining or losing something of value [38].In this paper, risk is associated with 'security losses' of the locational voltage of the power system.Risk assessment evaluates the potential voltage stability loss for system operators.
Various methods are used to measure risks.Specifically, CVaR is an index to represent the expected shortfall of the profit distribution at a specified confidence level in finance and banking [39].Risk measures used for decision making in electricity markets are discussed in [40,41], and the advantages of CVaR are explained.It is considered to be a more consistent measure of risk than Value-at-Risk (VaR) [28].
Taking into account the uncertainties associated with electric load and wind power generation integrations, the system may face the risk that the static voltage approaches the voltage collapse point V cp , which is the nose of the PV curve (as shown in Figure 2).
When the system situation reaches the voltage collapse point, a series of uncontrollable cascading events, such as a blackout, may come out step by step.Besides, the voltage collapse value can be obtained by the continues power flow.The question is how to evaluate the risk of the static security for the power system in an uncertain environment.We define the static voltage stability loss, which is shown as follows: where Ψ(x, α) is the cumulative distribution function (CDF) of node voltage stability loss for extreme conditions with respect to demand level x.Furthermore, it is the probability of the relative distance from the node voltage operation point to the collapse point f (x, y) without exceeding a threshold α.
Besides, ρ(y) is the probability density of variance y for f (x, y).
The β-VaR and β-CVaR values for the loss random variable following x and any specified probability level β in (0, 1) will be denoted by VaR β (x) and CVaR β (x): and: In Equation ( 5), the VaR function, which is the percentile of the loss distribution with confidence level β, is the smallest number, such that Ψ(x, α) = β.In Equation ( 6), the CVaR function is defined as the conditional expected loss associated with x relative to the loss higher than or equal to VaR with probability level β (as shown in Figure 3).It is difficult to calculate CVaR because it is hard to find an analytical representation of VaR, which has been involved in CVaR's definition.Without reference to VaR, a much simpler way of defining CVaR is:

Loss function Probability density
CVaR where [ ] + = max( , 0).However, instead of obtaining an analytical expression for the implementation of the approach, it is sufficient to generate random samples y = {y m |y 1 , y 2 , • • • , y M } from the node voltage magnitude distribution via QMC simulation: where M is the number of QMC simulations.

Index of Voltage Stability Margin Considering the Risk
The steady state relation between injected real power and voltage at a particular bus is described by the PV curve [42] (as shown in Figure 2).It provides sufficient information on the active power margin and voltage magnitude due to small disturbances.The PV curve helps system operators in reliability assessment.
For the purpose of locating the weakest bus of a system, the traditional local voltage stability margin (V LVSM ) is measured as the relative distance of the voltage magnitude from the operating point (V no ) to the voltage collapse point (V cp , the nose of the PV curve).Furthermore, V LVSM admits a value between zero and one.It is defined as: As we consider more and more renewables injected in the power system, the state-of-the-art operational status estimation is deemed to be a probability distribution instead of a deterministic point.Via the results of the probabilistic power flow by the QMC algorithm (Section 4), the probability density functions of the node voltage magnitude can be achieved by modeling the stochastic variation of system uncertainties, such as wind power generation and load fluctuation.We provide a powerful risk-based voltage stability margin index (V RBVSM ) to measure the security risk of voltage stability loss in the near future.The voltage stability loss is defined as V no − V cp , and V RBVSM is the CVaR β of node voltage stability loss at a specified probability level β.Therefore, the tail risk can be taken into account.
Furthermore, the CVaR-based voltage margin index V RBVSM is a linear sensitive index.It measures the relative distance from the current operational point to the voltage collapse point.This index is of great value for warning about dangers in the proximate future, while it will not change suddenly or sharply until access to the voltage collapse point.It provides ample time for the system operators to make preventive control decisions.Besides, the index is useful for identification of the weakest bus of a system, which helps the system operators focus on these weakest buses, and it can also be used for the application of security-constrained optimal power flow to quantify the system security risk.
Besides V RBVSM , which can be used to measure the risk of voltage stability loss for each local buses, we also provide another voltage stability index, the probabilistic voltage stability margin (V PVSM %).As if there is not only one bus that is unstable in the system, the indices can help the system operators to identify the weakest area in the power system for the state-of-art forecast.The n% V PVSM % is defined as the n% percentage of buses that are voltage stable buses for the near future.
where N s is the number of voltage stable buses; each indicator function of voltage stable bus N s (i) for bus i is calculated in (11); N b is the total number of buses in the power system.They are all estimated for a particular time interval in the near future.Equations ( 10) and (11) are not the only indices, but help the operators make decisions for preventive control and planning ahead of time.

A Quasi-Monte Carlo Simulation for Real-Time Probabilistic Voltage Forecast
In this section, we provide the key theoretical foundation of the probabilistic sampling technique QMC for the real-time static voltage forecast assessment.Instead of generating the pseudo-random samples, the QMC simulation focuses on obtaining a set of low-discrepancy sampling points.

Quasi Monte Carlo Simulation
Quasi-Monte Carlo methods, using a deterministic set of points from low-discrepancy sequences, aim to outperform the classical Monte Carlo method [43][44][45][46].The goal for almost all of the MC simulations is to converge to the required accuracy as rapidly as possible, with the least sample simulations.Although pseudo-random sampling points typically used by classical MC are uniformly distributed, geometrically, they suffer from clusters and empty spaces, which are unequally separated over the sampling region.Rather than randomly generating the samples, the QMC tries to select the point set homogeneously equidistantly distributed over the space.
The MC integration error |Q − Q n | decreases asymptotically to zeros as convergence rate O(n − 1 2 ) for general integrable f with increasing n [47].QMC is used to improve on this convergence behavior, which targets the behavior n − 1 2 and not the proportionality constant σ( f ) as the standard variance reduction techniques.
There is another way to analyze the integration error, based on Niederreiter's development of the concept of discrepancy, which is used to measure the uniformity for finite point sets.
Remark 1 (Discrepancy [44]).Let C s denote the s -dimensional unit cube.For a Jordan measurable set J ⊆ interval C s and a sequence of N points {x n } in C s , let Vol(J) be the volume of J, and let I J (x) be the characteristic function of J. Define the n J as the number of points lying inside J, Then, we define the discrepancy as: If we only look at hyper-rectangles with one corner at zero, J = [0, a) where a ∈ C s , the star discrepancy is defined as: Reference [48] shows that for uniform points x i ∼ U[0, 1) s , the star discrepancy for classical MC is 2 ).However, QMC performs using the samples generated from a deterministic sequence; the low-discrepancy sequence, which has been constructed with star discrepancy bound ), yielding a convergence rate of O( (logN) s N ), shows the surprising attribute of asymptotically superior discrepancy, in comparison with the pseudo-random sequence used by classical MC.
Finally, Koksma-Hlawka gives the answer to the question of the relationship of the integration error in MC and the discrepancy of a point set.
Theorem 1 (Koksma-Hlawka [44]).If function f has bounded variation in the sense of Hardy and Krause, then the Monte Carlo error is bounded as follows.
where V( f ) = 1 0 |d f | is the variation of f and D * N is the star discrepancy of the sequence {x n }.
By Theorem 1, the Koksma-Hlawka inequality function guides us to search for effective point sets and sequences by making precise the role of star discrepancy.It is particularly attractive because using the points with lower star discrepancy D * N might reduce the integration error.Obviously, we see the terrific possibility of linear (N −1 ) convergence of QMC integration error, instead of the slow (N −0.5 ) MC integration error [24].
In the field of computational finance, QMC has had widespread success and obvious application [43].Besides, References [24,49] handle it for the research of statistical static timing analysis, small signal stability analysis and other engineering problems.All of the results show significant gains over classical MC.In this paper, we will try to use it for the calculation of probabilistic power flow.

The Formulation of the Real-Time Voltage Forecast Problem
In an uncertain environment, with stochastic and intermittent variables injected into the power system, the forecast of node voltage magnitude is not truly deterministic.Most quantities of practical interest present it as probability distributions.Hence, the QMC simulation technique is a state-of-the-art prediction strategy for probabilistic voltage forecast analysis.
The real-time statistic voltage and its posterior probability mass function at time T forecast from initial time t can be obtained by the probabilistic power flow, which is defined as, where j is the number of buses, V i(T) is the discrete value of voltage that may occur at time T and M is the number of QMC runs.The indicator function 1{•} is set to be one, while the voltage at bus j is V i(T) .The computation of the above conditional probabilities depends on a probabilistic model of the power injections, such as wind power generation The relationships of voltage state V i , W and D are described in Equation ( 17).
The technique of marginal CDF transformation [23] maps the problem into the independent standard uniform space.Hence, the voltage state can be regarded as a function of independent standard uniform random variables: where u ∼ U(0, 1), and g i is a piecewise constant function, when we assume that the system topology is fixed in each time interval.The step function g V i , having only a finite number of pieces, can be written as: where

is the number of intervals over the domain of Û.
A k is the k-th interval of the piecewise constant function g i .Their relationships are described in the following indicator function 1{•}: The intervals A k are assumed to have the following properties: (1) The intervals are disjoint.For all a = b, A a ∩ A b = ∅; (2) The union of the intervals is a multi-dimensional integration over the entire m-dimensional unit hypercube [0, 1] m .
For example, the Lebesgue integration of the piecewise constant function g V i ( Û) can be calculated as: where l(A k ( Û)) is the length of the interval A k ( Û), which satisfy Û ∈ A k .
For QMC simulation, it can be obtained as: where {U h , (h = 1, 2, • • • , u)} is the low-discrepancy sequence in QMC simulation.In a typical MC simulation setting, in order to obtain accuracy and flexibility, we suffer from the cost of speed: a huge number of random samples {U h } is needed to guarantee a higher accuracy.In fact, the convergence rate of MC integration error is as high as N −1/2 for general integral function f as we run multiple N-point MC runs.Additionally, the proof is based on the central limit theorem [24].This means that, in number theoretic investigations of deterministic sampling patterns, the integral result will be more accurate when the point set samples are more uniformly distributed in each unit hypercube, which can be achieved when N increases in MC simulation.Consequently, the QMC technique is a choice to gain extraordinary accuracy by improving the sampling generation properties.The low-discrepancy sequences {U h } are generated by QMC, instead of the pseudo-random sequences, to satisfy the accurate evaluation of the Lebesgue integration [0,1] m g V i ( Û) d( Û) with a lower computation burden.
According to [47], the convergence of QMC integration error increases to as high as N −1 .
Based on the sample sets of voltage state V i in each time interval [t, t − 1), the probability of node voltage in each interval can be achieved.Furthermore, the whole process of the risk-based static voltage evaluation assessment will be adopted in the next Section 4.3.

The Process of the Proposed Algorithm
The process of the risk-based voltage stability assessment is shown as below: (1) Obtain the wind power forecast error distribution based on historical data; (2) Obtain the PV curve and voltage collapse value through continues power flow; (3) Generate sample sets for the forecast time interval via the QMC algorithm; (4) Calculate V RBVSM to evaluate the tail risk of the voltage stability loss for each bus; (5) Calculate V PVSM % to evaluate the voltage stability level for the entire system; (6) Take the value of the weakest buses for further dispatch adjustment.

Results and Discussion
In this section, we present some results to compare the performance of the proposed probabilistic static voltage assessment algorithm with the certainty equivalent method.We first build the probabilistic wind power model, then test our algorithm in the IEEE 39-bus system to obtain its behavior in the basic situation and at last test how static voltage reliability is affected in the IEEE 118 bus system.The proposed algorithm's effectiveness and advantages are demonstrated with the test results.The simulations are programmed based on Matpower 4.0 in Matlab (MathWorks, Natick, MA, USA) , and the test are fulfilled on an laptop with an Intel Core i7-4790 at 3.1-GHz and 32 GB memory (Intel, Santa Clara, CA, USA).

Statistical Distribution of Wind Power Forecast Error Analysis
The data analyzed in this section are obtained from the amended Renewable Energy Sources (RES) Act (EEG) in Germany [50].The dataset contains a two-year time series of 15-min projection and the actual value of wind power for 10 wind generations in a wind farm, with around one million data available for analysis, which correspond to the years 2014 and 2015.
Figure 4a,b shows the PDF of forecast error corresponding to the prediction horizon of 1 and 5 h.In the short-term forecast horizon, a particular normal distribution is no longer fitted for the actual forecast error, which is usually used in the probabilistic analysis.Specifically, the shorter the forecast horizon is, the higher the kurtosis is than the others.
The EM algorithm is employed to the Gaussian mixture model, which generates the maximum likelihood estimation of parameters in a statistical model.The clusters of the Gaussian mixture wind power forecast error model for wind generations G N1 and G N2 located geomagnetically different for a 1-h and 5-h forecast horizon are presented in Figure 5a,b.Furthermore, the PDF of the 1 h-ahead forecast error Gaussian mixture distribution for G N1 is shown in Figure 6.The cluster for 1 h-ahead forecast gathers more closely than the clusters for the 5 h-ahead forecast, which shows that the confidence level of 1 h-ahead forecast error is smaller than that of the 5 h-ahead forecast error.The parameters for the function of the 1 h-ahead forecast error Gaussian mixture distribution for 10 wind generators are shown in Table 1.G N1 has the largest variance in the average level, which means it has the biggest volatility compared to the other wind generations.Figure 7 describes the 3D clusters of the Gaussian mixture model of wind power forecast error for G N1 , G N2 and G N3 for a 1-h forecast horizon.It shows that in the view of G N1 , the point sets are distributed more discretely than in other directions.

Case Study: IEEE 39-Bus System
The performance of the proposed algorithm is tested on the IEEE 39 bus New England system [51] in this section.This system consists of 39 buses, 46 lines, 10 units terminal and 29 loads.The wind farms are located at Buses 14, 15 and 30.The rated capacity of each wind farm is 500 WM.Wind generations are fitted for the Gaussian mixture distribution, from which the parameters are obtained in Section 5.1.The network topology of the modified IEEE 39 bus system is given in Figure 8. Figure 9 shows the comparison between the MC and QMC 1000 times simulation for the probabilistic forecast error model for wind generation G N1 located at Bus 14.We can figure out that the result of QMC in Figure 9a has less barbs to fit for the Gaussian mixture distribution, which shows that it is closer to the estimated model and has a better performance approach to the actual distribution.However, the results of MC in Figure 9b are worse, which need to sacrifice the efficiency of the computational cost to gain a better performance.The shape of the magnitude voltage distribution is highly different at different locational buses (in Figure 10), from the perspective of the width of the voltage distribution.The width of the voltage distribution at Bus 15 is three-times larger than the width of the voltage distribution at Bus 30.The peak value of the magnitude at Bus 30 is the most centralized one overall on the whole map of the voltage distribution in the 39-bus system.This is because the randomness comes from the renewable generation at Buses 14, 15 and 30; the voltage at these buses has a longer tail to simulate the worst cases of the small probability events.Besides, the peak value of the magnitudes at Buses 14, 15, 19, 30 is about 0.98675, 0.96188, 1.01094, 1.01143 V, respectively.Besides, the peak value of the phase at Bus 14 is about −52.13 radians.Most of the probability is clustered with respect to this range of phase value.However, the peak value of the phase at Buses 15, 19 and 30 is about −45.59, −42.45, −41.01 radians, respectively.In addition to this, the shape of the tail is also different for different buses.Even if there is no magnitude or phase out-of-limit event, the risk for each locational voltage is extremely different.The nose curve of the voltage magnitudes at Buses 14 and 20 is compared in Figure 11.The voltage collapse points are 0.710 and 0.915 V, respectively.Then, the index of VaR, CVaR and V RBVSM can be calculated to evaluate the safety-distance from the predicted operational point to the collapse point.The comparison results are shown in Table 2, the sequence of buses is in an ascending order according to the index of V RBVSM to evaluate the weakest voltage bus.The smaller the value of V RBVSM is, the larger the risk of voltage loss is.Hence, Bus 36 is the weakest bus, which may face the voltage out-of-limit event immediately.Additionally, Bus 8 is the strongest bus, which is the safest bus in the IEEE 39-bus system.We should pay more attention to the weakest buses in the system.Several precautions, such as wind curtailment, should be adopted ahead of time.

Case Study: IEEE 118-Bus System
A modified IEEE 118-bus system is used to study the proposed method.The IEEE 118-bus system has 54 thermal generators, 186 branches and 91 load buses.We assume that three wind turbines are located at Buses 15, 54 and 96.The system topology is given in [52].
The computational cost of the proposed algorithm is approximately 40.17 s via 200 times QMC simulation.However, the computational cost of the method complemented with the same convergence rate via MC runs is about 8000.52 s, which runs the MC simulation 40,000 times.The convergence rate compared between QMC and MC is shown in Figure 12a.
We make a comparison of our proposed method (Method 1) and a conventional risk-based approach (Method 2) [53].The work in [53] gives a process of voltage instability risk calculation for a given contingency: The indices of VaR, CVaR and V RBVSM of the voltage magnitudes are shown in Table 3 to evaluate the weakest 20 buses in the IEEE 118-bus system.The risk of voltage out-of-limit events at these buses is larger than the risk of voltage loss at any other buses.However, using Method 2, they only give one risk index value for measuring voltage instability in a particular case.The weakest bus is Bus 24, the same as the result using Method 1. Furthermore, the risk index of Bus 24 is equal to 11. Besides, the risk index for other buses is equal to zero.The index proposed in Method 2 is not as sensitive as the index using Method 1.As shown in Table 3, the ranking of the 20 weakest buses is listed to support the system operators to pay attention to these buses.Even though Method 2 can identify the weakest bus and using the risk index for measuring the weakest bus as the final result, it is hard for Method 2 to identify the weakest area or to help the system operators to perform preventive control as if the voltage instability may happen in a different area.We assume that the nodal load of the IEEE 118-bus system is fitted for the California 'duck curve'-nodal load curve as shown in Figure 12b.In order to evaluate the magnitudes at the weakest Bus 24, we obtain the index of V RBVSM as shown in Figure 13 according to the time varying wind power output forecast for each 5 min.As we know, when the value of V RBVSM is closer to zero, the probability of a voltage collapse event happening at that particular point is higher.In this scenario, the nodal load reaches the peak value at the last 1 h, as the value of V RBVSM decreases to the lowest level in the same time period.

Summary and Conclusions
In the operation of today's electricity power system, the on-line security of the operation is still the most significant issue, especially in an uncertain environment.This paper develops a risk-based statistic voltage evaluation mechanism.Different from the conventional deterministic evaluation method, the powerful risk-based voltage statistic index V RBVSM provides the relative safety distance to the voltage collapse point on the system risk level considering the tail risk of the voltage stability loss.Using this voltage evaluation index results in more secure operating conditions when the precaution has been made ahead of serious events happening.The new voltage evaluation approach was illustrated and compared in both the IEEE 39-bus system and the IEEE 118-bus system.In both cases, the proposed method provides a more conservative result depending on the forecast estimation for the near future.While the parameter of the wind power output forecast error distribution has been studied, the real-time probabilistic voltage forecast can be simulated based on the QMC sampling techniques to speed up the whole approach.Furthermore, the voltage index proposed in this paper offers an alternative approach based on the real-time probabilistic voltage forecast method.The approach avoids the effect of small probability serious events on voltage spikes while providing control to enhance the system security level.The proposed method is of great value in the application of real-time economic dispatch, system security prevention and control in the future.

ForecastFigure 1 .
Figure 1.Comparison of the histogram of 15-min forecast error with the Gaussian mixture distribution using historical data on [36].

Figure 2 .
Figure 2. The relationship between transmitted power and receiving end voltage (PV curve).

Figure 3 .
Figure 3.The comparison of VaR and CVaR.

ForecastFigure 4 .Figure 5 .
Figure 4. (a) A histogram of 1 h-ahead measured wind forecast error; and (b) a histogram of 5 h-ahead measured wind forecast error.

ForecastFigure 6 .
Figure 6.The Gaussian mixture wind power forecast error distribution for a 1-h horizon.

Figure 7 .
Figure 7.The three-dimensional clusters of the Gaussian mixture model of wind power forecast error for a 1-h horizon.

Figure 9 .
Figure 9.The comparison between the (a) QMC and (b) MC simulation of the Gaussian mixture wind forecast error model for wind generation G N1 at Bus 14 (1000 times).

Figure 10 .
Figure10.The voltage magnitudes distribution at different buses in the 39-bus system (voltage magnitudes/p.u. vs. probability density).

Figure 11 .
Figure 11.The nose curve of voltage magnitudes at Buses 14 and 20.

( 1 )Figure 12 .
Figure 12.(a) The convergence rate for QMC and MC sampling (0.005 is the rate of convergence); and (b) the duck curve of nodal load in the California power system.

Figure 13 .
Figure13.The index of V RBVSM at Bus 24 according to the duck curve nodal load at risk level β = 95%.

Table 1 .
The parameters of the 1 h-ahead forecast error Gaussian mixture distribution for 10 wind generations.

Table 2 .
The indices of VaR, CVaR and V RBVSM for evaluating the weakest 10 buses at the risk level β = 95%.

Table 3 .
The indices of VaR, CVaR and V RBVSM for evaluating the weakest 20 buses at the risk level β = 95%.