The model at the planning investment level is to solve the optimal DRGs allocation scheme under the disaster attack scenario. The objective function consider both the investment and construction cost. The objective function is as follows: C i n v = ∑ n ∈ Ω N c p v x n p v + c w t x n w t (2) where c p v and c w t are the cost coefficients of PV units and WTs, respectively; x n p v and x n w t are the decision variables for optimal allocation at node n ; Ω N the set of nodes in the system.

The constraints at the planning investment level are mainly the allocation budget constraints of DRGs, including the cardinality budget and the monetary budget. The mathematical form of the constraints is as follows: ∑ n ∈ Ω N x n p v ≤ G p v , ∑ n ∈ Ω N x n w t ≤ G w t ∑ n ∈ Ω N c p v x n p v + c w t x n w t ≤ C b u d g e t (3) where G p v and G w t denote the available number of PV units and WTs respectively; C b u d g e t is the budget cost.

An ideal DRGs allocation scheme can not only reduce the output of carbon-emitting units and improve the low-carbon performance of the system operation, but also effectively prevent the risk of insufficient power supply and increased carbon emissions caused by disaster scenarios. Under a given allocation scheme and determined scenario conditions, the model at the simulated operation level minimizes the load loss in the disaster scenarios by means of NR, while also considering the economic efficiency and carbon emission mitigation. Therefore, the optimization objective includes the operation cost C E C O o p and the degree of load loss C L S o p under the multi-scenario disasters. In addition, the optimal allocation scheme should be able to mitigate the additional carbon emissions caused by disasters, so the optimization objective also includes the carbon emission cost C C O 2 o p . The objective function is as follows; C o p = C L S o p + C E C O o p + C C O 2 o p (4) C L S o p = ∑ s ∈ Φ φ s ∑ t ∈ T ∑ l ∈ Γ V o L L l L l , t , s p r e − L l , t , s (5) C E C O o p = ∑ s ∈ Φ φ s C s G + C s R E (6) C s G = ∑ t ∈ T ∑ g ∈ N G a g P g , t , s + b g (7) C s R E = ∑ t ∈ T ∑ m ∈ N R E λ R E P m , t , s (8) C C O 2 o p = ∑ s ∈ Φ φ s ⋅ p C O 2 ∑ t ∈ T ∑ g ∈ N G δ g P g , t , s + ∑ l ∈ Γ δ l s u b L l , t , s s u b (9)

where Φ denotes the set of disaster scenarios; φ s is the probability of occurrence for disaster scenario s ; N G , N R E and Γ are the set of carbon emission units, DRGs and loads respectively; C s G and C s R E denote the output cost of carbon emission units and the operation and maintenance cost of DRGs in the scenario s respectively; a g and b g are the cost coefficients; V o L L l is the value of lost load (VoLL) l ; L l , t , s p r e is the expected load value of load l during period t in the disaster scenario s ; L l , t , s is the actual load value of load l during period t in the disaster scenario s ; P m , t , s is the output of DRGs; P g , t , s is the output power of carbon emission units; λ R E is the operation and maintenance cost coefficient of DRGs; p C O 2 is the penalty cost per unit mass of CO2; δ g is the emission coefficient of carbon emission units; L l , t , s s u b is the output of power substitution on the consumer side; δ l s u b is the emission coefficient of energy substitution.

The constraints of the above model are as follows:

In order to evaluate the impact of the response scheme on the carbon emissions of the disaster-stricken distribution system, this paper first defines the total carbon emission (TCE) as an evaluation index to evaluate the carbon emission of the system under a given scheme. As previously analyzed, carbon emissions from rigid loads with production plan postponement need to be considered in post-disaster carbon accounting, so the TCE calculation formula is as follows: T C E = D ∑ s ∈ Φ φ s ∑ t ∈ T ∑ g ∈ N G δ g P g , t , s + ∑ l ∈ Γ δ l s u b L l , t , s s u b + ∑ l ∈ Γ δ g L l , t , s t r a n s (22) where D is the frequency of disasters in a year; L l , t , s t r a n s denotes the loss of rigid load with production plan postponement. This paper considers that the transfer load is powered by the substation after disaster restoration, so its carbon emission coefficient is δ g .

On this basis, in order to comprehensively consider the impact of disaster response schemes on power supply reliability and carbon emissions, a load loss emission ratio (LER) is defined to characterize the effect of disaster response measures. The smaller the value is, the smaller the load loss value corresponding to the unit carbon emission under the scheme is. That is, more loads are restored at the cost of smaller carbon emissions. The calculation formula of LER is as follows: L E R x = ∑ s ∈ Φ ∑ t ∈ T ∑ l ∈ Γ D L l , t , s p r e − L l , t , s T C E x (23) where T C E x is the total carbon emissions under response scheme x . It is worth noting that the value of LER is not the smaller the better. It is also necessary to comprehensively evaluate the scheme according to the actual load loss value and the total carbon emission.

The occurrence and development of natural disasters have certain regional and directional characteristics, and as a form of attack, its impact on the distribution system has certain uncertainty. Therefore, when analyzing the impact of natural disaster damage on the distribution system, it is necessary to fully consider the spatial and temporal characteristics of natural disaster development, clarify the disaster path in the distribution system, determine the scope of disaster impact on the distribution system, and evaluate the severity of disaster impact on the components in the distribution system.

Due to the strong correlation between the occurrence of natural disasters and geographical location, this paper combines the topology with the geographical information system (GIS), divides the system into square grids based on geographical areas (9), and draws a disaster frequency heat map as shown in Figure 2 based on the historical data of disasters. The red area in Figure 2 represents the disaster-prone area, and the deeper the red color, the higher the frequency of disasters, while the blue area represents the area where rare disasters occur. It can be seen from Figure 2 that the distribution system in the red area is more vulnerable to natural disasters and damage, resulting in fault.

According to the disaster frequency heat map, the average frequency of a certain type of disaster in each square grid can be calculated, so as to obtain the high frequency occurrence square grid of such disaster in the grid area. The occurrence and development of disasters usually have certain dynamic spatial and temporal distribution characteristics, so it is necessary to consider the impact of disaster development path. This paper assumes that the disaster moves geographically along the development path and affects the distribution system along the way, causing the power line on the moving path to fail. Figure 3 shows the path of disaster development in the distribution system. The disaster moves along different paths from the occurrence point, and the lines in the square grid on the path will be affected by the disaster. Because the square grids on the disaster path are in different geographical areas, the degree of disaster impact is not the same. This paper assumes that in the square grid near the disaster point, the power grid is more seriously affected, while the grid far away from the disaster point is relatively less affected by the disaster (18). According to the disaster path and the disaster intensity in the square grid, the fault probability curve can be used to calculate the fault probability of the line in the grid. The fault probability curve can usually be expressed as (Panteli et al., 2017): P s = 0 , I < I c r i t i c a l P s I , I c r i t i c a l ≤ I < I c o l l a p s e 1 , I ≥ I c o l l a p s e (24) where P s denotes the line fault probability in the grid related to the disaster intensity in scenario s ; I denotes the intensity of disaster; I c r i t i c a l is the critical disaster intensity that causes the line to fail; I c o l l a p s e is the intensity of the disaster that makes the line certain fault. Figure 4A shows the line fault probability curve under the variety of disaster intensity.

Because it is difficult to accurately measure and calculate the disaster intensity, in order to simplify the model and facilitate the solution, this paper uses the distance between the square grid on the disaster path and the disaster occurrence point to replace the disaster intensity as the input of the fault probability curve, and sets the line fault probability in each square grid to be the same. Therefore, the fault probability curve is discretized, and the discretized curve is shown in Figure 4B. In the figure, the disaster distance of the abscissa is represented by the square grid order along the disaster path.

The distribution system itself has a certain degree of resistance to natural disasters, which determines that the damage of disasters to the power system is not without an upper limit. Therefore, this should be taken into account when generating line fault constraints in disaster scenarios. The disaster scenario is modelled by using the line fault set considering the budget of disaster attack. For a certain type of disaster scenario, the following constraints should be met: ∑ i , j ∈ Ω R u i j s ≤ K b u d g e t (25) where u i j s indicates whether line i - j is affected by disasters and has faults. Ω R denotes the set of lines on the disaster path; K b u d g e t is the budget of disaster attack, the maximum number of fault lines under the disaster. Ignoring the sequence of the impact of disasters on the line, the model can be regarded as a typical N-K fault problem. The constraint Eq. 25 can be rewritten as follows: ∑ i , j ∈ Ω B 1 − u i j s ≥ N b r c h − K b u d g e t (26) where N b r c h is the total number of the branches in the system, containing the branches where the reconfiguration switches are located.

After considering the impact of disasters and NR, the final state of the line in the distribution system can be expressed as: z i j s = z i j 0 + U s w s ⋅ 1 − u i j s (27) where z i j 0 is the initial state of the line in the distribution system; U s w s is the state variable of the reconfiguration switch; It should be noted that when the line controlled by the reconfiguration switch is on the path affected by the disaster, if the disaster causes damage to the line, the line cannot operate normally even if the reconfiguration switch is closed.

For the DAD model proposed in this paper, in order to ensure the robustness of the allocation results, the most severe impact of natural disasters as attackers on the distribution system should be considered. Based on the N-K criterion, a method for selecting the worst-case fault scenario considering the line fault probability is proposed. This method evaluates the severity of the impact of a line fault on the distribution system by solving the risk of load loss (Nikkhah et al., 2018). The worst-case fault scenario selection process is shown in Figure 5.

The calculation formula of the risk of load loss (RoLL) mentioned in the above process is as follow: R o L L l = C L S o p ⋅ P l s (28) where R o L L l denotes the load loss risk value of line l ; C L S o p denotes the load loss after line l fault; P l s denotes the fault probability of line l under disaster scenario s .

After obtaining disaster information and related data, using the aforementioned disaster scenario modeling and worst-case fault scenario selection methods, it is possible to determine the specific attack scenario α in which the disaster acts as an attacker in the distribution system, corresponding to the second stage of the DAD model proposed in Section 3. Therefore, the min-max-min problem in this paper can be transformed into a mixed integer programming (MIP) optimization problem with both inner and outer layers in the form of min, so that it can be combined into a single-layer optimization model, which can be written as follow: min x , y s , z s   c T x + ∑ s ∈ Φ φ s d T y s s . t .   G x ≥ g K s y s ≥ h s F s z s ≥ m s Q s z s = r s I s y s + J s z s ≥ l s U s y s + T x y s = v s (29) where he specific expression of the optimization decision variables are: x = x n p v , x n w t T y s = L l , t , s , P g , t , s , P m , t , s , P b , t , s , P s w , t , s , θ i , t , s , F i j s , W i s T z s = U s w s , r i s T (30)

The first row of constraints in Eq. 29 corresponds to Eq. 3; the second row corresponds to the last row of Eqs 12, 13, 15, 17, 20; the third row corresponds to Eqs 19, 21; the fourth row corresponds to Eqs 18, 27; the fifth row corresponds to Eqs 11, 14, and the first five rows of Eq. 20; and the sixth row corresponds to Eq. 10.

It can be seen that there is a case where integer decision variables are multiplied by continuous decision variables in constraint Eq. 10, and it is difficult to solve this nonlinear constraint. Therefore, the following linear equivalent transformation is conducted for this constraint:

First, define the auxiliary variable w s , let: w s = A m P m , t , s (31)

Then, the following auxiliary constraints are introduced: w s ≤ P m , t , s w s ≥ P m , t , s − E − A m P m , t , s max 0 ≤ w s ≤ A m P m , t , s max (32)

Finally, the constraint Eq. 10 becomes: A g P g , t , s + w s = A l L l , t , s + A b , s P b , t , s + A s w P s w , t , s (33)

As can be seen from Eq. 29, the model is a hierarchical multi-scenario optimization problem. The investment level decision variable a is an ex-ante decision variable that is independent of the scenario, while the operation level decision variables b and c are related to the scenario and need to be calculated based on the specific disaster scenario. For the multi-scenario optimization problem, the progressive hedging algorithm (Rockafellar and Wets, 1991) is employed to solve it. PHA is a decomposition algorithm for multi-scenario optimization problems. Its main idea is to use orthogonal projection and augmented Lagrange multiplier method to decompose the original problem into sub-problems in multiple scenarios for iterative solution. The advantage of PHA is that it exhibits global convergence when dealing with convex optimization problems. For the optimization problem shown in Eq. 29, the PHA iteration steps are as follows:

1) Initialization. Set k = 0 , ω s k = 0 .

If γ k < ε , the iteration is terminated and then take x * = x ¯ k + 0.5 I E ; otherwise return to step 4.

In the above iteration process, k denotes the number of iterations. ω s k denotes the auxiliary multiplier of the k -th iteration under the scenario s . x s k denotes the investment level solution of the sub-problem solved by the k -th iteration under scenario s . x ¯ k is the weighted average of the solutions of the investment level under all scenario s ∈ Φ after the k -th iteration. ρ is a penalty parameter, and the value can refer to the reference (Watson and Woodruff, 2011). x * denotes the optimal solution of investment level. I E is a column vector whose elements are all 1. ⋅ denotes the rounding operation for matrix elements.

The modified IEEE 33-node distribution system is used for case study, as shown in Figure 6. In the system, node 16 and node 31 are connected with a micro turbine (MT) respectively. In addition, there are power substitutions at nodes 14, 22, and 24, and rigid load with load transfer at nodes 10, 17, and 27. The disaster scenario data and disaster geographic information are derived from the actual statistical data of a disaster-prone year in a city in China. Using the method proposed in Section 4.1 to process disaster data can obtain typical disaster scenario probabilities, disaster paths, and line fault probabilities on disaster paths, as shown in Table 1. Due to the obvious seasonal characteristics of the occurrence of disasters, the k-means clustering method is used to process the historical data of wind power, photovoltaic and load in the disaster-prone seasons. This paper mainly focuses on the summer with frequent floods and the spring with frequent wildfires, and obtains the typical intraday curves of wind power, photovoltaic and load under disaster scenarios, as shown in Figure 7. The VoLL of each node of the test system is shown in Table 2. The allocation budget of DRGs includes 2 photovoltaic units and 1 WT with a capacity of 400 kW, and the power generation costs are 0.35¥/(kW·h) and 0.28¥/(kW·h) respectively. The cost of purchasing power from the upstream power grid is 0.78¥/(kW·h), and other parameters are shown in Table 3. In this paper, the duration of a single disaster is considered to be 24 h, so the simulation time scale is set to 24 h.

In order to obtain the most severe line faults combination in the disaster scenario, traverse each line on the disaster path, set the line disconnected and calculate the RoLL value, and the results shown in Table 4 can be obtained. Set the disaster attack budget value K b u d g e t , and based on the results in Table 4, select the same number of lines with large RoLL values in the scenario as the worst-case fault scenario for the disaster scenario.

In order to verify the effectiveness of the model and disaster response measures proposed in this paper, five schemes are set up for comparison. The scheme settings are shown in Table 5. Set the disaster attack budget K b u d g e t to 2, and according to Table 4, the worst-case fault scenarios are: lines 5–6 and 7–8 (in scenario 1); lines 31–32 and 23–24 (in scenario 2); lines 6–7 and 26–27 (in scenario 3); lines 3–4 and 19–20 (in scenario 4). The results solved under the five different disaster response schemes are shown in Table 6.

By comparing the results of the scheme 1 and scheme 2, it can be seen that by adopting the NR, the load loss caused by the distribution system failure under the disaster is reduced by 81.94%, and the total operation cost is reduced by 5.34 × 10⁴¥. This indicates that NR can effectively restore most of the power supply after a distribution system failure, ensuring power supply reliability. Due to the NR restoring a large amount of load, resulting in an increase in the output of carbon emission units, the carbon emissions in scheme 2 have significantly increased compared to scheme 1. From the LER value, it can be seen that scheme 2 is superior to scheme 1.

From the results of scheme 3 and scheme 4, it can be concluded that compared with scheme 1, after the allocation of PV unit and WT, the load loss under disasters is reduced by 5,171.12 and 6,978.41 kW·h respectively, and the total carbon emission is reduced by 5.28% and 6.68% respectively. This is due to the fact that the output of DRGs has restored the lost load in some load islands, and the total carbon emission has decreased due to the replacement of the output of MT and upstream power grid. This indicates that the allocation of DRGs in the pre-disaster prevention stage can effectively improve the load restoration capability and reduce the carbon emissions. Compared with the result under random allocation, the optimal allocation of DRGs can further reduce the load loss of 1,807.29 kW·h and the total carbon emission of 1.40% under disaster scenarios. The LER value of scheme 4 is less than that of scheme 3, and it can be concluded that the optimal allocation effect of scheme 4 is better than that of scheme 3.

As can be seen from Table 6, scheme 5 can significantly reduce the load loss caused by disasters and greatly reduce the total cost by optimizing the allocation of DRGs before disasters and adopting NR after disasters. Compared with operation scheme 2, scheme 5 further reduces the total carbon emissions. However, compared with schemes 3 and 4, the total carbon emissions increase. The reason is that NR leads to the restoration of power supply for a large number of loads, and the output of carbon emission units increases. According to the LER value, scheme 5 is the optimal scheme.

In order to demonstrate the effect of NR on restoring lost load, calculation and analysis are carried out under a single disaster scenario. Based on disaster scenario 1, the test is carried out under scheme 1 and scheme 2 respectively. The results are shown in Figure 8 and Table 7. Comparing Figures 8A, B, it can be seen that without NR, two load islands are generated in the distribution system due to line faults. Although the islands are equipped with corresponding black-start power supplies, due to the limited capacity of the unit, most nodes still have load loss, and the reliability of power supply is low. After adopting the network reconstruction method, the switches between the tie line 12–22 and 25–29 is closed, so that the entire distribution system maintains a radial complete connection state, avoiding the generation of load islands. In this case, the load loss of the entire power grid is reduced. Although due to the transmission capacity limit of tie lines, the lost load has not been fully restored, compared to the case where NR is not taken, the load loss has been reduced by 92.04%, and the fault restoration ability of the distribution system has been greatly improved. From the LER value, adopting NR can restore more power at a lower cost of carbon emissions.

In order to compare the effects of different allocation schemes, the calculation under a single disaster scenario is performed based on disaster scenario 4, and the results under three different allocation schemes are obtained, as shown in Figure 9 and Table 8.

As can be seen from Figure 9A, when allocation scheme 1 is adopted, there are four DGs in island 1 for load restoration. However, because there is only one power dividing point of node 7, the PV units located at node 7 and node 26 cannot deliver power to node 27 and its downstream nodes. In fact, the total load of nodes 31 and 32 at the end of the branch is heavy, and the VoLL of the load is large. Therefore, under this allocation scheme, the renewable energy is not fully utilized, and the important load is not sufficiently powered, resulting in more load loss and higher total carbon emissions.

In Figure 9B, compared with scheme 1, scheme 2 changes the allocation position of WT from node 20 to node 24. Because the output of WT is not blocked by the fault line, it can replace part of the output of the upstream power grid. From Table 8, it can be found that the total carbon emission under allocation scheme 2 is lower than that under allocation scheme 1. Although the change in the location of the WT has little impact on island 1, it has caused island 2 to lose the power source for restoration, leading to island 2 becoming a load island. The load loss of the scheme is greatly increased compared with the allocation scheme 1, and the LER value is the largest among the three schemes, which shows that the scheme is not desirable.

The optimal allocation scheme is adopted in Figure 9C, so that there are two power dividing points of node 7 and node 32 in different periods of a day in island 1. It can be seen from the figure that the scheme has more diverse power delivery paths than the first two allocation schemes. During periods when the load of nodes 31 and 32 is small, the PV unit located at node 33 can bypass node 32 to provide power supply to other loads, and at this time node 32 does not act as a power dividing point. When the load at the end of the branch is large, node 32 becomes a power dividing point, and is powered by the power supplies on both sides. This means that under this scheme, PV units are allocated more reasonably and the load power supply mode is more flexible. It can be seen from the results in Table 8 that the optimal allocation scheme has smaller load loss than the first two schemes, and its total carbon emissions are relatively low, and the LER value is the smallest, which is an ideal scheme.

For a given disaster response scheme, changes in disaster severity will lead to differences in disaster response effects. On the basis of scheme 5, the results under the multi-scenario disasters can be obtained by changing the disaster attack budget K b u d g e t , as shown in Figure 10 and Table 9. It can be seen from Figure 10 that as the severity of the disaster increases, the load loss and operating costs continue to increase. From Table 9, it can be seen that the total carbon emissions show a downward trend, indicating that disasters have a certain inhibitory effect on the carbon emissions. However, this inhibitory effect is at the expense of load loss. From Table 9, it can also be found that the LER value is at a high level when the disaster severity is high, which means that the load lost for each unit of carbon emissions generated is large. When the value of K b u d g e t is not greater than 2, the load loss and total cost of the system are relatively low, and the LER value is also small, indicating that the disaster response scheme at this time has a good effect. When the value is greater than 2, the load loss and total cost rise rapidly, and the LER value increases sharply, indicating that the disaster is severe enough that current allocation budget and operational methods have gradually failed to meet the disaster response requirements. Therefore, it is necessary to increase the allocation budget or take more effective measures [such as line hardening (Wang et al., 2019b), etc.] to cope with more severe disaster scenarios. When the value is greater than 4, the load loss and total cost growth slow down, indicating that the marginal benefit of disaster attacks has decreased, and the system is suffering the most severe disaster damage.

According to the impact of the severity of the disaster on the distribution system, it can be divided into three ranges. When K b u d g e t ∈ 0,2 , it is called the inhibition interval, within which the existing disaster prevention and response measures play an effective role in inhibiting the disaster. When K b u d g e t ∈ 2,4 , it was called a growth range, and the existing disaster prevention and response measures within this range gradually failed to meet the disaster response requirements. Therefore, the impact of disaster increased sharply with the severity of the disaster. When K b u d g e t ∈ 4,5 , it is called the destroy range. In this range, the damage degree of the disaster on the distribution system tends to be maximized, and existing disaster prevention and response measures have little effectiveness.

The allocation budget mainly includes the number and capacity of allocated units. On the basis of the scheme 5, taking the value of K b u d g e t as 2, changing the capacity of the allocated units and the number of PV units, the results shown in Figure 11 can be obtained. It can be seen from the figure that with the increase of unit capacity and allocation number, the load loss and carbon emissions of the system have significantly decreased, and the LER value also decreases. It is worth noting that with the increase in unit capacity and number, the downward trend of load loss and LER values becomes slower. This result shows that the marginal benefits of disaster prevention and response brought about by the increase in the allocation budget are decreasing. When making disaster prevention and response decisions in practical, it is necessary to take into account the severity of the disaster and the risk tolerance level of the distribution system, while weighing the marginal benefits and allocation costs, and comprehensively formulate an appropriate allocation scheme.