Risk Assessment for Distribution Systems Using an Improved PEM-Based Method Considering Wind and Photovoltaic Power Distribution

The intermittency and variability of permeated distributed generators (DGs) could cause many critical security and economy risks to distribution systems. This paper applied a certain mathematical distribution to imitate the output variability and uncertainty of DGs. Then, four risk indices—EENS (expected energy not supplied), PLC (probability of load curtailment), EFLC (expected frequency of load curtailment), and SI (severity index)—were established to reflect the system risk level of the distribution system. For the certain mathematical distribution of the DGs’ output power, an improved PEM (point estimate method)-based method was proposed to calculate these four system risk indices. In this improved PEM-based method, an enumeration method was used to list the states of distribution systems, and an improved PEM was developed to deal with the uncertainties of DGs, and the value of load curtailment in distribution systems was calculated by an optimal power flow algorithm. Finally, the effectiveness and advantages of this proposed PEM-based method for distribution system assessment were verified by testing a modified IEEE 30-bus system. Simulation results have shown that this proposed PEM-based method has a high computational accuracy and highly reduced computational costs compared with other risk assessment methods and is very effective for risk assessments.


Introduction
The extensive penetration of renewable-type distributed generators (DGs) (e.g., wind and PV) in distribution networks could bring many benefits to the grid, as they are alternative to conventional generators [1,2].However, the randomness of these DGs could cause some critical risks to security and economy aspects of distribution systems, such as power quality and stability, fault level, and the value of load curtailment, which impose challenges when planning distribution systems [3,4].Therefore, it is becoming increasingly important to assess the risks associated with the variability and uncertainty of DGs.
In the past decade, many researches concentrated on how to assess the impacts of random DGs on distribution systems.Many risk indices have been established in [5][6][7][8].These papers aimed to establish the risk indices by using the product of probability and severity, but ignored many significant aspects such as the value of load curtailment.Many other risk indices have also been presented.In [9], three risk indices-including LOLP (loss of load probability), EENS (expected energy not supplied), and ECOST (expected customer interruption cost)-were used for reliability and price risk assessment.In [10,11], the indices of SAIFI (system average interruption frequency index), SAIDI

Distribution of Wind and Photovoltaic DGs
In distribution systems, the output of DGs contributes to whether the load can be supplied when malfunction occurs.Therefore, the output randomness of DGs should be imitated by mathematical formulas.In this paper, wind generators and photovoltaic generators are only considered as DGs.

Output Uncertainty of Wind Generators
Large amounts of research in the past decades have demonstrated that wind speed v is the main stochastic factor that determines the output power of wind generators.Generally, a Weibull distribution can be used to imitate the stochastic wind speed v [6].The probability density function (PDF) for wind speed v is described as (1): where k and c are, respectively, the shape parameter and the scale parameter of the wind speed distribution.According to the historical data of wind speed v, these two indices can be estimated.
Based on the PDF of wind speed v, the output of a wind turbine generator can be acquired [28]: where P w is the active output power, P N is the rated output power of wind turbine generator, v N is rated wind speed, v in is cut-in wind speed, and v out is cut-out wind speed.

Output Uncertainty of Photovoltaic Generators
In most occasions, illumination intensity is thought to be the major factor that affects the active output of photovoltaic generators.Because of cloud cover and other insolation-reducing phenomena, the illumination intensity I can also be represented by a random variable.In general, illumination intensity I approximately follows a beta distribution with shape indexes α and β [29]: where I max is the maximum intensity during a certain interval.The two shape indexes α and β can be evaluated by the mean value and the variance of illumination intensity.
Many studies have showed that the active output of photovoltaic generators could be described as [30]: where A is photoelectric array area and η is the photoelectric transformation efficiency.Combining Equations ( 3) and ( 4), the PDF of photovoltaic generators' output power P v can be acquired: where P max v is the maximum value of P v , which can be calculated by According to (5), P v also follows beta distribution with shape indexes α and β.

Risk Assessment Indices
It should be noted that risk indices such as EENS, PLC, EFLC, and SI in [13][14][15] are more considered with the penetration of DGs.Therefore, in this paper, these four risk indices were used to reflect the system risk level of distribution systems.The calculation methods of these four risk indices are illustrated by the following [13]: (1) EENS (unit: MWh/y) represents the energy that is not supplied with the penetration of DGs in the distribution system, which can be computed by (6): where i is the load level, N L is the set for load levels, T i is the duration of i, and Q i is the set for system state s at load level i. p T (s) and C 0 (s) are, respectively, the occurrence probability and total load curtailment of system state s.
(2) PLC can be calculated by: (3) EFLC (unit: times/y) can be calculated by: where T is total duration of T i and m(s) is the set of component j at system state s.For component j, θ j represents transfer rate of component j, which is equal to repair rate µ j in this paper.( 4) SI (unit: min/y) represents the equivalent duration under an entire system outage during peak conditions, which can be calculated: where L is the value of peak load in 1 year.In these four risk indices, SI is closely related to EENS.In order to calculate the index of EENS, the total load curtailment at state s C 0 (s) should be calculated beforehand; this can be calculated by optimal power flow (OPF), introduced in Section 3.3.For all of these four risk indices, the occurrence probability of system state s p T (s) is needed.Therefore, an enumeration method is applied for system state selection, which is introduced in Section 4.2.

Improved Point Estimate Methods
Point estimate method (PEM), which was firstly proposed by Hong in 1998, focuses on the statistical information that is provided by the first few central moments of the input random variables [25].For each variable, K points are used, and K is a parameter named concentrations and depends on the Hong's method.These points and the function-which relates input and output variables-are used to obtain the information about the uncertainty associated with problem output random variables.
This paper relies on a 2m and 2m + 1 type scheme of PEM for computing the four risk indices described in Section 3.1, which gives a good trade-off between the solution accuracy and the computational efforts.m denotes the number of random variables, which are wind speed v or illumination intensity I.The 2m type scheme of PEM uses the first two concentrations for each input random variable (i.e., K = 2).For the 2m + 1 type scheme of PEM, K = 3.
Generally speaking, PEM is used to calculate the estimation value (Z j ) of F(X) (Z = F(X) = F(X 1 , X 2 , . . ., X m )) when the stochastic numerical characteristics of variables X 1 , X 2 , . . ., X m are known [26,27].The point estimation principle can be depicted as: where j is the power of F(X), m is the number of DGs, µ Xi and ξ i,k are, respectively, the expectation and standard location of variable X i , which represents illumination intensity I or wind speed v in this paper.ω l,k represents the weight of X i at x l,k , which can be calculated by (11): where σ Xi represents the variance of X i .Standard location ξ i,k and weight ω l,k could be determined by (12): In traditional PEM, K ∑ k=1 ω i,k is set as 1 m , which means the impacts of X i on Z are the same.This assumption can simplify the calculation but may cause a big computational error as there are different impacts of X i on Z.In the improved PEM, K ∑ k=1 ω i,k is set as α i , which can be determined by the analytic hierarchy process (AHP) described in [31] or the contribution of X i on Z.This equation appears to be more reasonable, as it considers the different impacts of X i on Z.
When K = 2 and ξ i,3 = 0, the equation of ( 12) can be solved: 2 ) The calculation results based on (13) is named 2m improved PEM.When K = 3 and ξ i,3 = 0, 2m + 1 improved PEM can be constructed, and the solving results of Equation ( 12) are: The AHP is a rough estimation of weights and may be complex in operation.Therefore, the contribution of X i on Z which can be approximated by the capacity proportion of DG i was used to calculate the value of α i .
According to Equations ( 10)-( 14), the j-th raw moment of the output random vector Z can be acquired.Then, the mean and the standard deviation of Z can be estimated by Equation (15): For each risk assessment index, the mean and the standard deviation can be respectively acquired by Equation (15).For a variable with a certain mathematical distribution, the mean is the most probable value of the variable [26].Therefore, the mean of Z, µ Z , is used as the estimation value of each risk index in this paper.

Optimal Power Flow Algorithm in Distribution Systems
In distribution systems, load curtailment and generation dispatch are extensively applied to enable the system to go from urgent state to normal state.In order to obtain the system security indices in distribution systems, an optimal power flow (OPF) algorithm was presented in this paper to assess the total load curtailment C 0 (s) in state s.Equations ( 16) and (17) were used in the optimal power flow algorithm to simulate the curtailment and dispatch work and ensure the security of distribution system [32]: where l is the number of load buses in distribution systems.∂z ∂x (z = I ij , P s ; x = P Gk , P Lh ) are sensitivity coefficients.N b is the number of buses; ∆P Gk and ∆P Lh are respectively the power variations of generator k and load h; C GK and C Lh are respectively the control cost of generator k and load h; and I ij is the current through overload line ij.
According to the results of this optimal power flow algorithm in distribution systems, the total load curtailment in state s, C 0 (s), can be calculated by Equation ( 18): where ∆P d is the value of load curtailment caused by direct structure change of the distribution system.In line with the values of C 0 (s) and p T (s), the risk assessment indices in Equations ( 6)-( 9) can be calculated.

Structure for Risk Assessment
To assess the four risk indices in distribution systems, PEM was used and the uncertainties of distribution generators were sufficiently considered.The flowchart of the improved PEM-based method for risk assessment of distribution systems is shown in Figure 1.
In Step A, active outputs of photovoltaic generators and wind generators are calculated according to the stochastic wind speed v and illumination intensity I.In step B, random operation state s of the distribution system was firstly selected.Then, optimal power flow algorithm in state s was applied to compute the value of load curtailment based on the improved point estimate method.Then, the four risk indices listed in Section 3.1 were computed until all the random operation states s are considered.
where l is the number of load buses in distribution systems.According to the results of this optimal power flow algorithm in distribution systems, the total load curtailment in state s, C0(s), can be calculated by Equation ( 18): where ΔPd is the value of load curtailment caused by direct structure change of the distribution system.In line with the values of C0(s) and pT(s), the risk assessment indices in Equations ( 6)-( 9) can be calculated.

Structure for Risk Assessment
To assess the four risk indices in distribution systems, PEM was used and the uncertainties of distribution generators were sufficiently considered.The flowchart of the improved PEM-based method for risk assessment of distribution systems is shown in Figure 1.
In Step A, active outputs of photovoltaic generators and wind generators are calculated according to the stochastic wind speed v and illumination intensity I.In step B, random operation state s of the distribution system was firstly selected.Then, optimal power flow algorithm in state s was applied to compute the value of load curtailment based on the improved point estimate method.Then, the four risk indices listed in Section 3.1 were computed until all the random operation states s are considered.

Procedure for Calculating Risk Indices
In [13], according to DGs' probability model, the indices of EENS, PLC, EFLC, and SI were calculated by a proposed hierarchical risk assessment method, which is the combination of a Monte Carlo simulation and an enumeration method.However, this proposed hierarchical method is relatively complex in application.In [6], 2m + 1 PEM was utilized to assess the overload risk of a transmission line.This 2m + 1 PEM has lower accuracy compared to a Monte Carlo simulation.However, the computational burden of 2m + 1 PEM is greatly reduced.Consequently, the improved PEM is applied to compute the four risk indices in this paper.For comparison, the calculation results of the four risk indices based on the hierarchical risk assessment method in [13] were used as the exact values, and the error in accuracy of the improved PEM-based method for risk assessment was analyzed.
For convenience, symbol G is used to represent any parameter of the four risk indices including EENS, PLC, EFLC, and SI.Namely, to calculate the parameter G is to calculate these four risk indices.The designed procedure for computing the symbol G based on PEM is summarized as Figure 2.

Procedure for Calculating Risk Indices
In [13], according to DGs' probability model, the indices of EENS, PLC, EFLC, and SI were calculated by a proposed hierarchical risk assessment method, which is the combination of a Monte Carlo simulation and an enumeration method.However, this proposed hierarchical method is relatively complex in application.In [6], 2m + 1 PEM was utilized to assess the overload risk of a transmission line.This 2m + 1 PEM has lower accuracy compared to a Monte Carlo simulation.However, the computational burden of 2m + 1 PEM is greatly reduced.Consequently, the improved PEM is applied to compute the four risk indices in this paper.For comparison, the calculation results of the four risk indices based on the hierarchical risk assessment method in [13] were used as the exact values, and the error in accuracy of the improved PEM-based method for risk assessment was analyzed.
For convenience, symbol G is used to represent any parameter of the four risk indices including EENS, PLC, EFLC, and SI.Namely, to calculate the parameter G is to calculate these four risk indices.The designed procedure for computing the symbol G based on PEM is summarized as Figure 2.
, , ( , , , , )  In Figure 2, 2m and 2m + 1 improved-PEM were used, respectively.The output power of DG i is calculated by Equation ( 2) or (4).Further, OPF was applied to calculate the value of C 0 (s), which can be referred from Equations ( 16) and (17).The detailed procedure of the improved PEM-based method for computing risk indices G is summarized as follows.
Step (1): Enumeration method is used to list the states of the distribution system.The occurrence probability p T (s) in system state s can be calculated by Equation ( 19): where N f is the total number of failure components, N n is the total number of normal components.For component i, µ i is the repair rate, λ i is the outage rate, and these two parameters can be obtained by statistics.
Step (2): Initialize and generate the constructed points of wind speed v as well as illumination intensity I. Initialize the initial value of load curtailment calculation E(Z) = 0 and the initial value of G = 0.
Step (3): Construct the parameters of improved point estimate method, which contains the following substeps: 1) Select the input variable X i (i = 1, 2, . . ., m).For wind generating unit i, X i represents the stochastic wind speed v; for photovoltaic generating units i, X i represents the stochastic illumination intensity I. m is the total number of wind generators and photovoltaic generators.Then, the output power of DG i is calculated by Equation ( 2) or (4).
3) Based on the output of DG i, optimal power flow in state s is applied to compute the value of C 0 (s) according to Equations ( 15) and (16).
5) Repeat sub-steps 2-4 until variable i reaches the total number m, and the calculation of G in state s is completed.
Step (5): Output the values of the four risk indices EENS, PLC, EFLC, and SI.

Case Studies
In order to verify the rationality of the proposed PEM-based method for risk assessment of a distribution system, IEEE 30-bus system with DGs is applied, as shown in Figure 3.For all nodes, the upper voltage bound is 1.06 p.u. and lower voltage bound is 0.94 p.u. in this paper [6].All the cases are achieved by MATLAB in an Advanced Micro Devices 64 Dual Core 3.3 GHz PC.For photovoltaic farms, shape indexes α = 15.42,β = 4.3, A = 2.17 m, η = 13.53%.For wind farms, rated wind speed v N = 15 m/s, cut-in wind speed v in = 4m/s, cut-out wind speed v out = 20 m/s, shape index k = 6.23, and scale index c = 10.43.

Calculation of Risk Indices Based on PEM
In this modified IEEE 30-bus system, two wind farms and two photovoltaic farms are connected to nodes 7, 14, 24, and 29.The capacity proportion of these four farms is 20%:30%:20%:30%.However, the ratio of DGs' total generation capacity to the total load Kd varies from 0 to 1.With the known capacity proportion of DGs, the value of αi can be obtained: α1 = 0.2, α2 = 0.3, α3 = 0.2, α4 = 0.3.
According to the designed procedure for computing risk indices in Figure 2 and the parameters of DGs, the calculation results of EENS, PLC, EFLC, and SI, based on 2m and 2m + 1 improved PEM, are listed in Tables 1 and 2, respectively.As can be seen from Tables 1 and 2, with the increase of total generation capacity Kd, all of these risk indices decrease gradually.When the generation capacity Kd = 0, which means that DGs are not permeated, the value of EENS is about 3156.7 MWh/year.However, EENS decreases by nearly a half and drops to about 1582.5 MWh/year when the generation capacity Kd rises to 1. Analogously,

Calculation of Risk Indices Based on PEM
In this modified IEEE 30-bus system, two wind farms and two photovoltaic farms are connected to nodes 7, 14, 24, and 29.The capacity proportion of these four farms is 20%:30%:20%:30%.However, the ratio of DGs' total generation capacity to the total load K d varies from 0 to 1.With the known capacity proportion of DGs, the value of α i can be obtained: According to the designed procedure for computing risk indices in Figure 2 and the parameters of DGs, the calculation results of EENS, PLC, EFLC, and SI, based on 2m and 2m + 1 improved PEM, are listed in Tables 1 and 2, respectively.As can be seen from Tables 1 and 2, with the increase of total generation capacity K d , all of these risk indices decrease gradually.When the generation capacity K d = 0, which means that DGs are not permeated, the value of EENS is about 3156.7 MWh/year.However, EENS decreases by nearly a half and drops to about 1582.5 MWh/year when the generation capacity K d rises to 1. Analogously, the value of PLC decreases from about 0.0062293 to 0.0045285 when generation capacity K d rises from 0 to 1.
The variation tendency of SI is consistent with the variation tendency of EENS since the index of SI is calculated by EENS, which is shown in Section 3.1.Similarly, the variation tendency of EFLC is consistent with the variation tendency of PLC.In addition, the calculation results of risk indices based on 2m improved PEM are very close with the calculation results of risk indices based on 2m + 1 improved PEM.
From Tables 1 and 2, EENS and SI decrease smoothly with the increment of generation capacity K d .In addition, the slope of EENS and SI decreases with the growth of generation capacity K d .Therefore, the decreasing trend of EENS and SI is not as tangible as before, decreasingly when generation capacity K d increases.
As shown in Tables 1 and 2, the decrease of EENS, PLC, ELIC, and SI slows when generation capacity K d increases to about 0.6.Namely, a large value of K d that promotes the decrease of risk indices of distribution systems is not significantly remarkable.However, the consumption of wind and photovoltaic energy becomes an increasing problem with the increasing generation capacity K d .It is with reluctance that too much wind and photovoltaic power is abandoned in distribution systems.Therefore, it is not necessary to increase the value of K d to a certain large degree.

Deviation and Computational Cost Comparison
In order to illustrate the effectiveness of this improved PEM-based method for risk assessment in distribution systems, the calculation results of EENS, PLC, EFLC, and SI are shown in Table 3. Comparing the calculation results of risk indices in Tables 1 and 2 with Table 3, it can be seen that they are very close in value.However, the calculation results of these improved PEM-based methods are all less than that of the hierarchical risk assessment method since PEM is an order approximation, which ignores the higher-order terms.Besides, the calculation results of 2m + 1 improved PEM are closer to the exact value than the calculation results of 2m improved PEM, as 2m improved PEM is a third-order approximation but 2m + 1 improved PEM is a fourth-order approximation, in fact.
According to the analysis in Section 3.1, the variation tendency of SI is consistent with the variation tendency of EENS, and the variation tendency of EFLC is consistent with the variation tendency of PLC.Therefore, the deviation comparison of EENS and PLC are only needed to be analyzed, which are shown in Figures 4 and 5.
In Figures 4 and 5, errors δ 1 and δ 2 are respectively calculated by comparing the results of EENS and PLC in Tables 1 and 2 with Table 3.As can be seen from Figures 4 and 5, the maximum value of computational error δ 1 is about 0.99% and the maximum value of computational error δ 2 is about 0.37%.The computational errors δ 1 and δ 2 are all less than 1%, which greatly verifies the effectiveness of this improved PEM-based method for risk assessment of distribution systems.Also, the computational error δ 2 of PLC is less than the computational error δ 1 of EENS.This is because optimal power flow algorithm (based on Equations ( 15) and ( 16)) is used in the calculation procedure for EENS.
Compared with hierarchical risk assessment method-which is a combination of the Monte Carlo method and enumeration method-in fact, the computational burden of this PEM-based method is greatly reduced.When generation capacity K d = 0.5, the computational costs of risk indices are as listed in Table 4.
Sustainability 2017, 9, 491 11 of 14 method is greatly reduced.When generation capacity Kd = 0.5, the computational costs of risk indices are as listed in Table 4.As can be seen from Table 4, the computational costs of this improved PEM-based method for risk assessment in distribution systems are much less than the computational costs of the hierarchical method.For EENS, the computational costs of 2m improved PEM and 2m + 1 improved PEM are 25.845 s and 27.742 s, respectively, but the computational cost of the hierarchical method is 67.293 s.This is because the Monte Carlo method, with many computational hurdles, is used in the hierarchical method.On the contrary, only a few numerical operations are needed in the improved PEM.The calculation of risk indices has to be evaluated only K times for each input variable, which is stochastic wind speed v or stochastic illumination intensity I at the K points.
In addition, the computational costs of 2m + 1 improved PEM-based method are almost the same compared to 2m improved PEM-based method, as the only difference between these two method is greatly reduced.When generation capacity Kd = 0.5, the computational costs of risk indices are as listed in Table 4.As can be seen from Table 4, the computational costs of this improved PEM-based method for risk assessment in distribution systems are much less than the computational costs of the hierarchical method.For EENS, the computational costs of 2m improved PEM and 2m + 1 improved PEM are 25.845 s and 27.742 s, respectively, but the computational cost of the hierarchical method is 67.293 s.This is because the Monte Carlo method, with many computational hurdles, is used in the hierarchical method.On the contrary, only a few numerical operations are needed in the improved PEM.The calculation of risk indices has to be evaluated only K times for each input variable, which is stochastic wind speed v or stochastic illumination intensity I at the K points.
In addition, the computational costs of 2m + 1 improved PEM-based method are almost the same compared to 2m improved PEM-based method, as the only difference between these two  As can be seen from Table 4, the computational costs of this improved PEM-based method for risk assessment in distribution systems are much less than the computational costs of the hierarchical method.For EENS, the computational costs of 2m improved PEM and 2m + 1 improved PEM are 25.845 s and 27.742 s, respectively, but the computational cost of the hierarchical method is 67.293 s.This is because the Monte Carlo method, with many computational hurdles, is used in the hierarchical method.On the contrary, only a few numerical operations are needed in the improved PEM.The calculation of risk indices has to be evaluated only K times for each input variable, which is stochastic wind speed v or stochastic illumination intensity I at the K points.
In addition, the computational costs of 2m + 1 improved PEM-based method are almost the same compared to 2m improved PEM-based method, as the only difference between these two methods are a few numerical operations.These simulations have demonstrated that the improved PEM-based method has high accuracy and highly reduced computational costs, therefore, it is extremely applicable for risk assessment of distribution systems.

Influence of DGs on System Risk
Many researches have shown that many other aspects such as type, location, dispersion, and capacity proportion of DGs have great influences on the system risk indices.In this part, the influences of these aspects are analyzed.
In Section 5.1, two wind farms and two photovoltaic farms are connected to nodes 7, 14, 24, and 29.Suppose that four wind generating units are respectively connected to nodes 7, 14, 24, and 29 and their capacity proportion is still 20%:30%:20%:30%.This case is named case 1 in this paper.According to the developed 2m + 1 improved PEM-based method for risk assessment in Figure 2, the four system risk indices of EENS, PLC, EFLC, and SI can be acquired when generation capacity K d = 0.2, which can be seen in Table 5.As can be seen from Table 5, the index of EENS decreases from 2385.4 MWh/y to 2131.7 MWh/y and the values of other indices are also cut down distinctly.Therefore, for DGs' type, wind farms' power support is better than that of photovoltaic farms.These results reflect that the randomness of the output power of photovoltaic farms is even bigger than the active output power of wind farms.
To illustrate the influence of DGs' location, suppose that two wind farms and two photovoltaic farms are respectively connected to nodes 3, 14, 19, and 30 in this part, and their capacity proportion is still 20%:30%:20%:30%.This case is named case 2 in this paper.Likewise, the calculation results of system risk indices based on 2m + 1 improved PEM, when generation capacity K d = 0.2, are outlined in Table 6.Comparing the calculation results of Table 6 with Table 2, it is observed that the index of EENS decreases from 2385.4 MWh/y to 2193.6 MWh/y and the values of other indices are also cut down distinctly.This comparison shows that the locations of nodes 3, 14, 19, and 30 are more suitable for optimal siting of DGs than nodes 7, 14, 24, and 29.However, these locations may be not the most optimal siting for DGs, and how to ascertain the optimal siting of DGs is worth further study.
To illustrate the influence of DG dispersion, suppose that three wind farms and three photovoltaic farms are respectively connected to nodes 3, 7, 14, 24, 29, and 30 in this part, and their capacity proportion is 15%:15%:20%:15%:15%:20%.This case is named case 3 in this paper and the calculation results of system risk indices are outlined in Table 7.In case 3, the total proportion of wind farms and photovoltaic farms is still 50%:50% but the dispersion degree of DGs in case 3 is enhanced.The results in Tables 2 and 7 show that system risk indices decrease along with dispersion degree increment.Therefore, when installed, DGs should be as highly dispersed as possible for the optimal siting of DGs for application.
The influence of DGs' capacity proportion reflects the influence of DGs' type on system risk indices.As wind farms' power support is better than that of photovoltaic farms, the values of system risk indices are cut down when the wind generating units' capacity proportion is increased.
These test results illustrate that appropriate DGs' siting and sizing have great influences on the system risk indices.For the sizing of DGs, system risk indices increase along with the DG's capacity, however, the corresponding investment costs will increase.For the appropriate siting of DGs, different locations may suit different types of DGs, and DGs' dispersion can decrease the value of system risk indices.

Conclusions
DGs connected to a distribution system could cause some critical risks to distribution systems from security and economy aspects.In order to reasonably assess the risks of distribution systems with penetrating DGs, four risk indices-EENS, PLC, EFLC, and SI-were used in this paper to reflect the system risk level in distribution systems.The output uncertainties of DGs were depicted by some certain mathematical distributions.Then, an improved PEM-based method was proposed to calculate these four system risk indices.In this improved PEM-based method, enumeration method was used to list the states of distribution systems, and an improved PEM was presented to deal with the uncertainties of DGs.
Test results of case studies have shown that this proposed PEM-based method is highly effective for assessing the risk of a distribution system with DGs, as it has a high computational accuracy and highly reduced computational costs compared with other risk assessment methods.Simulation results also demonstrate that total generation capacity, type, location, dispersion, and capacity proportion of DGs have great influences on the system risk indices.The determining method of appropriate DGs' siting and sizing, considering the operation risk of a distribution system with DGs, is worthy of further research.

Figure 1 .
Figure 1.Diagram of the improved point estimate method (PEM)-based method for risk assessment of distribution lines in a distribution system.

Figure 1 .
Figure 1.Diagram of the improved point estimate method (PEM)-based method for risk assessment of distribution lines in a distribution system.

Figure 2 .
Figure 2. Flowchart of computing procedure for risk indices.

Figure 2 .
Figure 2. Flowchart of computing procedure for risk indices.

Figure 4 .
Figure 4.The computational error δ1 of EENS in pace with Kd.

Figure 5 .
Figure 5.The computational error δ2 of PLC in pace with Kd.

Figure 4 .
Figure 4.The computational error δ 1 of EENS in pace with K d .

Figure 4 .
Figure 4.The computational error δ1 of EENS in pace with Kd.

Figure 5 .
Figure 5.The computational error δ2 of PLC in pace with Kd.

Figure 5 .
Figure 5.The computational error δ 2 of PLC in pace with K d .
Nb is the number of buses; ΔPGk and ΔPLh are respectively the power variations of generator k and load h; CGK and CLh are respectively the control cost of generator k and load h; and Iij is the current through overload line ij.

Table 1 .
Calculation results of risk indices based on 2m improved PEM.
Kd: generation capacity; EENS: expected energy not supplied; PLC: probability of load curtailment; EFLC: expected frequency of load curtailment; SI: severity index.

Table 2 .
Calculation results of risk indices based on 2m + 1 improved PEM.

Table 1 .
Calculation results of risk indices based on 2m improved PEM.
d : generation capacity; EENS: expected energy not supplied; PLC: probability of load curtailment; EFLC: expected frequency of load curtailment; SI: severity index.

Table 2 .
Calculation results of risk indices based on 2m + 1 improved PEM.

Table 3 .
Calculation results of risk indices based on hierarchical risk assessment method.

Table 4 .
Calculation results of risk indices based on hierarchical risk assessment method.

Table 4 .
Calculation results of risk indices based on hierarchical risk assessment method.

Table 4 .
Calculation results of risk indices based on hierarchical risk assessment method.

Table 5 .
Calculation results of risk indices in case 1.

Table 6 .
Calculation results of risk indices in case 2.

Table 7 .
Calculation results of risk indices in case 3.