# Reliability Assessment of Active Distribution System Using Monte Carlo Simulation Method

Authors: Shaoyun Ge; Li Xu; Hong Liu; Mingxin Zhao
Journal: Journal of Applied Mathematics
 (2014)

---

## Abstract

In this paper we have treated the reliability assessment problem of low and high DG penetration level of active distribution system using the Monte Carlo simulation method. The problem is formulated as a two-case program, the program of low penetration simulation and the program of high penetration simulation. The load shedding strategy and the simulation process were introduced in detail during each FMEA process. Results indicate that the integration of DG can improve the reliability of the system if the system was operated actively.


---

## Body

## 1. Introduction

As an effective complement to centralized power generation, distributed generation (DG) technologies are becoming mature. With flexible operation mode and environment-friendly features, more and more distributed generations were integrated into the distribution network, which brings a series of impacts on the distribution system and challenges the original reliability assessment method [1]. As two basic methods of reliability assessment for the distribution system with DG, the studies of analytical method [2–8] and Monte Carlo method [9–15] have made some progress. These studies indicate that the adoption of DG was actually equal to increasing the number and capacity of the backups of distribution system, thus improving the reliability of the system.The traditional distribution system reliability evaluation is based on failure mode and effect analysis method, which start with failure of a component followed by analysis of the failure effect of users or load points in the system and then the calculation of reliability indices associated with the user, such as SAIDI and SAIFI. However, in the aforementioned case of the high penetration level of DG, even if all the components in the system are in normal state, DG’s intermittent output may still cause system power shortage. Therefore, just considering the fault condition cannot take the power shortage risk of system reliability into account.When DG was integrated into distribution system, both low penetration level case and high penetration level case should be considered. In Section2, the boundary conditions of the study were confined. In Section 3, the DGs mathematical modeling was introduced, and the Monte Carlo simulation methods designed for low penetration case and high penetration case were proposed in detail. In Section 4, a study case was carried out to exemplify the methods mentioned above. At last, the conclusion was drawn in Section 5.


## 2. System Configuration

While such technical issues are fundamental for a real-world system, in the analysis presented in this paper, emphasis is placed on the energy transactions between system components on relatively long time scales (hours and beyond), and, hence, all dynamic phenomena and associated components have been ignored.Penetration of DG is defined as follows only from a view of reliability. When all nonsource components (bus, feeders, switches, distribution transformers, etc.) are in the nonfault state and no load in the distribution system would shed due to the lack of DGs’ output, the distribution system is defined to have a low penetration; otherwise it has a high penetration.Reliability indicators of the active distribution network are the same as traditional distribution network in a low penetration of DGs. Reliability indicators of the active distribution network include reliability indicators of the generation network in addition to a high penetration, such as LOLE (loss of load expectation) and LOEE (loss of expected energy).


## 3. Modeling and Simulation

### 3.1. Simulation Principles

Modeling of the system accounts for the energy transactions between the various components, while all electric transients are ignored, so that long-term simulations (e.g., annual) are possible. The simulation code developed includes appropriate component models and reproduces faithfully the operating policy presented in this section. A one-hour simulation time step is used, and scheduling decisions are made at the beginning of each simulation hour.


### 3.2. Wind Turbines (WTG)

Real-time wind speeds are generated in each period of a year with autoregressive moving average model in [16], based on which real-time output of wind turbines is calculated by output curve of wind turbines [17].


### 3.3. Photovoltaic (PV)

Real-time output of PV panels is calculated by real-time solar radiation and output curve [18] like output model of wind turbines.


### 3.4. Battery Energy Storage Systems (BESS)

In this paper we choose the lead-acid battery as the BESS. A variant of the kinetic battery model (KiBaM) [19] is used to simulate the lead-acid BESS. KiBaM considers the battery as a two-well system, where the first well represents directly available capacity, while the second one contains chemically bound energy, which may become available at a limited rate. This model takes into account capacity reduction with increased charge or discharge currents, as well as the recovery effect.


### 3.5. Simulation of Low Penetration

#### 3.5.1. Sequential Sampling for Nonsource Component

All nonsource component was simulated by the two-state Markov model shown as in Figure1, the fault transfer rate was λ, and the recovery rate was μ, so the fault possibility was Pd=λ/(λ+μ), the regular state possibility was Pu=μ/(λ+μ), time between failures (TTF) and fault repair time (TTR) of components are exponentially distributed, and the Monte Carlo simulation method for state duration sampling was adopted for its sequential sampling [20].Figure 1
The two-state Markov model.


#### 3.5.2. Sampling for DG

When not considering the limitation of upper system capacity, only when the nonsource components in system is in fault, the outage of DG and BESS can 6 exacerbate the power off risk of the user. Compared with normal state, the nonsource components fault duration is very short, and DG and BESS states can be considered unchanged during the fault, so their state can be sampled nonsequentially.Single wind turbine is still using two-state model of Figure1. For wind turbines set consisting of m wind turbines unit, m random number between [0,1] is first generated on behalf of its sample state, and Ri>Pd represents that the unit is working.PV array is consisting of a set of photovoltaic panels, to reflect a portion of PV array’s fault, on the basis of two-state model; the derating state is added, and its probability isPe. Sampling method was shown in Figure 2; first, the probability of all the states was placed cumulatively on number of axes between [0,1], and the random number R represents which state the PV array is in. For example, in Figure 3, the photovoltaic array is in derating state.Figure 2
Sampling process for the three-state element.Figure 3
A distribution system with WTG, PV, and BESS.The sample of the BESS state consists of two parts: one is sample of the running state, and the other is sample of the state of charge (SOC). Sampling methods of running state are the same with PV array, and derating state represents part of BESS’s failure. When sampling the SOC, SOC was converted from a random number in accordance with the SOC changing law under integration mode.


#### 3.5.3. Failure Mode and Effect Analysis (FMEA) Process

FMEA analysis is the basis of reliability assessment. The DG’s integration will have an impact on the traditional FMEA process. In this paper, feeder segment concept was used to analyze the FMEA process of distribution system with DG. Feeder segment is a set of equipment in which all equipment has a common entrance [21] and all the load points in the same segment have the same failure consequences. For an example of the distribution system as Figure 3, the system can be divided into 12 feeder segments (S1~S12).Assuming the upper system’s capacity is adequate, when the system is in a normal state, all load points can be supplied by upper system and DG at the same time. When DG is in fault, the circuit breaker of DG can act immediately to isolate DG, and no load will suffer power loss. And only when the nonsource components fail, the load point will face the risk of power loss. At this point, the active islanding can narrow outage range and reduce outage time. However, the realization of the active islanding needs the adjusting of traditional protection. Therefore, the protection logic in paper [22] was used in this paper: for any faulty component of a feeder, only its nearest surrounding breaker/isolator switch will act for fault isolation.Still use Figure3 system as example, when the components in S5 (feeders, distribution transformers) was in fault, the acting circuit breakers should be B3, B4, B5, B6, B8, B9, B10, B12, and B13. Then, the isolation switches D2, D3, and D4 will act to isolate S5. According to the differences in each feeder segment’s power restoration result, the system can be divided into six kinds of segment as follows.Fault Segment ({S5}). Outage time for all load points in the region is the fault repair time.Unaffected Segment ({S1}). Load point in the region will not be affected.Upstream Quarantine Segment ({S2},{S3}). Located at upstream of fault segment, before fault isolation, load point in this region cannot restore power. After fault isolation, S3 as an example, isolating switch D1 and circuit breaker B5 will close, and both the DG and the system will supply load in S3. All the load point’s outage time in this region is fault isolation time.Upstream Active Islanding Segment ({S4}). When a fault occurs, due to the operation of the circuit breaker (B6), the region can immediately run in island mode, the island time is fault isolation time, the load point outage situation was decided by the island's power balance status, and when necessary, load shedding will be carried out, which would be described later. After fault isolation was over, the area can be reaccessed into systems.Downstream Quarantine Segment ({S6},{S7},{S9}). Before fault isolation, the region cannot run in island mode, and all load points suffer outage. After fault isolation, the region can be linked with other downstream segments to form a big island. For example, after fault isolation, feeder segments S7~S12 can form a big island, and island operation time was the difference between fault repair time and fault isolation time.Downstream Active Islanding Segment ({S8},{S10~S12}). Same as the upstream active islanding segment, the island operation time is divided into two parts: one is the seamless operation time before fault isolation; the other is the big island operation time after fault isolation.In addition, if the fault component is switch, we can simply merge feeder segments which connected by switch into one feeder segment and then reuse the FMEA analysis method above.


#### 3.5.4. Load Shedding Strategy

Load shedding determination conditions are as follows: at any moment of battery discharge while islanding, the greatest power all the battery can release was no less than the net exchange power in island. Consider(1)∑i=1nPout,1i(t)ηd≥Pex(t),Pex(t)≥0,t∈[tst,tend].Pout,1i(t) is the maximum power which can be released by the battery pack i during time interval t and n is the number of battery packs; Pex(t)=L(t)+Sl(t)-DG(t), in which L(t) is the total load; DG(t) is the sum of all DG outputs; Sl(t) is the active power loss within island and, to simplify the flow calculation, Sl(t) can be set at about 5% of the current load [23].In order to improve simulation speed, this paper takes the following heuristic load shedding strategies. First, assuming that all load points are not shed, when the case does not meet formula (1), the smallest load was shed first; repeat this process until it satisfies formula (1). During load shedding process, the dynamic switching was not considered, so the shedding state was unchanged during islanding time. When formula (1) was satisfied, each battery power was released by the proportional ratio of Pout,1i(t) and power to be absorbed was by proportional ratio of max⁡(Pin,1i(t),-Pin,2i(t),-Pin,3i(t)). It should be noticed that, for the load points in downstream quarantine segment and downstream active islanding segment, if a load point was shed in both parts of the two island times, the frequency of power-off is 1, and the power outages time was the sum of two parts.


#### 3.5.5. Simulation Process

(1)
Input the network structure; use the method in Section3.5.3, automatically traversing all the nonfaulty component’s feeder segments classification by computer; establish FMEA table.
(2)
Set initial analog clock value for 0; at this time all system components are in normal state. For each nonsource system component, generate a random number; use the state transition model to generate TTF for each component.
(3)
Find the component with the smallest TTF (TTFmin⁡), and advance the analog clock to TTFmin⁡. If some components are of the same TTF, randomly select one of them.
(4)
Generate a new random number and convert it into TTR for this component. At the same time, generate power isolation and recovery time (RT).
(5)
By using the FMEA table generated in step 1, check the feeder segments under this component’s outage. Fault segment load point outage time was TTR; upstream quarantine segment load point outage time was RT; record the time duration of power outages and the insufficiency of power supply.
(6)
For each upstream active island segment, sample DG’s state and BESS’s SOC atTTFmin⁡; calculate real-time output of DG and load during [TTFmin⁡, TTFmin⁡ + RT]; use the method in Section 3.5.4 to shed load; record the load point outage frequency and duration.
(7)
Similar to step 6, in [TTFmin⁡, TTFmin⁡ + RT], simulate the downstream active island; in [TTFmin⁡ + RT, TTFmin⁡ + TTR], simulate the downstream big island and record the load point outage frequency and duration.
(8)
For the faulty component in step 3, regenerate a random number and convert it intoTTFN for a new running time. At this moment, the component’s TTF is updated to TTFmin⁡ + TTR + TTFN.
(9)
If the analog clock is less than predetermined simulation time, return to step 3; otherwise, statistically count the outage frequency, duration, and insufficiency of each load point, and then calculate the reliability indices.


### 3.6. Simulation of High Penetration

#### 3.6.1. Feature of High Penetration System

High penetration distribution network with DG will have both power generation and distribution properties, so-called active distribution network, and the key of its reliability assessment is how to properly combine power generation and distribution properties together. To solve this problem, this paper uses the sequential Monte Carlo sampling method based on the system state transition sampling to conduct the reliability assessment for distributed network in high DG penetration level, which can simultaneously calculate the reliability index of power generation and distribution.Compared to the commonly used state duration sampling method, system state transition sampling method is to sample the state of the entire system rather than of each specific component; there is no need to store the run-outage state cycling sequence of each component within the system, so it can greatly save data storage space and improve the simulation speed, which is more suitable for system with a large number of components, especially the systems including much power generation equipment.


#### 3.6.2. System State Transition Sampling Method

The system state transition sampling method was focused on the state of the whole system rather than that of one component in the system. Assume that the system is composed ofm components and consists of n states. The state space of the system is G={S1,…,Sn}. Suppose that the present system state is Sk and the state transition rate of each component in the system is λi(i=1,…,m). Since the state duration of each component follows an exponential distribution, the duration of the system state Sk also follows the exponential distribution. The duration of the system state Sk has the following probability density function:
(2)f(t)=∑i=1mλiexp⁡⁡(-∑i=1mλit).It can be seen that the state transition rate to the next state from the system stateSk is the sum of each component’s state transition rate in the system state Sk.State transition of any component in the system can lead to a system state. Consequently, system hasm possible states starting from state Sk. The probability that the system reaches one of these possible states is given by the following equation:
(3)Pj=λj∑i=1mλi.λj is the state transition rate of the jth component. It is obvious that
(4)∑j=1mPj=1.Two uniformly distributed random numbersR1 and R2 between [0,1] are generated to sample the duration of the system state Sk and the system state which is to be reached. The duration of the system state Sk is given by the following equation:
(5)Tk=-ln⁡R1∑i=1mλi.The next system state can be determined by the sampling shown in Figure4. The probabilities of m possible future states are successively placed in the interval [0,1]. If R2 falls into the segment corresponding to Pj, the transition of the jth component leads system to the next system state.Figure 4
Explanation of system state transition sampling.


#### 3.6.3. Analysis and Simulation of Generation State

According to different states of the components within the system, the state space of the distribution system with distributed generationG can be divided into two subsets G1 and G2. In the set G1, all nonsource components are in normal state, and in the set G2, a nonsource element is in fault state.In the traditional distribution system, when all nonsource components are in normal state, the system load can be supplied without considering the higher level power supply capacity. But the distribution system with DG is a multiple supply network, in which the load will be supplied by both the DG and the higher level power. Because DG can supply part of the load, the substation capacity will be reduced in order to save cost of grid construction. So it is necessary to consider the impact of the higher level power supply capacity. In this case, whether the total generation of the system including DG can supply all the load needs to be assessed. It is in fact a part of the generation system reliability evaluation.The reliability of generation system can be evaluated by simple simulation. The reliability indices are determined by superposing the chronological load curve on the generation capacity curve. Assume that the supply from higher level isPs(t) and the DG provides PDG(t) during the tth period. The sum of the real time load and the loss is PL(t) (assume that the loss is 5% of the load to simplify the power flow calculation [24]). It means that the system is suffering power shortage when Ps(t)+PDG(t)<PL(t) and the ENS can be calculated by ENS=PL(t)-Ps(t)-PDG(t). The indices of generation system such as LOLE and LOEE are calculated with sequential simulation during all periods. The method will be improved with the consideration of battery charging and discharging strategy. Two dispatch strategies of battery used in the paper are as follows [25].Load-Following Strategy. Under the load-following strategy, renewable power sources charge the BESS but the generators do not. The BESS discharge when DG could not provide enough supply to the load. The generators serve as the backup.Cycle Charge Strategy. When DG output was surplus, the BESS was charged by excess power; when DG output was insufficient, if the maximum output of DG and BESS still cannot meet all the loads’ need, the system will first supply the remaining load and then charge the BESS. In addition, to prevent the battery from being in a low SOC, once the battery starts charging, it needs to be charged to a specified SOCset.Figures5 and 6 show the simulation process in load-following strategy and in cycle charging strategy, respectively. Pdi(t) is the maximum output of the ith BESS during the tth period. Pci(t) is the maximum active power that the battery could absorb. Considering the difference in the capacity and SOC between batteries, the exchanging power of each battery is proportional to its Pdi(t) or Pci(t). For example, when the DG is charging all the batteries, the ith battery can get
(6)Pexi(t)=(PDG(t)-PL(t))×Pci(t)∑i=1mPci(t).Figure 5
Simulation process in load-following strategy.Figure 6
Simulation process in cycle charging strategy.However, the powerPini(t) that it really absorbs is min⁡{ηcPexi(t),Pci(t)}. ηc is charge efficiency.


#### 3.6.4. Simulation Process

The time sequential Monte Carlo simulation using system state transition sampling evaluates the reliability of distribution system containing distributed wind energy, solar energy, and BESS. The step for simulation is 1 hour. The procedure used consists of the following steps.(1)
The simulation process starts from the normal system state in which all generating units and transmission components are in the upstate at timet0=0.(2)
Calculate the duration of the current system stateTk by (4) with the transition rate λi of each component. Get λi of a three-state component by adding its two state transition rates. If the present system state is a state in which at least one nonpower component is in fault state, set the transition rate of other components to zero.(3)
If all the nonpower components are in the normal state in the current state, call the generation system simulation process, and recordLLDk and ENSk during Tk. If there is a nonpower component in fault state, call the distribution system simulation process shown in Section 3.5 and record the indices λki, γki, and ENSki.(4)
Sett0=t0+Tk. If the time is shorter than the year limit, sample the next system state by the method in Figure 1 and go to step 2. If t0 is longer than N, stop simulating and calculate the following indices.Loss of load expectation (h/year) is(7)LOLE=1N∑k=1PLLDk.
Loss of expected energy (MWh/year) is(8)LOEE=1N∑k=1PENSk.
System average interruption frequency index (times/node* year) is
(9)SAIFI=(1/N)∑k=1Q∑i=1nλkiCi∑i=1nCi.
System average interruption duration index (h/node* year) is
(10)SAIDI=(1/N)∑k=1Q∑i=1nγkiCi∑i=1nCi.
Expected energy not supplied (MWh/year) is(11)EENS=1N∑k=1Q∑i=1nENSki.
P is the simulation times of the generation system state, Q is the simulation times of the distribution system state, Ci is the amount of users of each load point, and n is the amount of load points.


## 3.1. Simulation Principles

Modeling of the system accounts for the energy transactions between the various components, while all electric transients are ignored, so that long-term simulations (e.g., annual) are possible. The simulation code developed includes appropriate component models and reproduces faithfully the operating policy presented in this section. A one-hour simulation time step is used, and scheduling decisions are made at the beginning of each simulation hour.


## 3.2. Wind Turbines (WTG)

Real-time wind speeds are generated in each period of a year with autoregressive moving average model in [16], based on which real-time output of wind turbines is calculated by output curve of wind turbines [17].


## 3.3. Photovoltaic (PV)

Real-time output of PV panels is calculated by real-time solar radiation and output curve [18] like output model of wind turbines.


## 3.4. Battery Energy Storage Systems (BESS)

In this paper we choose the lead-acid battery as the BESS. A variant of the kinetic battery model (KiBaM) [19] is used to simulate the lead-acid BESS. KiBaM considers the battery as a two-well system, where the first well represents directly available capacity, while the second one contains chemically bound energy, which may become available at a limited rate. This model takes into account capacity reduction with increased charge or discharge currents, as well as the recovery effect.


## 3.5. Simulation of Low Penetration

### 3.5.1. Sequential Sampling for Nonsource Component

All nonsource component was simulated by the two-state Markov model shown as in Figure1, the fault transfer rate was λ, and the recovery rate was μ, so the fault possibility was Pd=λ/(λ+μ), the regular state possibility was Pu=μ/(λ+μ), time between failures (TTF) and fault repair time (TTR) of components are exponentially distributed, and the Monte Carlo simulation method for state duration sampling was adopted for its sequential sampling [20].Figure 1
The two-state Markov model.


### 3.5.2. Sampling for DG

When not considering the limitation of upper system capacity, only when the nonsource components in system is in fault, the outage of DG and BESS can 6 exacerbate the power off risk of the user. Compared with normal state, the nonsource components fault duration is very short, and DG and BESS states can be considered unchanged during the fault, so their state can be sampled nonsequentially.Single wind turbine is still using two-state model of Figure1. For wind turbines set consisting of m wind turbines unit, m random number between [0,1] is first generated on behalf of its sample state, and Ri>Pd represents that the unit is working.PV array is consisting of a set of photovoltaic panels, to reflect a portion of PV array’s fault, on the basis of two-state model; the derating state is added, and its probability isPe. Sampling method was shown in Figure 2; first, the probability of all the states was placed cumulatively on number of axes between [0,1], and the random number R represents which state the PV array is in. For example, in Figure 3, the photovoltaic array is in derating state.Figure 2
Sampling process for the three-state element.Figure 3
A distribution system with WTG, PV, and BESS.The sample of the BESS state consists of two parts: one is sample of the running state, and the other is sample of the state of charge (SOC). Sampling methods of running state are the same with PV array, and derating state represents part of BESS’s failure. When sampling the SOC, SOC was converted from a random number in accordance with the SOC changing law under integration mode.


### 3.5.3. Failure Mode and Effect Analysis (FMEA) Process

FMEA analysis is the basis of reliability assessment. The DG’s integration will have an impact on the traditional FMEA process. In this paper, feeder segment concept was used to analyze the FMEA process of distribution system with DG. Feeder segment is a set of equipment in which all equipment has a common entrance [21] and all the load points in the same segment have the same failure consequences. For an example of the distribution system as Figure 3, the system can be divided into 12 feeder segments (S1~S12).Assuming the upper system’s capacity is adequate, when the system is in a normal state, all load points can be supplied by upper system and DG at the same time. When DG is in fault, the circuit breaker of DG can act immediately to isolate DG, and no load will suffer power loss. And only when the nonsource components fail, the load point will face the risk of power loss. At this point, the active islanding can narrow outage range and reduce outage time. However, the realization of the active islanding needs the adjusting of traditional protection. Therefore, the protection logic in paper [22] was used in this paper: for any faulty component of a feeder, only its nearest surrounding breaker/isolator switch will act for fault isolation.Still use Figure3 system as example, when the components in S5 (feeders, distribution transformers) was in fault, the acting circuit breakers should be B3, B4, B5, B6, B8, B9, B10, B12, and B13. Then, the isolation switches D2, D3, and D4 will act to isolate S5. According to the differences in each feeder segment’s power restoration result, the system can be divided into six kinds of segment as follows.Fault Segment ({S5}). Outage time for all load points in the region is the fault repair time.Unaffected Segment ({S1}). Load point in the region will not be affected.Upstream Quarantine Segment ({S2},{S3}). Located at upstream of fault segment, before fault isolation, load point in this region cannot restore power. After fault isolation, S3 as an example, isolating switch D1 and circuit breaker B5 will close, and both the DG and the system will supply load in S3. All the load point’s outage time in this region is fault isolation time.Upstream Active Islanding Segment ({S4}). When a fault occurs, due to the operation of the circuit breaker (B6), the region can immediately run in island mode, the island time is fault isolation time, the load point outage situation was decided by the island's power balance status, and when necessary, load shedding will be carried out, which would be described later. After fault isolation was over, the area can be reaccessed into systems.Downstream Quarantine Segment ({S6},{S7},{S9}). Before fault isolation, the region cannot run in island mode, and all load points suffer outage. After fault isolation, the region can be linked with other downstream segments to form a big island. For example, after fault isolation, feeder segments S7~S12 can form a big island, and island operation time was the difference between fault repair time and fault isolation time.Downstream Active Islanding Segment ({S8},{S10~S12}). Same as the upstream active islanding segment, the island operation time is divided into two parts: one is the seamless operation time before fault isolation; the other is the big island operation time after fault isolation.In addition, if the fault component is switch, we can simply merge feeder segments which connected by switch into one feeder segment and then reuse the FMEA analysis method above.


### 3.5.4. Load Shedding Strategy

Load shedding determination conditions are as follows: at any moment of battery discharge while islanding, the greatest power all the battery can release was no less than the net exchange power in island. Consider(1)∑i=1nPout,1i(t)ηd≥Pex(t),Pex(t)≥0,t∈[tst,tend].Pout,1i(t) is the maximum power which can be released by the battery pack i during time interval t and n is the number of battery packs; Pex(t)=L(t)+Sl(t)-DG(t), in which L(t) is the total load; DG(t) is the sum of all DG outputs; Sl(t) is the active power loss within island and, to simplify the flow calculation, Sl(t) can be set at about 5% of the current load [23].In order to improve simulation speed, this paper takes the following heuristic load shedding strategies. First, assuming that all load points are not shed, when the case does not meet formula (1), the smallest load was shed first; repeat this process until it satisfies formula (1). During load shedding process, the dynamic switching was not considered, so the shedding state was unchanged during islanding time. When formula (1) was satisfied, each battery power was released by the proportional ratio of Pout,1i(t) and power to be absorbed was by proportional ratio of max⁡(Pin,1i(t),-Pin,2i(t),-Pin,3i(t)). It should be noticed that, for the load points in downstream quarantine segment and downstream active islanding segment, if a load point was shed in both parts of the two island times, the frequency of power-off is 1, and the power outages time was the sum of two parts.


### 3.5.5. Simulation Process

(1)
Input the network structure; use the method in Section3.5.3, automatically traversing all the nonfaulty component’s feeder segments classification by computer; establish FMEA table.
(2)
Set initial analog clock value for 0; at this time all system components are in normal state. For each nonsource system component, generate a random number; use the state transition model to generate TTF for each component.
(3)
Find the component with the smallest TTF (TTFmin⁡), and advance the analog clock to TTFmin⁡. If some components are of the same TTF, randomly select one of them.
(4)
Generate a new random number and convert it into TTR for this component. At the same time, generate power isolation and recovery time (RT).
(5)
By using the FMEA table generated in step 1, check the feeder segments under this component’s outage. Fault segment load point outage time was TTR; upstream quarantine segment load point outage time was RT; record the time duration of power outages and the insufficiency of power supply.
(6)
For each upstream active island segment, sample DG’s state and BESS’s SOC atTTFmin⁡; calculate real-time output of DG and load during [TTFmin⁡, TTFmin⁡ + RT]; use the method in Section 3.5.4 to shed load; record the load point outage frequency and duration.
(7)
Similar to step 6, in [TTFmin⁡, TTFmin⁡ + RT], simulate the downstream active island; in [TTFmin⁡ + RT, TTFmin⁡ + TTR], simulate the downstream big island and record the load point outage frequency and duration.
(8)
For the faulty component in step 3, regenerate a random number and convert it intoTTFN for a new running time. At this moment, the component’s TTF is updated to TTFmin⁡ + TTR + TTFN.
(9)
If the analog clock is less than predetermined simulation time, return to step 3; otherwise, statistically count the outage frequency, duration, and insufficiency of each load point, and then calculate the reliability indices.


## 3.5.1. Sequential Sampling for Nonsource Component

All nonsource component was simulated by the two-state Markov model shown as in Figure1, the fault transfer rate was λ, and the recovery rate was μ, so the fault possibility was Pd=λ/(λ+μ), the regular state possibility was Pu=μ/(λ+μ), time between failures (TTF) and fault repair time (TTR) of components are exponentially distributed, and the Monte Carlo simulation method for state duration sampling was adopted for its sequential sampling [20].Figure 1
The two-state Markov model.


## 3.5.2. Sampling for DG

When not considering the limitation of upper system capacity, only when the nonsource components in system is in fault, the outage of DG and BESS can 6 exacerbate the power off risk of the user. Compared with normal state, the nonsource components fault duration is very short, and DG and BESS states can be considered unchanged during the fault, so their state can be sampled nonsequentially.Single wind turbine is still using two-state model of Figure1. For wind turbines set consisting of m wind turbines unit, m random number between [0,1] is first generated on behalf of its sample state, and Ri>Pd represents that the unit is working.PV array is consisting of a set of photovoltaic panels, to reflect a portion of PV array’s fault, on the basis of two-state model; the derating state is added, and its probability isPe. Sampling method was shown in Figure 2; first, the probability of all the states was placed cumulatively on number of axes between [0,1], and the random number R represents which state the PV array is in. For example, in Figure 3, the photovoltaic array is in derating state.Figure 2
Sampling process for the three-state element.Figure 3
A distribution system with WTG, PV, and BESS.The sample of the BESS state consists of two parts: one is sample of the running state, and the other is sample of the state of charge (SOC). Sampling methods of running state are the same with PV array, and derating state represents part of BESS’s failure. When sampling the SOC, SOC was converted from a random number in accordance with the SOC changing law under integration mode.


## 3.5.3. Failure Mode and Effect Analysis (FMEA) Process

FMEA analysis is the basis of reliability assessment. The DG’s integration will have an impact on the traditional FMEA process. In this paper, feeder segment concept was used to analyze the FMEA process of distribution system with DG. Feeder segment is a set of equipment in which all equipment has a common entrance [21] and all the load points in the same segment have the same failure consequences. For an example of the distribution system as Figure 3, the system can be divided into 12 feeder segments (S1~S12).Assuming the upper system’s capacity is adequate, when the system is in a normal state, all load points can be supplied by upper system and DG at the same time. When DG is in fault, the circuit breaker of DG can act immediately to isolate DG, and no load will suffer power loss. And only when the nonsource components fail, the load point will face the risk of power loss. At this point, the active islanding can narrow outage range and reduce outage time. However, the realization of the active islanding needs the adjusting of traditional protection. Therefore, the protection logic in paper [22] was used in this paper: for any faulty component of a feeder, only its nearest surrounding breaker/isolator switch will act for fault isolation.Still use Figure3 system as example, when the components in S5 (feeders, distribution transformers) was in fault, the acting circuit breakers should be B3, B4, B5, B6, B8, B9, B10, B12, and B13. Then, the isolation switches D2, D3, and D4 will act to isolate S5. According to the differences in each feeder segment’s power restoration result, the system can be divided into six kinds of segment as follows.Fault Segment ({S5}). Outage time for all load points in the region is the fault repair time.Unaffected Segment ({S1}). Load point in the region will not be affected.Upstream Quarantine Segment ({S2},{S3}). Located at upstream of fault segment, before fault isolation, load point in this region cannot restore power. After fault isolation, S3 as an example, isolating switch D1 and circuit breaker B5 will close, and both the DG and the system will supply load in S3. All the load point’s outage time in this region is fault isolation time.Upstream Active Islanding Segment ({S4}). When a fault occurs, due to the operation of the circuit breaker (B6), the region can immediately run in island mode, the island time is fault isolation time, the load point outage situation was decided by the island's power balance status, and when necessary, load shedding will be carried out, which would be described later. After fault isolation was over, the area can be reaccessed into systems.Downstream Quarantine Segment ({S6},{S7},{S9}). Before fault isolation, the region cannot run in island mode, and all load points suffer outage. After fault isolation, the region can be linked with other downstream segments to form a big island. For example, after fault isolation, feeder segments S7~S12 can form a big island, and island operation time was the difference between fault repair time and fault isolation time.Downstream Active Islanding Segment ({S8},{S10~S12}). Same as the upstream active islanding segment, the island operation time is divided into two parts: one is the seamless operation time before fault isolation; the other is the big island operation time after fault isolation.In addition, if the fault component is switch, we can simply merge feeder segments which connected by switch into one feeder segment and then reuse the FMEA analysis method above.


## 3.5.4. Load Shedding Strategy

Load shedding determination conditions are as follows: at any moment of battery discharge while islanding, the greatest power all the battery can release was no less than the net exchange power in island. Consider(1)∑i=1nPout,1i(t)ηd≥Pex(t),Pex(t)≥0,t∈[tst,tend].Pout,1i(t) is the maximum power which can be released by the battery pack i during time interval t and n is the number of battery packs; Pex(t)=L(t)+Sl(t)-DG(t), in which L(t) is the total load; DG(t) is the sum of all DG outputs; Sl(t) is the active power loss within island and, to simplify the flow calculation, Sl(t) can be set at about 5% of the current load [23].In order to improve simulation speed, this paper takes the following heuristic load shedding strategies. First, assuming that all load points are not shed, when the case does not meet formula (1), the smallest load was shed first; repeat this process until it satisfies formula (1). During load shedding process, the dynamic switching was not considered, so the shedding state was unchanged during islanding time. When formula (1) was satisfied, each battery power was released by the proportional ratio of Pout,1i(t) and power to be absorbed was by proportional ratio of max⁡(Pin,1i(t),-Pin,2i(t),-Pin,3i(t)). It should be noticed that, for the load points in downstream quarantine segment and downstream active islanding segment, if a load point was shed in both parts of the two island times, the frequency of power-off is 1, and the power outages time was the sum of two parts.


## 3.5.5. Simulation Process

(1)
Input the network structure; use the method in Section3.5.3, automatically traversing all the nonfaulty component’s feeder segments classification by computer; establish FMEA table.
(2)
Set initial analog clock value for 0; at this time all system components are in normal state. For each nonsource system component, generate a random number; use the state transition model to generate TTF for each component.
(3)
Find the component with the smallest TTF (TTFmin⁡), and advance the analog clock to TTFmin⁡. If some components are of the same TTF, randomly select one of them.
(4)
Generate a new random number and convert it into TTR for this component. At the same time, generate power isolation and recovery time (RT).
(5)
By using the FMEA table generated in step 1, check the feeder segments under this component’s outage. Fault segment load point outage time was TTR; upstream quarantine segment load point outage time was RT; record the time duration of power outages and the insufficiency of power supply.
(6)
For each upstream active island segment, sample DG’s state and BESS’s SOC atTTFmin⁡; calculate real-time output of DG and load during [TTFmin⁡, TTFmin⁡ + RT]; use the method in Section 3.5.4 to shed load; record the load point outage frequency and duration.
(7)
Similar to step 6, in [TTFmin⁡, TTFmin⁡ + RT], simulate the downstream active island; in [TTFmin⁡ + RT, TTFmin⁡ + TTR], simulate the downstream big island and record the load point outage frequency and duration.
(8)
For the faulty component in step 3, regenerate a random number and convert it intoTTFN for a new running time. At this moment, the component’s TTF is updated to TTFmin⁡ + TTR + TTFN.
(9)
If the analog clock is less than predetermined simulation time, return to step 3; otherwise, statistically count the outage frequency, duration, and insufficiency of each load point, and then calculate the reliability indices.


## 3.6. Simulation of High Penetration

### 3.6.1. Feature of High Penetration System

High penetration distribution network with DG will have both power generation and distribution properties, so-called active distribution network, and the key of its reliability assessment is how to properly combine power generation and distribution properties together. To solve this problem, this paper uses the sequential Monte Carlo sampling method based on the system state transition sampling to conduct the reliability assessment for distributed network in high DG penetration level, which can simultaneously calculate the reliability index of power generation and distribution.Compared to the commonly used state duration sampling method, system state transition sampling method is to sample the state of the entire system rather than of each specific component; there is no need to store the run-outage state cycling sequence of each component within the system, so it can greatly save data storage space and improve the simulation speed, which is more suitable for system with a large number of components, especially the systems including much power generation equipment.


### 3.6.2. System State Transition Sampling Method

The system state transition sampling method was focused on the state of the whole system rather than that of one component in the system. Assume that the system is composed ofm components and consists of n states. The state space of the system is G={S1,…,Sn}. Suppose that the present system state is Sk and the state transition rate of each component in the system is λi(i=1,…,m). Since the state duration of each component follows an exponential distribution, the duration of the system state Sk also follows the exponential distribution. The duration of the system state Sk has the following probability density function:
(2)f(t)=∑i=1mλiexp⁡⁡(-∑i=1mλit).It can be seen that the state transition rate to the next state from the system stateSk is the sum of each component’s state transition rate in the system state Sk.State transition of any component in the system can lead to a system state. Consequently, system hasm possible states starting from state Sk. The probability that the system reaches one of these possible states is given by the following equation:
(3)Pj=λj∑i=1mλi.λj is the state transition rate of the jth component. It is obvious that
(4)∑j=1mPj=1.Two uniformly distributed random numbersR1 and R2 between [0,1] are generated to sample the duration of the system state Sk and the system state which is to be reached. The duration of the system state Sk is given by the following equation:
(5)Tk=-ln⁡R1∑i=1mλi.The next system state can be determined by the sampling shown in Figure4. The probabilities of m possible future states are successively placed in the interval [0,1]. If R2 falls into the segment corresponding to Pj, the transition of the jth component leads system to the next system state.Figure 4
Explanation of system state transition sampling.


### 3.6.3. Analysis and Simulation of Generation State

According to different states of the components within the system, the state space of the distribution system with distributed generationG can be divided into two subsets G1 and G2. In the set G1, all nonsource components are in normal state, and in the set G2, a nonsource element is in fault state.In the traditional distribution system, when all nonsource components are in normal state, the system load can be supplied without considering the higher level power supply capacity. But the distribution system with DG is a multiple supply network, in which the load will be supplied by both the DG and the higher level power. Because DG can supply part of the load, the substation capacity will be reduced in order to save cost of grid construction. So it is necessary to consider the impact of the higher level power supply capacity. In this case, whether the total generation of the system including DG can supply all the load needs to be assessed. It is in fact a part of the generation system reliability evaluation.The reliability of generation system can be evaluated by simple simulation. The reliability indices are determined by superposing the chronological load curve on the generation capacity curve. Assume that the supply from higher level isPs(t) and the DG provides PDG(t) during the tth period. The sum of the real time load and the loss is PL(t) (assume that the loss is 5% of the load to simplify the power flow calculation [24]). It means that the system is suffering power shortage when Ps(t)+PDG(t)<PL(t) and the ENS can be calculated by ENS=PL(t)-Ps(t)-PDG(t). The indices of generation system such as LOLE and LOEE are calculated with sequential simulation during all periods. The method will be improved with the consideration of battery charging and discharging strategy. Two dispatch strategies of battery used in the paper are as follows [25].Load-Following Strategy. Under the load-following strategy, renewable power sources charge the BESS but the generators do not. The BESS discharge when DG could not provide enough supply to the load. The generators serve as the backup.Cycle Charge Strategy. When DG output was surplus, the BESS was charged by excess power; when DG output was insufficient, if the maximum output of DG and BESS still cannot meet all the loads’ need, the system will first supply the remaining load and then charge the BESS. In addition, to prevent the battery from being in a low SOC, once the battery starts charging, it needs to be charged to a specified SOCset.Figures5 and 6 show the simulation process in load-following strategy and in cycle charging strategy, respectively. Pdi(t) is the maximum output of the ith BESS during the tth period. Pci(t) is the maximum active power that the battery could absorb. Considering the difference in the capacity and SOC between batteries, the exchanging power of each battery is proportional to its Pdi(t) or Pci(t). For example, when the DG is charging all the batteries, the ith battery can get
(6)Pexi(t)=(PDG(t)-PL(t))×Pci(t)∑i=1mPci(t).Figure 5
Simulation process in load-following strategy.Figure 6
Simulation process in cycle charging strategy.However, the powerPini(t) that it really absorbs is min⁡{ηcPexi(t),Pci(t)}. ηc is charge efficiency.


### 3.6.4. Simulation Process

The time sequential Monte Carlo simulation using system state transition sampling evaluates the reliability of distribution system containing distributed wind energy, solar energy, and BESS. The step for simulation is 1 hour. The procedure used consists of the following steps.(1)
The simulation process starts from the normal system state in which all generating units and transmission components are in the upstate at timet0=0.(2)
Calculate the duration of the current system stateTk by (4) with the transition rate λi of each component. Get λi of a three-state component by adding its two state transition rates. If the present system state is a state in which at least one nonpower component is in fault state, set the transition rate of other components to zero.(3)
If all the nonpower components are in the normal state in the current state, call the generation system simulation process, and recordLLDk and ENSk during Tk. If there is a nonpower component in fault state, call the distribution system simulation process shown in Section 3.5 and record the indices λki, γki, and ENSki.(4)
Sett0=t0+Tk. If the time is shorter than the year limit, sample the next system state by the method in Figure 1 and go to step 2. If t0 is longer than N, stop simulating and calculate the following indices.Loss of load expectation (h/year) is(7)LOLE=1N∑k=1PLLDk.
Loss of expected energy (MWh/year) is(8)LOEE=1N∑k=1PENSk.
System average interruption frequency index (times/node* year) is
(9)SAIFI=(1/N)∑k=1Q∑i=1nλkiCi∑i=1nCi.
System average interruption duration index (h/node* year) is
(10)SAIDI=(1/N)∑k=1Q∑i=1nγkiCi∑i=1nCi.
Expected energy not supplied (MWh/year) is(11)EENS=1N∑k=1Q∑i=1nENSki.
P is the simulation times of the generation system state, Q is the simulation times of the distribution system state, Ci is the amount of users of each load point, and n is the amount of load points.


## 3.6.1. Feature of High Penetration System

High penetration distribution network with DG will have both power generation and distribution properties, so-called active distribution network, and the key of its reliability assessment is how to properly combine power generation and distribution properties together. To solve this problem, this paper uses the sequential Monte Carlo sampling method based on the system state transition sampling to conduct the reliability assessment for distributed network in high DG penetration level, which can simultaneously calculate the reliability index of power generation and distribution.Compared to the commonly used state duration sampling method, system state transition sampling method is to sample the state of the entire system rather than of each specific component; there is no need to store the run-outage state cycling sequence of each component within the system, so it can greatly save data storage space and improve the simulation speed, which is more suitable for system with a large number of components, especially the systems including much power generation equipment.


## 3.6.2. System State Transition Sampling Method

The system state transition sampling method was focused on the state of the whole system rather than that of one component in the system. Assume that the system is composed ofm components and consists of n states. The state space of the system is G={S1,…,Sn}. Suppose that the present system state is Sk and the state transition rate of each component in the system is λi(i=1,…,m). Since the state duration of each component follows an exponential distribution, the duration of the system state Sk also follows the exponential distribution. The duration of the system state Sk has the following probability density function:
(2)f(t)=∑i=1mλiexp⁡⁡(-∑i=1mλit).It can be seen that the state transition rate to the next state from the system stateSk is the sum of each component’s state transition rate in the system state Sk.State transition of any component in the system can lead to a system state. Consequently, system hasm possible states starting from state Sk. The probability that the system reaches one of these possible states is given by the following equation:
(3)Pj=λj∑i=1mλi.λj is the state transition rate of the jth component. It is obvious that
(4)∑j=1mPj=1.Two uniformly distributed random numbersR1 and R2 between [0,1] are generated to sample the duration of the system state Sk and the system state which is to be reached. The duration of the system state Sk is given by the following equation:
(5)Tk=-ln⁡R1∑i=1mλi.The next system state can be determined by the sampling shown in Figure4. The probabilities of m possible future states are successively placed in the interval [0,1]. If R2 falls into the segment corresponding to Pj, the transition of the jth component leads system to the next system state.Figure 4
Explanation of system state transition sampling.


## 3.6.3. Analysis and Simulation of Generation State

According to different states of the components within the system, the state space of the distribution system with distributed generationG can be divided into two subsets G1 and G2. In the set G1, all nonsource components are in normal state, and in the set G2, a nonsource element is in fault state.In the traditional distribution system, when all nonsource components are in normal state, the system load can be supplied without considering the higher level power supply capacity. But the distribution system with DG is a multiple supply network, in which the load will be supplied by both the DG and the higher level power. Because DG can supply part of the load, the substation capacity will be reduced in order to save cost of grid construction. So it is necessary to consider the impact of the higher level power supply capacity. In this case, whether the total generation of the system including DG can supply all the load needs to be assessed. It is in fact a part of the generation system reliability evaluation.The reliability of generation system can be evaluated by simple simulation. The reliability indices are determined by superposing the chronological load curve on the generation capacity curve. Assume that the supply from higher level isPs(t) and the DG provides PDG(t) during the tth period. The sum of the real time load and the loss is PL(t) (assume that the loss is 5% of the load to simplify the power flow calculation [24]). It means that the system is suffering power shortage when Ps(t)+PDG(t)<PL(t) and the ENS can be calculated by ENS=PL(t)-Ps(t)-PDG(t). The indices of generation system such as LOLE and LOEE are calculated with sequential simulation during all periods. The method will be improved with the consideration of battery charging and discharging strategy. Two dispatch strategies of battery used in the paper are as follows [25].Load-Following Strategy. Under the load-following strategy, renewable power sources charge the BESS but the generators do not. The BESS discharge when DG could not provide enough supply to the load. The generators serve as the backup.Cycle Charge Strategy. When DG output was surplus, the BESS was charged by excess power; when DG output was insufficient, if the maximum output of DG and BESS still cannot meet all the loads’ need, the system will first supply the remaining load and then charge the BESS. In addition, to prevent the battery from being in a low SOC, once the battery starts charging, it needs to be charged to a specified SOCset.Figures5 and 6 show the simulation process in load-following strategy and in cycle charging strategy, respectively. Pdi(t) is the maximum output of the ith BESS during the tth period. Pci(t) is the maximum active power that the battery could absorb. Considering the difference in the capacity and SOC between batteries, the exchanging power of each battery is proportional to its Pdi(t) or Pci(t). For example, when the DG is charging all the batteries, the ith battery can get
(6)Pexi(t)=(PDG(t)-PL(t))×Pci(t)∑i=1mPci(t).Figure 5
Simulation process in load-following strategy.Figure 6
Simulation process in cycle charging strategy.However, the powerPini(t) that it really absorbs is min⁡{ηcPexi(t),Pci(t)}. ηc is charge efficiency.


## 3.6.4. Simulation Process

The time sequential Monte Carlo simulation using system state transition sampling evaluates the reliability of distribution system containing distributed wind energy, solar energy, and BESS. The step for simulation is 1 hour. The procedure used consists of the following steps.(1)
The simulation process starts from the normal system state in which all generating units and transmission components are in the upstate at timet0=0.(2)
Calculate the duration of the current system stateTk by (4) with the transition rate λi of each component. Get λi of a three-state component by adding its two state transition rates. If the present system state is a state in which at least one nonpower component is in fault state, set the transition rate of other components to zero.(3)
If all the nonpower components are in the normal state in the current state, call the generation system simulation process, and recordLLDk and ENSk during Tk. If there is a nonpower component in fault state, call the distribution system simulation process shown in Section 3.5 and record the indices λki, γki, and ENSki.(4)
Sett0=t0+Tk. If the time is shorter than the year limit, sample the next system state by the method in Figure 1 and go to step 2. If t0 is longer than N, stop simulating and calculate the following indices.Loss of load expectation (h/year) is(7)LOLE=1N∑k=1PLLDk.
Loss of expected energy (MWh/year) is(8)LOEE=1N∑k=1PENSk.
System average interruption frequency index (times/node* year) is
(9)SAIFI=(1/N)∑k=1Q∑i=1nλkiCi∑i=1nCi.
System average interruption duration index (h/node* year) is
(10)SAIDI=(1/N)∑k=1Q∑i=1nγkiCi∑i=1nCi.
Expected energy not supplied (MWh/year) is(11)EENS=1N∑k=1Q∑i=1nENSki.
P is the simulation times of the generation system state, Q is the simulation times of the distribution system state, Ci is the amount of users of each load point, and n is the amount of load points.


## 4. Study Case

In this paper, our analysis is the use of the multibranch feeders in the transformed IEEE RBTS Bus 6 as examples, as shown in Figure7, including a segment bus (substation), 30 feeder segments, 26 nodes, 23 distribution transformers, 23 load points, 5 DGs, and several switches. Each DG includes several fans of the same type, a photovoltaic array, and a battery, fan type is Enercon E33, and the battery type is Hoppecke 24 OPzS 3000.Figure 7
Modified system under study.The length of each line segment and the peak of load point are the same as the literature [26], the total load peak is 11.653 MW, and the time-varying load models are the same as the literature [27]. The feeder failure rate is 0.065 times/year × km, the distribution transformer failure rate is 0.015 times/year, and the average repair times are 5 hours (repair rate is 0.2 times/hour), distributed exponentially. Single wind turbine failure rate is 4.6 times/year, and repair rate is 58.4 times/year. The transfer rate of PV array is the same as the battery, λ=1.4, μ=58.4, λd=3.2, μd=46.72, λf=11.68, μf=0; the unit is all times/year; the derating operation capacity is half of the full state. Bus and switches are completely reliable, fault isolation and recovery time takes constantly 1 h, and the number of users for each load point is one.


### 4.1. Case of Low Penetration

#### 4.1.1. Impact of DG Capacity

System Average Interruption Time (SAIDI) is the most representative reliability index. In the case of no BESS, the system SAIDI’s changing with capacity change of WTG and PV is shown in Figure8. Abscissa represents the WTG and PV’s rated capacity. As can be seen, DG can improve the reliability level of the system, especially in its initial increasing stage. When DG capacity reaches to a certain extent, the effect of reliability improvement is close to saturation. In this paper’s wind speed and light intensity data examples, the WTG’s enhancement of reliability of the system is better than that of the PV.Figure 8
SAIDI variation as function of WTG/PV capacity.


#### 4.1.2. Impact of BESS Capacity

Take a few different sets of “WTG + PV” capacity combination; SAIDI index with the changes of BESS capacity was shown in Figure9. Abscissa represents the BESS’s rated capacity. As can be seen, in the initial stage of BESS capacity’s increase, the SAIDI is lower than the minimum value in Figure 6, but, with increasing capacity of DG, the SAIDI’s decline rate in the initial stage is increasing, which reflects the ability of BESS’s smoothing renewable DG output and its potential to enhance the reliability. In addition, with further increases of DG capacity, the initial decline speed tends to be the same, and, with increasing of BESS capacity, SAIDI’s decline rate is getting slower and slower. We can see that DG and BESS’s capacity and ratio should be allocated rationally based on the actual needs.Figure 9
SAIDI variation as function of battery capacity.


### 4.2. Case of High Penetration

#### 4.2.1. Impact of Substation Capacity

We still take WTG, PV, and BESS in each DG; the capacity of WTG is 1005 kw, the capacity of the PV is 1005 kw, and the capacity of the battery is 13200 kwh, when calculated under different supply substation capacity constraints using this method; the system reliability indices are shown as in Table1. We can see that when the capacity of the supply substation is less than the total system load peak (11.653 MW), the generating indicators LOLE and LOEE cease to be 0, indicating that the system will face a shortage risk. With the decrease of the supply substation capacity, power generation and distribution reliability indices are increasing, and especially the growth of power generation indices is more sensitive.Table 1
Impacts of substation capacity.
Charge and discharge strategy
Substation capacity (MW)
LOLE
LOEE
SAIFI
SAIDI
EENS
Load-following
13
0
0
1.3815
6.6364
32.671
11
9.3812
3.4664
1.3965
6.7033
32.968
9
212.24
156.94
1.3971
6.7549
33.223
7
849.61
1174.6
1.3999
6.7548
33.259
5
2128
3932.2
1.5
7.9927
38.353
3
5010.5
10653
1.7695
10.0764
47.568
1
8069
24362
2.1862
12.9036
60.861
Cycle charge
13
0
0
1.1001
5.0396
24.621
11
0.7126
0.2931
1.1125
5.0485
24.804
9
20.43
14.24
1.1138
5.0589
24.873
7
163.79
260.84
1.1347
5.2858
25.751
5
989.89
2268.3
1.2073
5.7744
28.048
3
3913.2
9145
1.4518
7.7969
36.992
1
7910.8
24161
2.0388
12.1559
57.085From Table1 we can see that if we appropriately reduce the capacity of the supply substation, the system can still maintain a certain level of reliability. Take cycle charging strategy, for example; when supply substation capacity reduces to 9 MW, the system's total lack of power supply LOEE + EENS will be 39.113 MWh/year, which is still lower than 49.290 MWh/year which is the number before the DG was accessed. This shows that after DGs access the distribution system, when there is no increase in the system load, the capacity of the supply substations may appropriately reduce. In other words, in the case of load growth, it can slow down the expansion of supply substation capacity. However, if the supply substation capacity reduces further or only relies on the renewable DG to supply the load growth, the level of system reliability will reduce. For example, when the supply substation capacity is lowered to 7 WM, EENS in cycle charging strategy will surge to 286.591 MWh/year.


#### 4.2.2. Impact of Substation and DG Capacity

Table2 shows the LOEE + EENS calculated with changing supply substation capacity and DG capacity in cycle charging strategy. The DGs access different locations. Column in Table 2 stands for the supply distribution capacity while row stands for total DG capacity. Each “×1” means that DG consists of WTG for 1005 kW, PV for 1005 kW, and BESS for 13200 kWh. The capacity of DG in “×2” is twice that in “×1.” The rest may be deduced by analogy.Table 2
Impact of substation and DG capacity.
Capacity
×1
×2
×3
×4
×5
×6
11
25.097
18.295
17.402
17.03
16.911
16.701
9
39.113
23.467
19.442
18.114
17.371
16.876
7
286.591
96.966
66.349
42.081
23.958
17.283
5
2296.35
852.53
406.22
170.22
58.804
21.871From Table2 we can see that in order to achieve a certain level of system reliability the capacity of supply substation and DG can be comprehensively selected. For example, when the requirement of EENS is about 23 MWh/year, we can assume that the capacity of the supply substation is 11 MW and the DG capacity is “×1”; we can also assume that the supply substation capacity is 5 MW and DG capacity is “×6.” However, due to the volatility of the output of renewable DG, in order to maintain the reliability level of the system, more DGs are needed to compensate for the reduced superior substation capacity. Compared with the former, DGs in the latter became the main power supply in the system; meanwhile the supply substation capacity can decrease by 6 MW, but the total capacity of the wind turbines and photovoltaic in the system will increase approximately by 25 MW, while battery capacity will increase approximately by 330 MWh. Therefore, when DG is the main power supply of the system, in order to ensure a certain degree of reliability, the total rated capacity of DG should be much higher than the system total load peak, and the cost is higher. According to the actual situation, we should reasonably optimize the capacity of supply substation and DG in order to meet the requirements of reliability and economy at the same time.


## 4.1. Case of Low Penetration

### 4.1.1. Impact of DG Capacity

System Average Interruption Time (SAIDI) is the most representative reliability index. In the case of no BESS, the system SAIDI’s changing with capacity change of WTG and PV is shown in Figure8. Abscissa represents the WTG and PV’s rated capacity. As can be seen, DG can improve the reliability level of the system, especially in its initial increasing stage. When DG capacity reaches to a certain extent, the effect of reliability improvement is close to saturation. In this paper’s wind speed and light intensity data examples, the WTG’s enhancement of reliability of the system is better than that of the PV.Figure 8
SAIDI variation as function of WTG/PV capacity.


### 4.1.2. Impact of BESS Capacity

Take a few different sets of “WTG + PV” capacity combination; SAIDI index with the changes of BESS capacity was shown in Figure9. Abscissa represents the BESS’s rated capacity. As can be seen, in the initial stage of BESS capacity’s increase, the SAIDI is lower than the minimum value in Figure 6, but, with increasing capacity of DG, the SAIDI’s decline rate in the initial stage is increasing, which reflects the ability of BESS’s smoothing renewable DG output and its potential to enhance the reliability. In addition, with further increases of DG capacity, the initial decline speed tends to be the same, and, with increasing of BESS capacity, SAIDI’s decline rate is getting slower and slower. We can see that DG and BESS’s capacity and ratio should be allocated rationally based on the actual needs.Figure 9
SAIDI variation as function of battery capacity.


## 4.1.1. Impact of DG Capacity

System Average Interruption Time (SAIDI) is the most representative reliability index. In the case of no BESS, the system SAIDI’s changing with capacity change of WTG and PV is shown in Figure8. Abscissa represents the WTG and PV’s rated capacity. As can be seen, DG can improve the reliability level of the system, especially in its initial increasing stage. When DG capacity reaches to a certain extent, the effect of reliability improvement is close to saturation. In this paper’s wind speed and light intensity data examples, the WTG’s enhancement of reliability of the system is better than that of the PV.Figure 8
SAIDI variation as function of WTG/PV capacity.


## 4.1.2. Impact of BESS Capacity

Take a few different sets of “WTG + PV” capacity combination; SAIDI index with the changes of BESS capacity was shown in Figure9. Abscissa represents the BESS’s rated capacity. As can be seen, in the initial stage of BESS capacity’s increase, the SAIDI is lower than the minimum value in Figure 6, but, with increasing capacity of DG, the SAIDI’s decline rate in the initial stage is increasing, which reflects the ability of BESS’s smoothing renewable DG output and its potential to enhance the reliability. In addition, with further increases of DG capacity, the initial decline speed tends to be the same, and, with increasing of BESS capacity, SAIDI’s decline rate is getting slower and slower. We can see that DG and BESS’s capacity and ratio should be allocated rationally based on the actual needs.Figure 9
SAIDI variation as function of battery capacity.


## 4.2. Case of High Penetration

### 4.2.1. Impact of Substation Capacity

We still take WTG, PV, and BESS in each DG; the capacity of WTG is 1005 kw, the capacity of the PV is 1005 kw, and the capacity of the battery is 13200 kwh, when calculated under different supply substation capacity constraints using this method; the system reliability indices are shown as in Table1. We can see that when the capacity of the supply substation is less than the total system load peak (11.653 MW), the generating indicators LOLE and LOEE cease to be 0, indicating that the system will face a shortage risk. With the decrease of the supply substation capacity, power generation and distribution reliability indices are increasing, and especially the growth of power generation indices is more sensitive.Table 1
Impacts of substation capacity.
Charge and discharge strategy
Substation capacity (MW)
LOLE
LOEE
SAIFI
SAIDI
EENS
Load-following
13
0
0
1.3815
6.6364
32.671
11
9.3812
3.4664
1.3965
6.7033
32.968
9
212.24
156.94
1.3971
6.7549
33.223
7
849.61
1174.6
1.3999
6.7548
33.259
5
2128
3932.2
1.5
7.9927
38.353
3
5010.5
10653
1.7695
10.0764
47.568
1
8069
24362
2.1862
12.9036
60.861
Cycle charge
13
0
0
1.1001
5.0396
24.621
11
0.7126
0.2931
1.1125
5.0485
24.804
9
20.43
14.24
1.1138
5.0589
24.873
7
163.79
260.84
1.1347
5.2858
25.751
5
989.89
2268.3
1.2073
5.7744
28.048
3
3913.2
9145
1.4518
7.7969
36.992
1
7910.8
24161
2.0388
12.1559
57.085From Table1 we can see that if we appropriately reduce the capacity of the supply substation, the system can still maintain a certain level of reliability. Take cycle charging strategy, for example; when supply substation capacity reduces to 9 MW, the system's total lack of power supply LOEE + EENS will be 39.113 MWh/year, which is still lower than 49.290 MWh/year which is the number before the DG was accessed. This shows that after DGs access the distribution system, when there is no increase in the system load, the capacity of the supply substations may appropriately reduce. In other words, in the case of load growth, it can slow down the expansion of supply substation capacity. However, if the supply substation capacity reduces further or only relies on the renewable DG to supply the load growth, the level of system reliability will reduce. For example, when the supply substation capacity is lowered to 7 WM, EENS in cycle charging strategy will surge to 286.591 MWh/year.


### 4.2.2. Impact of Substation and DG Capacity

Table2 shows the LOEE + EENS calculated with changing supply substation capacity and DG capacity in cycle charging strategy. The DGs access different locations. Column in Table 2 stands for the supply distribution capacity while row stands for total DG capacity. Each “×1” means that DG consists of WTG for 1005 kW, PV for 1005 kW, and BESS for 13200 kWh. The capacity of DG in “×2” is twice that in “×1.” The rest may be deduced by analogy.Table 2
Impact of substation and DG capacity.
Capacity
×1
×2
×3
×4
×5
×6
11
25.097
18.295
17.402
17.03
16.911
16.701
9
39.113
23.467
19.442
18.114
17.371
16.876
7
286.591
96.966
66.349
42.081
23.958
17.283
5
2296.35
852.53
406.22
170.22
58.804
21.871From Table2 we can see that in order to achieve a certain level of system reliability the capacity of supply substation and DG can be comprehensively selected. For example, when the requirement of EENS is about 23 MWh/year, we can assume that the capacity of the supply substation is 11 MW and the DG capacity is “×1”; we can also assume that the supply substation capacity is 5 MW and DG capacity is “×6.” However, due to the volatility of the output of renewable DG, in order to maintain the reliability level of the system, more DGs are needed to compensate for the reduced superior substation capacity. Compared with the former, DGs in the latter became the main power supply in the system; meanwhile the supply substation capacity can decrease by 6 MW, but the total capacity of the wind turbines and photovoltaic in the system will increase approximately by 25 MW, while battery capacity will increase approximately by 330 MWh. Therefore, when DG is the main power supply of the system, in order to ensure a certain degree of reliability, the total rated capacity of DG should be much higher than the system total load peak, and the cost is higher. According to the actual situation, we should reasonably optimize the capacity of supply substation and DG in order to meet the requirements of reliability and economy at the same time.


## 4.2.1. Impact of Substation Capacity

We still take WTG, PV, and BESS in each DG; the capacity of WTG is 1005 kw, the capacity of the PV is 1005 kw, and the capacity of the battery is 13200 kwh, when calculated under different supply substation capacity constraints using this method; the system reliability indices are shown as in Table1. We can see that when the capacity of the supply substation is less than the total system load peak (11.653 MW), the generating indicators LOLE and LOEE cease to be 0, indicating that the system will face a shortage risk. With the decrease of the supply substation capacity, power generation and distribution reliability indices are increasing, and especially the growth of power generation indices is more sensitive.Table 1
Impacts of substation capacity.
Charge and discharge strategy
Substation capacity (MW)
LOLE
LOEE
SAIFI
SAIDI
EENS
Load-following
13
0
0
1.3815
6.6364
32.671
11
9.3812
3.4664
1.3965
6.7033
32.968
9
212.24
156.94
1.3971
6.7549
33.223
7
849.61
1174.6
1.3999
6.7548
33.259
5
2128
3932.2
1.5
7.9927
38.353
3
5010.5
10653
1.7695
10.0764
47.568
1
8069
24362
2.1862
12.9036
60.861
Cycle charge
13
0
0
1.1001
5.0396
24.621
11
0.7126
0.2931
1.1125
5.0485
24.804
9
20.43
14.24
1.1138
5.0589
24.873
7
163.79
260.84
1.1347
5.2858
25.751
5
989.89
2268.3
1.2073
5.7744
28.048
3
3913.2
9145
1.4518
7.7969
36.992
1
7910.8
24161
2.0388
12.1559
57.085From Table1 we can see that if we appropriately reduce the capacity of the supply substation, the system can still maintain a certain level of reliability. Take cycle charging strategy, for example; when supply substation capacity reduces to 9 MW, the system's total lack of power supply LOEE + EENS will be 39.113 MWh/year, which is still lower than 49.290 MWh/year which is the number before the DG was accessed. This shows that after DGs access the distribution system, when there is no increase in the system load, the capacity of the supply substations may appropriately reduce. In other words, in the case of load growth, it can slow down the expansion of supply substation capacity. However, if the supply substation capacity reduces further or only relies on the renewable DG to supply the load growth, the level of system reliability will reduce. For example, when the supply substation capacity is lowered to 7 WM, EENS in cycle charging strategy will surge to 286.591 MWh/year.


## 4.2.2. Impact of Substation and DG Capacity

Table2 shows the LOEE + EENS calculated with changing supply substation capacity and DG capacity in cycle charging strategy. The DGs access different locations. Column in Table 2 stands for the supply distribution capacity while row stands for total DG capacity. Each “×1” means that DG consists of WTG for 1005 kW, PV for 1005 kW, and BESS for 13200 kWh. The capacity of DG in “×2” is twice that in “×1.” The rest may be deduced by analogy.Table 2
Impact of substation and DG capacity.
Capacity
×1
×2
×3
×4
×5
×6
11
25.097
18.295
17.402
17.03
16.911
16.701
9
39.113
23.467
19.442
18.114
17.371
16.876
7
286.591
96.966
66.349
42.081
23.958
17.283
5
2296.35
852.53
406.22
170.22
58.804
21.871From Table2 we can see that in order to achieve a certain level of system reliability the capacity of supply substation and DG can be comprehensively selected. For example, when the requirement of EENS is about 23 MWh/year, we can assume that the capacity of the supply substation is 11 MW and the DG capacity is “×1”; we can also assume that the supply substation capacity is 5 MW and DG capacity is “×6.” However, due to the volatility of the output of renewable DG, in order to maintain the reliability level of the system, more DGs are needed to compensate for the reduced superior substation capacity. Compared with the former, DGs in the latter became the main power supply in the system; meanwhile the supply substation capacity can decrease by 6 MW, but the total capacity of the wind turbines and photovoltaic in the system will increase approximately by 25 MW, while battery capacity will increase approximately by 330 MWh. Therefore, when DG is the main power supply of the system, in order to ensure a certain degree of reliability, the total rated capacity of DG should be much higher than the system total load peak, and the cost is higher. According to the actual situation, we should reasonably optimize the capacity of supply substation and DG in order to meet the requirements of reliability and economy at the same time.


## 5. Conclusion

In this paper we have treated the reliability assessment problem of low and high DG penetration level of active distribution system using Monte Carlo simulation method. The problem is formulated as a two-case program, the program of low penetration simulation and the program of high penetration simulation. The load shedding strategy, and the simulation process was introduced in detail during each FMEA process.Results indicate that the integration of DG can improve the reliability of the system if the system was operated actively.In the case of low DG penetration, DG’s integration improves the reliability of system. In the case of high DG penetration, for economic considerations, the supply substation capacity of distribution networks should be appropriately reduced, resulting in the fact that the supply substation capacity cannot be assumed unlimitedly in the reliability assessment. In this case, the distribution network with DG appeared to be different from the traditional distribution network of some new features, which have the property of both distribution system and generation system. This method can be applied to the planning and coordination of substation, distribution network, and DG and also to the optimization design of the microgrid.



---

*Source: 421347-2014-05-25.xml*