During extreme natural disasters, the ice load on power lines is considered to gradually increase over time. When the ice layer becomes overloaded, significant bending occurs in the line. Excessive bending deformation can lead to line breakage. The ice overload or breakage of lines results in an increase in unbalanced loads on the poles. If the acting load exceeds the pole’s maximum load-bearing capacity, the crossarms may be damaged, leading to bending and twisting of the pole itself, eventually causing breakage and collapse.

The primary impact of extreme weather on distribution networks is manifested in the significant increase in line fault rates and the occurrence of multiple line faults. Given the multitude of distribution network lines and the vast number of fault scenarios corresponding to different meteorological conditions, each with varying line fault rates, it is feasible to select typical fault scenarios based on their randomness. These scenarios can then provide the line fault rates necessary for research on distribution network fault recovery. Information entropy is a method used to quantify the disorder in a system (Shannon, 1948). The selection of typical fault scenarios based on the fault rates of distribution network lines involves choosing scenarios with a high probability of occurrence. The information entropy is represented as Eq. 9: W = ∑ i ∈ Ω l − ⁡ log 2 ⁡ P j , t fault τ j , t   t ∈ T i c e (9)

Each fault scenario corresponds to an τ j , t vector, representing the entropy value of the system under that scenario. The occurrence probability of a typical fault scenario is relatively high, so its information entropy W should not be too large or too small and must satisfy the Eq. 10: W min ≤ ∑ j ∈ Ω l − ⁡ log 2 ⁡ P j , t fault τ j , t ≤ W max ∑ t τ j , t ≤ 1 t ∈ T i c e (10)

By calculating the entropy values of the distribution network under different scenarios and selecting typical fault scenarios based on reasonable entropy value limits, the fault lines corresponding to these entropy values can be identified. These identified fault lines are then used for the next step in the study of distribution network fault recovery measures.

The fault recovery process for distribution networks under extreme weather scenarios is illustrated in Figure 1.

Step 1: Input meteorological parameters for the extreme weather scenario and relevant data about the distribution network structure.

Step 2: Calculate the distribution of fault rates for each line under the current meteorological parameters using the ice load-based line fault rate model.

Step 3: Generate multiple fault scenarios by changing meteorological parameters, and calculate the corresponding entropy values using the information entropy formula.

Step 4: Determine whether the selected scenario is a typical fault scenario. If it is, select the high fault rate lines in that scenario and formulate a fault recovery strategy for them. if not, return to Step 2.

The objective is to minimize the sum of the load shedding amount, network loss cost, and DG operating cost, which is represented as: f = min x max u ∈ U min h [ λ 1 ∑ i ∈ Ω n 1 − y i w l P i , load + λ 2 ∑ i j ∈ Ω b I i j s q r R i j + ∑ i ∈ Ω n a P i , DG 2 + b P i , DG + c ] (11)

Analyzing Eq. 11, it's evident that the objective function of the proposed robust fault recovery strategy is a min-max-min problem: The first layer min problem uses x as the decision variable, aiming to generate a network topology with the least load shedding. The second layer max problem uses u as the decision variable to find the worst-case fluctuation scenario for DG injection power within the given uncertainty set U . The third layer min problem uses h as the decision variable.

The constraints of the robust fault recovery strategy include power flow equations, system operational safety constraints, SOP operation constraints, and radial network configuration constraints. This paper constructs the distribution network branch power flow equations based on the second-order cone programming model proposed in literature (Farivar and Low, 2013), which include node injected power balance Eq. 12 and branch voltage drop Eq. 13. The specific forms are as follows:

(1) Node Injected Power Balance Equation

During the fault recovery process, changes in the distribution network topology necessitate the relaxation of the branch flow voltage equation using the big-M method. When a branch i j is open, i.e., x i j = 0, the branch flow voltage equation need not be satisfied. Conversely, when a branch i j is closed, i.e., x i j = 1, the branch voltage drop equation must be strictly adhered to. The revised Eq. 13 is as follows: U i s q r − U j s q r ≤ M 1 − x i j + 2 P i j R i j + Q i j X i j − R i j 2 + X i j 2 I i j s q r , ∀ i j ∈ Ω b U i s q r − U j s q r ≥ − M 1 − x i j + 2 P i j R i j + Q i j X i j − R i j 2 + X i j 2 I i j s q r , ∀ i j ∈ Ω b I i j s q r U i s q r ≥ P i j 2 + Q i j 2 , ∀ i j ∈ Ω b (14)

(3) System Operational Safety Constraints

System operational safety constraints include node voltage Eq. 15 and branch capacity Eq. 16. Node voltage constraints ensure that the voltage at each node remains within a safe operating range during the fault recovery period. The constraint is represented as Eq. 15: y i U i min 2 ≤ U i s q r ≤ y i U i max 2 (15)

Branch capacity constraints ensure that the power in each branch does not exceed limits. The constraint is represented as Eq. 16: − P ¯ i j x i j ≤ P i j ≤ P ¯ i j x i j ∀ i j ∈ Ω b − Q ¯ i j x i j ≤ Q i j ≤ Q ¯ i j x i j ∀ i j ∈ Ω b I i j s q r ≤ x i j I ¯ i j s q r (16)

(4) SOP Operating Constraints

Instead of representing DG output with a single deterministic forecast value, a box-type interval model is constructed as Eq. 22: U = P i , DG ∈ P i , DG ref − Δ P i , DG , P i , DG ref + Δ P i , DG ∑ i ∈ Ω n P i , DG − P i , DG ref Δ P i , DG ≤ Γ DG (22) where Γ DG is the “uncertainty parameter” for the possible values of the uncertain quantity, which can be chosen from the set {0, −1, 1}.

To contrast with the robust fault recovery model proposed in this paper, a deterministic fault recovery model is established as follows. The objective function of the deterministic model is to minimize the load shedding amount, network loss, and DG operating costs, specifically formulated as Eq. 23: f = ⁡ min [ λ 1 ∑ i ∈ Ω n 1 − y i w load P i , load + λ 2 ( ∑ i j ∈ Ω b I i j s q r R i j + ∑ i ∈ Ω n a P i , DG 2 + b P i , DG + c ) ] (23)

The constraint conditions, aside from the node injected power balance constraint, are the same as those in the robust fault recovery model. In the deterministic recovery model, the uncertainty of DG output is not considered, and the node injected power balance Eq. 12 is modified to Eq. 24: ∑ i j ∈ Ω b P i j − ∑ k i ∈ Ω b P k i − R k i I k i s q r = P i ∑ i j ∈ Ω b Q i j − ∑ k i ∈ Ω b Q k i − X k i I k i s q r = Q i P i = P i , SOP + P i , DG − y i P i , load Q i = Q i , SOP + Q i , DG − y i Q i , load , ∀ i ∈ Ω n (24)

In summary, the objective function of the robust fault recovery model proposed in this paper is represented by Formula 11, with constraints including node injected power balance Eq. 12, branch voltage drop Eq. 14, node voltage Eq. 15, branch capacity Eq. 16, SOP operating Eqs 17–21, and DG uncertainty set (22). This model is a two-stage robust optimization form of a min-max-min problem with certain constraints and belongs to a large-scale combinatorial optimization problem, which is challenging to solve directly. The robust fault recovery model Eqs 11–22 and the deterministic model Eqs 13–21, 23, 24 provide a comprehensive framework for addressing fault recovery in power distribution networks under uncertain conditions.

To better illustrate the algorithm’s process, the two-stage robust optimization model Eqs 11–22 are integrated into the following compact form: min x max u ∈ U min h ∈ F x , u b T h (25) s . t . A x ≤ b (26) B x = d (27) F x , u = h ∈ S h : C h ≥ E x , D h = u G l h ≤ g l T h , ∀ l = 1 , … , m (28)

Model Eqs 25–28 present a two-stage optimization problem that is challenging to solve directly. To address this, the C&CG algorithm is employed to decompose the two-stage robust optimization model into a master problem and sub-problems, which are then solved iteratively.

The master problem decides the network topology under the condition that the output power of each DG unit, u l * , is already determined. The expression for the master problem is: MP : min x , η η (29) s . t . A x ≤ b (30) η ≥ b T h l , ∀ l ≤ k (31) E x ≤ C l , ∀ l ≤ k (32) D y l = u l * , ∀ l ≤ k (33) G i h l ≤ g i T h l , i = 1 , … , m , ∀ l ≤ k (34)

A real variable η is introduced to replace the max-min problem in Eq. 25, and Eq. 31 is added to obtain a relaxed problem as shown in Eqs 29–34.

The subproblem, solved after determining the network topology in the first phase, finds the worst-case fluctuation scenario of distributed generation u l * . The subproblem is described as Eqs 35–37: SP : Q x * = max u ∈ U min h ∈ F x * , u b T h (35) s . t . C h ≥ E x * (36) D h = u (37) G l x ≤ g l T x , ∀ i = 1 , … , m (38)

The subproblem is a max-min, two-layer optimization problem, which is difficult to solve directly. To address this, dual theory is applied to transform the inner min problem into its corresponding max problem, which is then combined with the outer max problem to form a single-layer max problem. The transformed subproblem expression is: Q x * = max u , π l , π 2 , w i , λ i E x * T π 1 + u T π 2 (39) s . t . C T π 1 + D T π 2 + ∑ i G i w i + g i λ i = b (40) w i 2 ≤ λ i , i = 1 , … , n (41) where π 1 , π 2 are the dual variables for Eqs 36, 37. w i , λ i is the dual variable for Eq. 38.

The objective function of the subproblem (39) contains the nonlinear, bilinear term u T π 2 , which needs to be processed for easier solving of the model. The big-M method can be used to linearize u T π 2 as follows: u T π 2 = ∑ s u s π 2 , s = ∑ s g s g s ≤ π 2 , s g s ≤ M u s g s ≥ π 2 , s − 1 − u s M g s ≥ 0 u s ∈ 0 , 1 (42) where g s is a new variable introduced, which represents the dot product of the corresponding elements of the variables u and π 2 .

The basic principle of robust optimization is to transform the original robust optimization model into a two-stage optimization problem, which is then solved using the C&CG algorithm. The algorithm iteratively solves the master and subproblems, adding the columns generated by the subproblem to the constraints of the master problem. The optimal result is output when both the upper and lower bounds meet the convergence precision. The model-solving process of the C&CG algorithm is illustrated in Figure 2.

The specific iterative process is as follows:

Step 1: Initialize data, set the upper and lower bounds of the objective, let L B = − ∞ , U B = + ∞ , k = 0, and set O as an empty set.

Step 2: Solve the master problem Eqs 29–34, obtaining the optimal solution x k + 1 * , η k + 1 * , h 1 * , . . . , h k * , and update the lower bound L B = η k + 1 * of the original problem.

Step 3: Substitute the given x = x k + 1 * from the master problem into the subproblem Eqs 39–42, solve the worst-case scenario u k + 1 * and the optimal objective function value Q ( x k + 1 * ) , and update the upper bound U B = ⁡ min U B , Q x k + 1 * of the original problem.

Step 4: Check if the convergence condition UB − LB ≤ ε is met. If yes, the iteration is complete, and the current optimal solution is output. If not, proceed as follows:

(1) If Q x k + 1 * < + ∞ , create a new decision variable h k + 1 , add the Eq. 43 to the master problem Eqs 29–34:

Update the constraint set accordingly. Set k = k + 1 and return to Step 2.

In the simulated scenario, the distribution network is under snowstorm conditions with an ambient temperature of −5°C, wind speed of 10 km/h, and precipitation rate of 40 mm/h. The center of the snowstorm, located at coordinates (−150 km, −150 km), moves towards the distribution network at a speed of 4.2 km/h in the 45° direction along the x-axis. The climate impact radius R is set at 130 km, with R values of 0.4 for line ice coverage fault rate, 0.18 for line fault rate, and 0.1 for pole ice coverage fault rate, increasing to 0.46 and 0.43, respectively, as conditions worsen. Over time, both the ice thickness and the line-pole fault rates change. Figure 4 illustrates the ice thickness and line-pole fault rates for this fault scenario.

Figure 4A shows that as the snowstorm center moves closer to the distribution network, the ice layer thickness on the lines rapidly increases. Figure 4B indicates that with the increase in ice layer thickness, the fault rates of lines and poles first experience a slow growth, then the line fault rate exhibits an exponential increase after 20 h, and the pole fault rate rises rapidly after 30 h. Line breakage may occur if the ice thickness exceeds the line’s load-bearing capacity.

By varying environmental temperature, precipitation rate, snowstorm center location, and movement speed, various extreme ice and snow disaster scenarios can be simulated. The entropy value corresponding to each scenario can be calculated using Eq. 9, with the probability distribution of the system’s information entropy shown in Figure 5. The higher the probability of occurrence of an entropy value, the more likely the scenario is considered a typical fault scenario.

For this simulation, the fault rates of each line in Distribution Network 1 are shown in Figure 6. The entropy value calculated for this scenario using Eq. 9 is 14.71. According to Figure 5, this entropy value has a high probability of occurrence, marking it as a typical extreme ice and snow fault scenario. Figure 6 reveals that Lines 1 (a1-a2) and 18 (a2-a19) in Distribution Network 1 have higher fault rates than other lines, indicating that these lines are relatively vulnerable within the system and more prone to faults due to ice and snow disasters. Consequently, the fault recovery strategies proposed in this paper will be applied specifically to these two high-fault-rate lines to formulate fault recovery plans.

For Distribution Network 1, Lines 1 and 18 were selected as fault lines, leading to the cessation of all loads from nodes a2 to a33, resulting in a substantial active power loss of up to 3715 kW. To analyze the impact of SOP integration on fault recovery, three SOP integration scenarios were tested:

Scenario 1: Fault recovery through internal SOP1 and network reconfiguration in Distribution Network 1. Scenario 2: Based on Scenario 1, connecting Distribution Network 2 through SOP2. Scenario 3: Based on Scenario 1, connecting Distribution Networks 2 and 3 through SOP2 and SOP3, respectively.

The fault recovery effects for these three scenarios are shown in Table 2.

Table 2 reveals that Scenario 1 has the lowest total load recovery rate at 31.89%. This is due to the disconnection of Distribution Network 1 from the main power source, with the required active power for loads being supplied only by internal DG and transmitted through internal SOP1, which also provides reactive power compensation. The fault recovery capacity in this scenario is limited. Comparing Scenarios 1 and 2, the total load recovery rate in Scenario 2 increases to 45.09%, indicating that connecting Distribution Network 1 with Network 2 through SOP2 allows flexible power transmission control, restoring more lost loads. Scenario 3 has the highest total load recovery rate at 55.31%, as the faulted distribution network is interconnected with Networks B and C through two SOPs, allowing two power sources to support power restoration in the outage area. However, not all loads are restored under this SOP capacity, as the capacity of SOP influences its power flow transfer ability and consequently affects the fault recovery outcome.

To analyze the impact of SOP capacity, Scenario 3 was tested under four different SOP capacities: 200, 400, 600, and 800 kVA, and the fault recovery conditions and SOP output powers were assessed. The influence of SOP capacity on fault recovery effectiveness is shown in Table 3.

As indicated in Table 3, as the integrated SOP capacity in the system increases, the fault recovery capability of the distribution network continuously strengthens. With a single SOP capacity of 200 kVA, the total load recovery rate is relatively low at 55.31%, with a primary load recovery rate of only 63.75%, indicating poor recovery effectiveness. At 400 kVA, the total load recovery rate improves to 72.01%, with primary loads fully restored, and a slight increase in the tertiary load recovery rate. At 600 kVA, nodes a7 and a25 are restored, with a high total load recovery rate of 88.69%, including 100% recovery for primary and secondary critical loads. With 800 kVA, the total SOP capacity accounts for 64.60% of the load power, with strong power flow transfer capability, resulting in the restoration of all lost loads.

The active power outputs of SOPs under four different capacities are shown in Figure 7.

Figure 7 shows that for all four SOP capacities, the active power of VSC a12 in Distribution Network 1 is negative, while VSC a22’s is positive, indicating that the active power flow on the SOP is from node a12 to node a22. After Distribution Network 1 is disconnected from the main power source, VSC a12 obtains electrical energy from nearby DGs at nodes a9 and a14. VSC a22, lacking nearby DG resources, receives active power transmitted through SOP1 from node a12 to restore the lost loads 19–22. For SOP2 connecting Distribution Networks 1 and 2, the active power of VSC a18 is positive, and that of VSC b18 is negative, indicating active power flow from Network 2’s node b18 to Network 1’s node a18. Similarly, SOP3 injects active power from Distribution Network 3 into Network A’s node a28 to provide electricity for the lost loads.

Considering that the integration of DG changes the system’s power flow distribution and affects the fault recovery outcome, this section analyzes the impact of DG integration on fault recovery. Based on the integration of three 200 kVA SOPs, a comparative analysis was conducted with different levels of DG injection power. The fault recovery effects under varying DG injection power levels are presented in Table 4.

Table 4 shows that with the increase in DG injection power, the total load recovery rate correspondingly rises from 55.31% to 100.0%. When the DG injection power is 150 kW, only 63.75% of primary load can be restored, indicating low power supply reliability. At a DG injection power of 225 kW, the recovery rate of primary load increases from 63.75% to 100%, and the recovery rates of secondary and tertiary loads also rise. With 300 kW of DG injection power, primary and secondary critical loads are fully restored, and the tertiary load recovery rate increases to 73.10%. At 375 kW, DG within the Distribution Network 1 accounts for 40.38% of the total load power, enabling the recovery of all loads. Thus, it is evident that DG injection power plays a positive role in fault recovery. DG injection power can locally meet the demands of nearby loads, enabling more loads to be restored.

To demonstrate the effectiveness of the robust fault recovery strategy for distribution networks proposed in this paper, a comparison of fault recovery effects between deterministic and robust strategies was conducted. In the deterministic recovery strategy, DG power remains at the forecasted value, while in the robust strategy, DG output fluctuates within the range of an uncertainty set. With a single SOP capacity of 200 kVA and the same parameter settings, the fault recovery effects of the two strategies are shown in Table 5.

As evident from Table 5, the recovery effect of the robust strategy is superior to that of the deterministic strategy. The deterministic strategy, which bases the network topology on the predicted values of DG, has weaker adaptability to uncertainties, resulting in a lower total load recovery rate of only 55.31%. In contrast, the robust strategy, considering the worst-case DG output samples, achieves an improved total load recovery rate of 61.78%. Furthermore, it ensures 100% recovery of primary loads, though the recovery rate of secondary loads is slightly reduced. It is apparent that the proposed robust strategy yields better recovery results than the deterministic strategy and possesses a superior ability to adapt to the uncertainties in DG output.

To verify the scalability of the proposed method in large-scale actual systems, an analysis was conducted using a 234-node system interconnected by SOPs. This system consists of three 78-node systems interconnected through SOPs (Ji et al., 2023). In Distribution Network A, 7 DGs were connected at nodes 5, 24, 29, 47, 52, 60, and 78, with each having a forecasted active power of 150 kW and a power factor set to 0.9. SOP1 connects to Distribution Network A at nodes a22 and a59, SOP2 connects Distribution Network A at node a18 and Distribution Network B at node b18, and SOP3 connects Distribution Network A at node a28 and Distribution Network C at node c28, with each having a capacity set to 1000kVA and a loss coefficient of 0.02. Two fault recovery strategies, deterministic recovery strategy and robust recovery strategy, were compared. Table 6 presents the recovery effects of the two fault recovery strategies under the same fault location and parameter settings.

From Table 6, it is evident that the robust strategy still outperforms the deterministic strategy in terms of load recovery effectiveness. The deterministic strategy exhibits a weaker adaptability to the uncertainty of DGs, with a total load recovery rate of 76.08%. In contrast, the robust strategy achieves an improved total load recovery rate of 81.38%, with critical loads being fully restored to 100%. Thus, in large-scale actual systems, the proposed robust strategy can achieve better recovery outcomes than the deterministic strategy.