Restoration Strategy for Urban Power Distribution Systems Considering Coupling with Transportation Networks Under Heavy Rainstorm Disasters

Abstract

With the increasing frequency of extreme weather events such as heavy rainstorm disasters, the stable operation of power systems is facing significant challenges. This paper proposes a two-stage restoration strategy for the distribution networks (DNs). First, a grid-based modeling approach is developed for urban DNs and transportation networks (TNs), capturing the dynamic evolution of heavy rainstorm disasters and more accurately modeling the impact on TNs and DNs. Then, a two-stage restoration strategy is designed for the DN by coordinating soft open points (SOPs) and mobile energy storage systems (MESSs). In the disaster progression stage, SOPs are utilized to enable the flexible reconfiguration and islanding of the DN, minimizing load loss. In the post-disaster recovery stage, the MESS and repair crew are optimally dispatched, taking into account the state of the TN to expedite power restoration. Finally, the experimental results demonstrate that the proposed method reduces load loss during restoration by 8.09% compared to approaches without precise TN and DN modeling.

1. Introduction

As global climate change intensifies, the increasing frequency of extreme weather events is exerting growing pressure on the stability and reliability of power and transportation systems [1,2,3]. In recent years, urban flooding caused by heavy rainstorm disasters has resulted in severe damage to power and transportation infrastructure [4]. For instance, in August 2023, Zhuozhou City, Hebei Province, China, was hit by a heavy rainstorm disaster, resulting in the shutdown of 1500 kV line, 3110 kV substations, 335 kV substations, and 6610 kV lines due to the disaster or emergency evacuation, affecting the safe and reliable use of electricity in 97 neighborhoods [5]. In February 2024, California was hit by a heavy rainstorm disaster, leading to road closures and power outages that affected over 20,0000 users [6]. In May 2024, a heavy rainstorm disaster in Guangxi, China, caused six 10 kV power lines to go offline, impacting the electricity supply of 4970 households [7].

As critical hubs for public systems, underground distribution equipment faces heightened risks during heavy rainstorm disasters, especially due to aging drainage systems and surface runoff [8]. During urban flooding, water accumulation can obstruct traffic, significantly impeding the mobility of power grid repair teams [9]. This may lead to extensive load loss, grid separation, or even system collapse. Therefore, it is imperative to accurately model the progression of heavy rainstorm disasters and develop feasible and effective restoration strategies to mitigate their impact.

In recent years, significant progress has been made in studying severe rainstorm disasters and their destructive effects on TNs and DNs. Heavy rainstorm disasters not only disrupt road traffic and reduce transportation flow but also severely damage power systems, thereby affecting normal urban operations. In-depth research on coupled damage and DN restoration, considering coupling with TNs during heavy rainstorm disasters, is crucial.

A rainfall pattern construction method based on dynamic time warping (DTW) is proposed in [10]. This method comprehensively recreates common urban rainfall processes by simulating spatiotemporal rainfall distribution characteristics. Rainfall development under various return periods and peak conditions is simulated in [11] using the classical Chicago rainfall pattern method, establishing foundations for urban water-logging mechanisms during heavy rainstorm disaster. The study conducted in [12] further investigated urban waterlogging processes under storm conditions and proposed a rainfall-based simulation method to effectively predict water accumulation patterns in urban low-lying areas.

The damage mechanisms of heavy rainstorm disasters on DN have received widespread attention. The study in [13] systematically analyzed the impacts of floods on power systems and proposed disaster response measures. Substation submergence caused by flooding frequently results in large-scale blackouts and significant economic losses. The study in [14] introduced a GIS-based probabilistic model to quantify power shortages and economic losses caused by substation flooding. These findings offer a theoretical basis and modeling tools for designing more resilient power system layouts.

Regarding TNs [15,16] examined how water depth reduces traffic flow and how road disruptions impact overall network efficiency. Research indicates that TN failures not only hinder human mobility but also cause cascading effects on other critical infrastructure, such as power systems. However, studies on the coupled failure relationship between TNs and power DNs during heavy rainstorm events remain limited.

After extreme disasters, utilizing transportation system resources for the timely restoration of the distribution network is an effective measure. Recently, distributed energy resources (DERs), including MESSs and SOPs, have introduced new restoration methods. The study in [17] highlights that using distributed energy resources (DERs), such as MESSs, can significantly enhance post-disaster restoration efficiency in DNs. In [18], a multi-period optimization model was proposed, considering network reconfiguration and various DERs. A coordinated operation method considering the interaction between power and transportation systems is proposed in [19]. Reference [20] develops a two-stage reconfiguration strategy based on the contribution of SOPs. The research in [21] introduces virtual load and islanding models to improve reliability cost assessment during restoration. The work in [22] suggests utilizing electric vehicles (EVs) as flexible mobile resources for load restoration.

When developing restoration strategies, it is essential to consider the characteristics of the disaster itself in order to propose targeted methods. The study in [23] mapped water depth to corresponding locations during the evolution of water inundation disasters and proposed improved strategies for enhancing network resilience. The research in [24] introduced disaster insurance as a complementary approach for risk sharing and resilience enhancement. Reference [25] considers the coupling between power distribution and heating systems and proposes a network reconfiguration method for the coordinated operation of the two systems. The authors of [26,27] proposed a two-stage restoration strategy, including pre-disaster resource allocation and post-disaster scheduling. Further, the work in [28] incorporated repair crews and developed a fault restoration model that considers road damage and the collaboration of dual repair teams.

A deeper consideration of the interdependence between DN and TN faults during the progression and evolution of disasters could yield significant insights for enhancing distribution network restoration strategies. This paper proposes a restoration strategy for DNs that considers coupling with TNs under heavy rainstorm disasters. Based on rainfall intensity’s compound effects on TNs and DNs, we present an indirect coupling model considering precipitation intensity and design a two-stage restoration method for DN reconfiguration and emergency resource dispatch to maximize restoration efficiency and system stability. The specific contributions of this study are as follows:

(1)

A failure model of DN and TN is developed based on urban grid partitioning under heavy rainstorm disasters. Using a two-dimensional surface hydrodynamic model, flooding scenarios are constructed to simulate the impacts of heavy rainstorm disasters on DN equipment and TN infrastructure. This approach quantifies the failure mechanisms of TNs and DNs, correlating rainfall intensity with infrastructure damage. It provides decision support for disaster restoration efforts.

(2)

A restoration model is developed to integrate DN reconfiguration with emergency resource scheduling. During the disaster progression stage, SOP technology is applied to enable flexible DN reconfiguration, minimizing load loss. In the restoration stage, the framework considers constraints imposed by TN disruptions and DN restoration needs. This model is validated in IEEE33-bus DN and Sioux-Falls24-nodeTN.

2. A Failure Model of DNs and TNs Under Heavy Rain Disasters

2.1. Framework of DN and TN Failure Models

This model primarily employs a grid partitioning method to physically model urban areas. Based on the dynamic evolution of flooding depth caused by heavy rain disasters, a system encompassing both DNs and TNs is established, as shown in Figure 1.

By mapping urban areas to the corresponding distribution and TNs through grid division, it can be seen that areas with significant flooding depth experience severe damage to the DN and greater load shedding during the disaster progression stage. In the restoration stage, more mobile resources are required for restoration in these areas. However, transportation roads in the same areas are also heavily impacted by flooding, reflecting the indirect coupling relationship between DNs and TNs under heavy rain disasters.

2.2. Urban Flooding Model Under Heavy Rain Disasters

To facilitate the modeling of urban flooding damage caused by heavy rain disasters, this paper adopts a grid partitioning approach. A city area is physically divided into several grids, each corresponding to a portion of the DN and TN in the real world. The following assumptions are also made:

• Uniform rainfall distribution within grids: due to the relatively small scale of grid partitioning compared to the rainfall range, it is assumed that the rainfall amount is equal within each grid.
• Uniform flood depth within grids: when water flows due to flooding, each grid is treated as a single unit with an equal depth of accumulated water.
• Average flood depth at grid boundaries: for roads in the TN that span grid edges, the flood depth is taken as the average of the flood depths of the two adjacent grids; since DN failures primarily relate to bus failures, there are no cases spanning grid edges.

Heavy rain typically features high rainfall intensity and distinct single peaks. The predominant international method for setting heavy rain patterns is the Chicago Rain Model [29]. This method, based on heavy rain intensity formulas, flexibly adjusts the peak position and total rainfall duration according to statistical data to accurately describe the rainfall process. The pre-peak rainfall intensity is given as follows:

$$I(t_a)={\displaystyle \frac{M\left[\frac{(1-n)t_a}{k}+l\right]}{{\left(\frac{t_a}{k}+l\right)}^{n+1}}}$$

(1)

where $t_a$ is the pre-peak rainfall time; $M$ is the rainfall intensity parameter, which is associated with the rainfall amount for a specific return period; $l$ is the correction parameter for rainfall duration; $k$ is the rainfall peak position coefficient, ranging from 0 to 1, and this used to indicate the position of the rainfall peak during the entire rainfall process; and $n$ is the decay coefficient of the heavy rain and is related to the return period.

The post-peak rainfall intensity is expressed as follows:

$$I(t_b)={\displaystyle \frac{M\left[\frac{\left(1-n\right)t_b}{1-k}+l\right]}{{\left(\frac{t_b}{1-k}+l\right)}^{n+1}}}$$

(2)

where $t_b$ is the post-peak rainfall time.

By integrating the rainfall intensity, the depth of accumulated water caused by rainfall can be calculated. However, in the real world, the depth of accumulated water is also influenced by water flow and urban drainage. This paper employs a two-dimensional surface hydrodynamic model to simulate the water flow process [30], considering the balance between rainfall amount, drainage capacity, and water flow rate to determine the variation in the accumulated water depth in each region.

The two-dimensional surface model mentioned in this paper assumes that accumulated water is in vertical force balance and is subjected to frictional force $f$ and the pressure from surrounding water $F_N$ in the horizontal direction, which are given as follows:

$$F_N=-mg\nabla h$$

(3)

$$f=-mg{k}_f^2{n}_m^2{d}^{-4/3}{\left|\mathit{v}\right|}^2{\displaystyle \frac{\mathit{v}}{\left|\mathit{v}\right|}}$$

(4)

where m is the mass of the accumulated water, g is the gravitational acceleration, h is the water surface elevation, $k_f$ is the unit conversion coefficient, $n_m$ is the Manning coefficient, d is the depth of accumulated water, and $\mathit{v}$ is the horizontal flow velocity vector.

Considering the initial velocity of the accumulated water to be zero, working based on the conservation of momentum, the following equation can be obtained:

$$(\mathit{f}+{\mathit{F}}_{\mathrm{N}})\Delta t=m\mathit{v}$$

(5)

where $\Delta t$ is the time step. By dividing the accumulated water area into small micro-elements and substituting them into the equation above, a univariate quadratic equation in v is derived, allowing the calculation of the horizontal velocity of the accumulated water. Integrating the velocity demonstrating the impact of water flow on the depth.

Additionally, urban drainage and infiltration are considered, with the infiltration volume of accumulated water defined as follows:

$$Q_{\mathrm{d}}=\mu S\sqrt{2gd}$$

(6)

where $\mu$ is the discharge coefficient, and S is the equivalent infiltration area.

Furthermore, the depth of accumulated water is related to the building coverage rate.

In order to facilitate the solution, it is necessary to assume that only the flow of water in the four principal directions of east, west, south and north is considered and that the change in the depth of standing water per unit of time is as follows:

$$\Delta d={\displaystyle \frac{1}{1-\theta }}\left({\displaystyle \frac{Q_{\mathrm{W}}+Q_{\mathrm{S}}+Q_{\mathrm{E}}+Q_{\mathrm{N}}}{\Delta x\Delta y}}-Q_d\Delta t+I(t)\Delta t\right)$$

(7)

where $Q_{\mathrm{W}}$, $Q_{\mathrm{S}}$, $Q_{\mathrm{E}}$, and $Q_{\mathrm{N}}$ represent the flow rates in the east, west, south, and north directions, respectively.

2.3. TN Failure Model Under Heavy Rain Disasters

For TNs, heavy rain-induced flooding often causes a decrease in vehicle speeds on roads. When the flood depth exceeds a certain threshold, roads may become impassable, altering the network’s topology. Considering the impact of flood depth on the TN, an affected TN model is established based on simulated flood depth data [31]. This model primarily analyzes the impact of flood depth on road capacity, especially in scenarios where flooding may lead to road surface blockages, reduced vehicle speeds, and traffic interruptions.

Specifically, when the flood depth on a road exceeds the critical flood depth $d_p$ that vehicles can traverse, the travel speed $V$ is reduced to zero, indicating that vehicles cannot pass through that road segment. When the flood depth is below the critical depth, the travel speed $V$ can be calculated using the following equations:

$$V=\left\{\begin{array}{ll}{\displaystyle \frac{V_0}{2}}\tanh \left(-{\displaystyle \frac{d+0.5d_p}{3}}\right)+{\displaystyle \frac{V_0}{2}}/3.6,\hfill & d<d_p\hfill \\ 0,\hfill & d\ge d_p\hfill \end{array}\right.$$

(8)

where $V_0$ is the designed travel speed under no flooding conditions, $d_p$ is the critical flood depth that vehicles can traverse, $V$ represents the travel speed of emergency vehicles at simulation time step t, and d is the flood depth on the road segment at that simulation time step.

2.4. DN Failure Model Under Heavy Rain Disasters

For the DN, heavy rain disasters typically cause two types of damage. On one hand, urban flooding can submerge distribution equipment such as distribution cabinets, leading to bus failures in the DN. On the other hand, heavy rain disasters are often accompanied by hurricanes, which can damage poles, lines, and other facilities. This study focuses on urban DNs; therefore, only the former type of damage is primarily considered.

Specifically, based on the model in Section 2.1, the accumulated water depth is calculated at a particular time. When the accumulated water depth is below the equipment’s designed water depth, it is assumed that the flood does not cause equipment failure, resulting in a failure probability of zero. When the accumulated water depth exceeds the equipment’s designed water depth, the failure probability of the equipment increases significantly and is exponentially related to the flood depth. When the accumulated water depth surpasses the distribution station’s flood resistance depth H_w, the failure probability becomes one. Working based on the dynamic variation in flood depth, the failure rate of equipment at time t can be defined by the following equations [32]:

$$\lambda (t)=\left\{\begin{array}{ll}0,\hfill & h_w(t)\le H_{\mathrm{bw}}\hfill \\ \epsilon \exp \left(\gamma {\displaystyle \frac{h(t)-H_{\mathrm{bw}}}{H_w}}\right),\hfill & H_{\mathrm{bw}}<h(t)\le H_w\hfill \\ 1,\hfill & h(t)>H_w\hfill \end{array}\right.$$

(9)

where h is the accumulated water depth in the region where the equipment is located at time t, $H_w$ is the design flood resistance depth of the distribution station, $H_{\mathrm{bw}}$ is the critical flood depth below which the equipment remains undamaged, $\epsilon$ is the reduction coefficient representing the initial failure rate, and $\gamma$ is the scale parameter used to adjust the rate at which the failure rate increases.

Using the cumulative failure probability formula, the failure probability of the equipment within a specific time period can be determined.

Assuming that the failure probabilities of all critical equipment within the station are independent, the failure probability of the distribution station within a time interval $\Delta t$ can be derived as follows:

$$P(t)=1-\exp [-{\displaystyle \sum _{i\in I}\lambda (t)\Delta t}]$$

(10)

This failure probability is subsequently used to generate disaster scenarios for DNs under heavy rainstorm disasters. These stem from the failure of distribution stations.

3. Two-Stage Disaster Restoration Model

3.1. Model Framework

The disaster progression is illustrated in Figure 2, where it is divided into the pre-disaster prevention period, the disaster progression period, and the post-disaster restoration period. This paper mainly focuses on the defense and restoration efforts after the disaster.

Firstly, based on the probabilistic failure model of the DN outlined in Section 2.4, it is assumed that the failure probabilities of buses are independent during heavy rain periods. A state sampling method is employed to extract the states of each bus, thereby enabling scenario generation for load shedding in the DN.

During the disaster progression stage, the state of the DN is continuously changing in real time and is adversely affected by heavy rain. Under such circumstances, mobile energy storage systems (MESSs) find it challenging to support the DN effectively. Therefore, during the disaster progression stage, intelligent SOPs within the DN are utilized to perform islanding and reconfiguration through the switching of SOPs, thereby reducing load loss during disaster occurrence.

In the post-disaster restoration stage, as the intensity of rainfall decreases, mobile resources such as MESSs and mobile power generation vehicles are deployed to appropriate buses within the DN to replace distributed generators (DGs). SOPs are again used to perform islanding and reconfiguration to restore the lost load.

3.2. DN Reconfiguration Under SOP Commitment

In this stage, SOPs are primarily used to automatically detect and isolate fault areas, adjust the power flow, and delineate islanded supply regions. This process aims to maximize the stable operation of the DN while minimizing operational costs. Specifically, the approach involves the following steps:

The goal is to minimize the amount of load shedding, which is formulated as follows:

$$\min C_1=c_{\mathrm{load}}{\displaystyle \sum _{t=1}^{T_1}{\displaystyle \sum _{i=1}^{N_{\mathrm{bus}}}\Delta }}t{P}_{i,\mathrm{LD}}^t$$

(11)

where $c_{\mathrm{load}}$ is the cost per unit of load shedding, $N_{\mathrm{bus}}$ is the set of all bus, $T_1$ is the first stage of the restoration time under heavy rainstorm disasters, $\Delta t$ is the length of the time interval, and $P_{i,\mathrm{LD}}^t$ is the active power of the load at bus i that can be shed at time t.

The constraints primarily include component constraints, topology constraints, and power flow constraints.

Component constraints include constraints on the loads and energy storage systems. Loads are categorized into controllable and uncontrollable loads. During the disaster, a portion of the controllable loads can be shed to achieve restoration. The constraints on energy storage systems are expressed as capacity constraints and charge/discharge power constraints:

$$\left\{\begin{array}{l}0\le P_{j,\mathrm{LDU}}^t\le P_{j,LU}^t,\forall t\in [1,T]\hfill \\ 0\le Q_{j,\mathrm{LDU}}^t\le Q_{j,LU}^t,\forall t\in [1,T]\hfill \end{array}\right.$$

(12)

$$\left\{\begin{array}{l}{P}_{j,\mathrm{LD}}^t=z^t{P}_{j,L}^t,\forall t\in [1,T]\hfill \\ Q_{j,\mathrm{LD}}^t=z^t{Q}_{j,L}^t,\forall t\in [1,T]\hfill \end{array}\right.$$

(13)

$$\left\{\begin{array}{l}{E}_{j,\mathrm{bat}}^{t+1}=E_{j,\mathrm{bat}}^t+P_{j,\mathrm{ch}}^t{\eta }^{\mathrm{ch}}\Delta T-P_{j,\mathrm{dis}}^t\Delta T/{\eta }^{\mathrm{dis}}\hfill \\ E_{j,\mathrm{bat},\max }\times 20\%\le E_{j,\mathrm{bat}}^{t+1}\le E_{j,\mathrm{bat},\max }\times 90\%\hfill \\ E_{j,\mathrm{bat},\max }\times 20\%\le E_{j,\mathrm{bat}}^t\le E_{j,\mathrm{bat},\max }\times 90\%\hfill \end{array}\right.$$

(14)

$$\left\{\begin{array}{l}0\le P_{j,\mathrm{dis}}^t\le D_{j,\mathrm{dis}}^t{P}_{j,\mathrm{dis},\max }\hfill \\ 0\le P_{j,\mathrm{ch}}^t\le D_{j,\mathrm{ch}}^t{P}_{j,\mathrm{ch},\max }\hfill \\ D_{j,\mathrm{dis}}^t+D_{j,\mathrm{ch}}^t\le 1,\forall t\in [1,T],\forall j\in [1,N_{\mathrm{bus}}]\hfill \end{array}\right.$$

(15)

where (12) represents the capacity constraints for uncontrollable loads, while (13) uses a binary variable z to indicate the state of controllable loads. $P_{j,\mathrm{LDU}}^t$ and $Q_{j,\mathrm{LDU}}^t$ denote the current active and reactive power of the uncontrollable loads, and $P_{j,LU}^t$ and $Q_{j,LU}^t$ represent the total active and reactive power of the uncontrollable loads. $P_{j,\mathrm{LD}}^t$ and $Q_{j,\mathrm{LD}}^t$ denote the current active and reactive power of the controllable loads, and $P_{j,L}^t$ and $Q_{j,L}^t$ represent the total active and reactive power of the controllable loads.

Above, (14) and (15) define the capacity constraints for energy storage systems, where $P_{j,\mathrm{ch}}^t$ and $P_{j,\mathrm{dis}}^t$ denote the charge and discharge power of the energy storage at bus j, η is the charge/discharge efficiency, $E_{j,\mathrm{bat}}^t$ is the energy stored at bus j at time t, and $E_{j,\mathrm{bat},\max }$ is the maximum energy capacity of the storage system. Equation (16) imposes the charge/discharge power constraints, where $P_{j,\mathrm{dis},\max }$ and $P_{j,\mathrm{ch},\max }$ represent the maximum discharge and charge power of the storage system, and $D_{j,\mathrm{ch}}^t$ and $D_{j,\mathrm{dis}}^t$ indicates the charge/discharge state of the storage system.

Topology constraints include virtual power flow constraints and other constraints to ensure proper islanding and reconfiguration. The virtual power flow constraints [25] have been previously introduced and will not be reiterated here. Other constraints ensure the delineation of different islanded regions, maintaining the division and merging of islands across different time periods to optimize the number of segmented regions. These are as follows:

$${\displaystyle \sum _{q=1}^{N_C}{M}_{i,q}^t}=1,\hspace{1em}\forall i\in [1,N_{\mathrm{bus}}],t\in [1,T]$$

(16)

$${\displaystyle \sum _{i=1}^{N_{\mathrm{bus}}}{M}_{i,q}^t}\ne 1,\hspace{1em}\forall q\in [1,N_C],t\in [1,T]$$

(17)

$$M_{i,q}^t=1,\hspace{1em}\forall N_q^t=1,t\in [1,T]$$

(18)

where $M_{i,q}^t$ is a binary variable indicating whether bus i is within island q at moment t, and $N_q^t$ is a binary variable indicating whether island q exists at moment t. Equation (16) ensures that each bus belongs to only one island, Equation (17) ensures that there are no isolated bus within an island, and Equation (18) ensures that if island q exists at time t, there must be at least one bus i that belongs to island q.

The DistFlow power flow constraints are employed by power flow constraints. In addition to satisfying voltage and current constraints, the binary variable ${\alpha }_{i,j}^t$ is introduced to represent the disconnection of lines due to bus failures. For faulted lines, the constraints primarily include the following:

$$\left\{\begin{array}{l}-{\alpha }_{i,j}^{t}C\le P_{i,j}^t\le {\alpha }_{i,j}^{t}C\hfill \\ -{\alpha }_{i,j}^{t}C\le Q_{i,j}^t\le {\alpha }_{i,j}^{t}C\hfill \\ {(P_{i,j}^t)}^2+{(Q_{i,j}^t)}^2\le S_{{\max }^2}\hfill \end{array}\right.\hspace{1em}\forall (i,j),\forall t\in [1,T]$$

(19)

where ${\alpha }_{i,j}^t$ represents the switch state of the line between bus i and j at time t, C is a large constant, $P_{i,j}^t$ and $Q_{i,j}^t$ are the active and reactive power flowing into bus j from bus i at time t, and $S_{\max }$ is the maximum apparent power of the line.

3.3. Multi-Period MESS Scheduling

As the heavy rain disaster transitions to the fault restoration stage, the intensity of rainfall diminishes. MESS and repair teams, acting as mobile emergency resources, are deployed to restore power supply to the DN. This stage primarily considers constraints related to the TN and repair teams. The reconfiguration results from the first stage are used as inputs for the second stage, where additional constraints for MESSs are incorporated to optimize their routing.

The goal is to maximize the restored load during the restoration stage. This is formulated as follows:

$$\min C_2=c_{\mathrm{load}}{\displaystyle \sum _{t=1}^{T_2}{\displaystyle \sum _{i=1}^{N_{\mathrm{bus}}}\Delta }}t{P}_i^t$$

(20)

where $c_{\mathrm{load}}$ is the cost per unit of load shedding, $N_{\mathrm{bus}}$ is the set of all bus, $T_2$ is the second stage of the restoration time under heavy rainstorm disasters, $\Delta t$ is the length of the time interval, and $P_{i,\mathrm{LD}}^t$ is the active power of the load at bus i that can be shed at time t.

The constraints primarily include TN constraints and DN operation constraints.

TN constraints:

$${\displaystyle \sum _i{x}_{i,j}^k}\le 1,j\in i_{MESS}$$

(21)

$${\displaystyle \sum _{i,j,k}{x}_{i,j}^k}\le N_{MESS}$$

(22)

where $x_{i,j}^k$ is a binary variable that indicates whether the k-th MESS has transmission from node i to node j, $N_{MESS}$ denotes the total number of MESSs, and $i_{MESS}$ denotes the node where the MESS is deployed. Equation (21) indicates that the MESSs must be from the deployed node where it is located, and (22) indicates that the MESS number must be within the number of deployed MESSs.

$${\mu }_i^k+t{r}_{i,j}^k+-M\left(1-x_{i,j}^k\right)\le {\mu }_j^k$$

(23)

$$p_j{x}_{i,j}^k\Delta t\le q_j^k\le Q_{\max }^k-r{L}_{i,j,h}+M\left(1-x_{i,j}^k\right)$$

(24)

$$Q_{\min }^k\le q_i^k\le Q_{\max }^k$$

(25)

where ${\mu }_j^k$ denotes the time at which the k-th MESS arrives at node j, and $t{r}_{i,j}^k$ denotes the travel time from node i to node j. Equation (23) describes the logical relationship using a large M. If the MESS travels to node j after completing the service at node i, the time of arrival at node j is at least the time at which it started the service at node i, plus the travel time from i to j. Otherwise, the constraint is relaxed by a sufficiently large positive number M. $q_j^k$ denotes the actual power of the k-th MESS when it arrives at node j, $p_i$ is the power restoration demand at node i, r is the power required to move the emergency energy storage device per unit path, and $Q_{\min }^k$ and $Q_{\max }^k$ denote the maximum capacity of the k-th MESS. Equation (24) ensures that the power of the MESS when it arrives at node j is not less than the demand at node j. The right-hand side expresses this logical relationship by introducing a large number of M. That is, if the MESS travels to node j after completing the service at node i, the power at the time of arrival at node j is the power after traveling from node i to node j minus the demand at node i and the power loss during travel. Otherwise, this constraint is relaxed by a sufficiently large positive number M. Equation (25) then ensures that the amount of power when MESSs arrive at any node is within the specified upper and lower bounds.

Repair constraints:

$${\alpha }_{i,j}^{t_0}\le {\alpha }_{i,j}^t,\forall (i,j)\in B,t_0\le t$$

(26)

$${\displaystyle \sum _B{\alpha }_{i,j}^{t+\Delta t}}-{\displaystyle \sum _B{\alpha }_{i,j}^t}\le m,\forall t\in [1,N_T],\forall (i,j)\in B$$

(27)

where B is the set of nodes where faulted lines are located. Equation (26) ensures that the status of the lines can only be maintained or repaired after the repair begins, and cannot be further damaged, while (27) limits the number of lines that can be repaired simultaneously to m.

The variable $t{r}_{i,j}^k$ is calculated using Dijkstra’s Algorithm 1 that considers the effect of actual road conditions, and the process is as follows:

	Algorithm 1 Minimum travel time (i, j, $t_{dep}$)
	Input: Starting node i, end node j, departure time $t_{dep}$,         road data matrix R, Depth of water accumulation matrix D. Output: Minimum travel time $t{r}_{i,j}^k$, path length L
1.	Initialization:
	Define Ni as the number of nodes (derived from R) Define time[u] as the currently known shortest travel time from start node to u For all nodes u: time[u] = ∞ time[start node] = $t_{dep}$ Define p[u] as the predecessor of node u in the shortest path For all nodes u: p[u] = NaN Define pr[u] is a boolean indicating whether node u has been processed For all nodes u: pr[u] = false Define Q is a priority queue sorted by time Insert start node into Q
2.	Main Loop:
	While Q is not empty:    Extract the node current node with the smallest time from Q    If pr[current node] = true, continue to the next iteration    For each neighbor node of current node: Let Li,n = R[current node, neighbor node] If Li,n = ∞, skip this neighbor node Let dp = D[current node, neighbor node] Compute speed V = Speed function(dp) travel time on edge ti,n = Li,n/V If time[neighbor node] > time[current node] + ti,n:    time[neighbor node] = time[current node] + ti,n    p[neighbor node] = current node    Insert neighbor node into Q    pr[current node] = true    If current node = end node, break the loop
3.	Check Path:
	If time[end node] = ∞, output error (no feasible path)    Else:    path = []    node = end node    while node ≠ NaN: Insert node at the front of path node = p[node]
4.	Backtrack Path:
	L = 0    For each adjacent pair (i,j) in path:    L = L + R[i,j]
	L = 0    For each adjacent pair (i,j) in path:    L = L + R[i,j]
	Return: Minimum travel time $t{r}_{i,j}^k$, path length L

4. Two-Stage Solution Method

The model proposed in this paper is a mixed-integer second-order cone programming problem, which can be solved using the commercial solver Gurobi. The specific process is illustrated in Figure 3. Based on the rainfall return periods parameters, the DN failure probability model and the TN failure model are derived from the simulation of the rainfall process. Subsequently, the restoration models for the disaster progression period and post-disaster recovery are solved in periods to obtain the optimal load shedding and MESS scheduling schemes.

5. Case Study Analysis

This case study utilizes the IEEE 33-bus system and the Sioux-Falls24-node network. Geographic information is based on randomly generated elevations. The urban area is divided into a 30×30 grid, resulting in an improved 33-bus DN and a 24-node TN with geographical connection lines, as illustrated in Figure 4.

In the IEEE 33-bus system, bus 1 serves as the source bus, while the remaining buses are load buses. Specifically, buses 7, 10, 15, and 25 are points of distributed resource integration. The total system load is 8580 kW. The unit load shedding cost is 14 CNY/kWh.

An extreme heavy rain disaster process simulation lasting 240 min was conducted for the area where the DN was located. The evaluation time interval Td was set to 5 min. Based on the models presented in Section 2, the rainfall time series data are obtained, as shown in Figure 5. The rainfall meteorological parameters were based on historical data from a specific city in China [33]. The parameters were M = 9.898, k = 0.25, l = 0.656, and n = 2. As depicted in Figure 5, the maximum rainfall occurred at t = 55, with a rainfall intensity of approximately 9 mm/min and a cumulative rainfall of about 152 mm within one hour, categorizing it as an extremely heavy rain event. The total cumulative rainfall over 240 min reached 352 mm.

Urban drainage parameters are derived from the data in [34]. Numbers of drainage wells are 4602, S = 4602 × 0.05, n_m = 0.06, μ= 0.65, θ = 0.5, and k_f = 1. Based on the proposed model in this paper, the water depth in each grid was calculated and mapped to the corresponding bus in DN. The water depths at buses with different elevations are shown in Figure 6.

For the TN and DN failure models, the parameters are set as ε = 5000, γ = 2 × 10⁻⁴, H_bw = 0.3, H_w = 1, and d_p = 0.3. The repair time for crew in DN is 1 h, while the scheduling time scale for MESS in TN is 15 min. Disaster scenarios were generated based on bus failure probabilities. Buses 7, 8, and 33 failed, causing failures in lines 6, 8, and 32. Considering ground damage, additional failures occurred in lines 12, 18, 21, 24, and 25. The first-stage DN reconfiguration and restoration results are shown in Figure 7.

In the first stage, the DN was reconfigured using SOPs. Based on the distribution of tie lines, the original DN was divided into four supply islands, with a total load loss of 1654 kW. Once the rainfall dropped below the set threshold, the second stage of restoration began.

When the average depth of waterlogging across the entire area dropped below 0.3 m, the second-stage restoration was implemented. In this case, the second stage began at 180 min. In the second stage, DN was equipped with two MESS units, each with a 600 kWh capacity. Each MESS had a maximum discharge power of 200 kW. The MESS energy consumption rate was 0.371 kWh/km, with a distance cost of 50 CNY/km. Repair crew constraints allowed the simultaneous repair of faulty lines. The restoration process was divided into four periods, considering MESS and repair crew availability. Figure 8 shows the MESS dispatch results for the second stage, while Figure 9 displays DN reconfiguration and repair outcomes.

Figure 10 illustrates the voltage profiles at various buses during both restoration periods. The voltages remained within the quality standard range of 0.95 p.u.–1.05 p.u. throughout the restoration process.

During period 1, line repairs and reconfiguration adjusted the network structure, forming two islands. MESS1 moved from TN node 15 through node 19 to node 20, corresponding to movement from bus 8 to bus 13 in the DN. MESS2 moved from TN node 15 through node 10 to node 11, corresponding to movement from bus 8 to bus 33 in the DN. With MESS support, the load loss decreased to 1022 kW, while maintaining voltage levels within the 0.95 p.u.–1.05 p.u. requirement.

In period 2, MESS1 traveled through TN nodes 21 and 24 to reach node 13, corresponding to DN bus 22. MESS2 remained stationary at DN bus 33. At this point, line repairs had restored network connectivity, though some loads remained unrestored due to tie line usage.

In period 3, MESS1 traveled through TN nodes 12, 3 and 1 to reach node 2, corresponding to DN bus 25. MESS2 remained stationary at DN bus 33. At this point, line repairs had restored network connectivity, though some loads remained unrestored due to tie line usage.

In period 4, MESS1 returned to its station at node 15 via TN nodes 6, 9, and 10, and MESS2 returned to its station at node 15 via TN nodes 10, achieving complete load restoration.

Four restoration strategies were evaluated under identical scenarios:

Scheme 1: a single-stage restoration approach utilizing an SOP reconfiguration.

Scheme 2: the method proposed in [35] incorporates line repairs without MESS support for the DN.

Scheme 3: a simplified approach without urban network meshing, only reducing MESS operating speed during heavy rainstorm disasters.

Scheme 4: the method proposed in this study.

Table 1. Restoration results for each period for the four schemes.

Scheme	Load Loss For Each Period
	Stage 1	Stage 2
		1	2	3	4	5
Scheme 1	1654 kW	-	-	-	-	-
Scheme 2	1654 kW	1426 kW	854 kW	542 kW	180 kW	0
Scheme 3	1654 kW	1395 kW	701 kW	432 kW	0	0
Scheme 4	1654 kW	1395 kW	532 kW	288 kW	0	0

shows the restoration results for each time period for the four schemes. It can be seen that Scheme 1 did not fully recover the load after reconfiguration by SOPs alone due to the lack of consideration of overhaul and MESS support; in Scheme 2, although the lost load was eventually fully recovered, the time for restoration was long due to the lack of consideration of MESS support, and the amount of lost load was on the large side at all stages of restoration; in Scheme 3, when considering MESS scheduling, due to the irrational cuts in the speed of movement, the MESS failed to reach the maximum point of lost load required for support, resulting in a large amount of lost load at all stages of restoration. Scheme 4, by gridding the city and accurately considering the influence of waterlogging depth on MESS in different areas, enables the MESS to play a better supporting role in DN reconstruction. Compared to the first three schemes, Scheme 4 reduces the load shedding during the entire recovery process by 156.5%, 20.34%, and 8.09%, respectively.

In addition, to verify the generality of the proposed method, tests were conducted on the IEEE 123-bus system and five modified Sioux-Falls 24-node networks. The MESS placement and network topology reconfiguration results in the two stages are shown in Figure 11. Compared to the model without accurate TN and DN modeling, the IEEE 123-bus case achieved a 5.03% reduction in load loss.

6. Conclusions

This paper proposes an emergency power supply strategy for DNs under the impact of heavy rain disasters, enhancing the post-disaster restoration capability of DNs. By comprehensively considering the combined effects of precipitation intensity on both transportation and DNs, an indirectly coupled failure model is established between TNs and DNs. Building on this model, a two-stage restoration method is designed for DN reconfiguration and emergency resource scheduling. Experimental results demonstrate that this method significantly improves the efficiency and system stability of DNs during post-disaster restoration. The conclusions are as follows:

By integrating distribution areas and TNs through a grid-based approach and employing a two-dimensional surface hydrodynamic model, urban flooding scenarios are successfully constructed. This allows for a quantitative analysis of the combined impact of heavy rain disasters on DNs and transportation facilities. By simulating infrastructure damage mechanisms under varying precipitation intensities, the evolution processes of transportation and power systems are accurately quantified under heavy rain disasters.

Tailored to the distinct characteristics of different disaster stages, the proposed two-stage restoration strategy significantly improves load restoration efficiency under heavy rain disasters. By coordinating resources such as SOPs, line repairs, and MESS, the strategy optimizes the restoration process, ensuring greater load restoration and enhanced system resilience.

However, there are still some limitations in this study. First, the proposed model does not consider the uncertainties inherent in disaster scenarios, such as variations in precipitation intensity or the availability of resources, which may affect the robustness of the strategy. Second, the study focuses only on a single type of disaster, limiting its applicability to more complex multi-hazard environments.