Optimal Coordination of Transportable Power Sources and Repair Crews for Service Restoration of Distribution Networks Considering Uncertainty of Traffic Congestion

Zhao Shi; Yan Xu; Dunjian Xie; Shiwei Xie; Amer M. Y. M. Ghias

网刊加载中。。。

使用Chrome浏览器效果最佳，继续浏览，你可能不会看到最佳的展示效果，

确定继续浏览么?

复制成功，请在其他浏览器进行阅读

Optimal Coordination of Transportable Power Sources and Repair Crews for Service Restoration of Distribution Networks Considering Uncertainty of Traffic Congestion PDF

- ORCID：
Zhao Shi (Student Member, IEEE)
✉
- ORCID：
Yan Xu (Senior Member, IEEE)
✉
- ORCID：
Dunjian Xie (Student Member, IEEE)
✉
- ORCID：
Shiwei Xie (Member, IEEE)
✉
- ORCID：
Amer M. Y. M. Ghias (Senior Member, IEEE)
✉

Nanyang Technological University, Singapore； the School of Electrical Engineering and Automation, Fuzhou University, Fuzhou 350108, China

Updated：2024-01-22

DOI：10.35833/MPCE.2023.000012

OUTLINE

Abstract

This paper proposes a new method for service restoration of distribution network with the support of transportable power sources (TPSs) and repair crews (RCs). Firstly, a coupling model of distribution networks and vehicle routing of TPSs and RCs is proposed, where the TPSs serve as emergency power supply sources, and the RCs are used to repair the faulted lines. Considering the uncertainty of traffic congestion, the probability distribution of the travel time spent on each road is derived based on the Nesterov user equilibrium model, and a two-stage stochastic program is formulated to determine the optimal routings of TPSs and RCs. To efficiently solve the proposed stochastic mixed-integer linear program (MILP), a two-phase scenario reduction method is then developed to scale down the problem size, and an adaptive progressive hedging algorithm is used for an efficient solution. The effectiveness of the proposed methods and algorithms has been illustrated in a modified IEEE 33-bus system.

Keywords

Service restoration; distribution network; traffic uncertainty; transportable power source; repair crew; stochastic program

I. Introduction

IN recent decades, extreme weather events and natural catastrophes occur frequently [

1]. The consequence of extensive power outages and corresponding economic loss highlight the importance of enhancing the power system resilience [2]. After a contingency, service restoration is to dispatch available resources to restore the system to its normal state rapidly and efficiently. Hence, an effective service restoration is crucial to the resilience enhancement of the distribution networks (DNs).

The transportable resources have a great potential to improve the power system resilience due to their high mobility and flexibility. Transportable power sources (TPSs) and repaire crews (RCs) are two of the main transportable resources that have been widely applied in DNs. In particular, TPSs include large-capacity batteries or small generator sets carried by trucks or vehicles that could be timely dispatched to restore critical loads [

3], [4]. Meanwhile, RCs can be used for coupling repair and restoration processes. To fully restore a power grid after disasters, it is essential to dispatch RCs properly to repair the specific damaged components and support the network restoration of DNs [5], [6]. Reference [7] coordinates the scheduling of transportable energy storage (TES) with dynamic microgrid formation to minimize the total cost of DN after the disaster. Co-optimizing RC routing and reconfiguration for DN restoration is studied in [8]. However, most existing studies on coupling traffic routing models with DNs only consider either TPSs or RCs. Since both resources have strong interdependence with the power networks, it is essential to coordinate TPSs and RCs together within the coupling model of traffic routing and DNs.

The uncertainties of multiple components could also pose challenges in the service restoration process. Reference [

9] discusses the influence of wind energy uncertainty on power system restoration planning. Reference [10] designs a real-time recovery method for power networks with prepositioning and assignment of repair resources considering the uncertain propagation of the disaster. Reference [11] models load demand and repair time uncertainty in a DN repair and restoration problem. The unpredictable nature of weather is discussed in [12], [13] for its influence on power system resilience. For the problem with the routing of TPSs and RCs, the fluctuations of traffic demand in the traffic network result in recurrent congestion and the variability of travel time [14]. The corresponding uncertain link travel time will then affect the detailed dispatching plan of the disaster recovery scheme [15]. In reality, it is hard to predict the traffic demand between each origin-destination (OD) accurately, especially in some rural areas. Therefore, it is essential to consider the stochastic traffic demand when coupling the service restoration with the resource routing problem. A few studies, like [16], incorporate travel time uncertainty in an optimal post-disruption repairing schedule. Still, it is inapplicable for the service restoration problem since it does not consider the coupling model of the DN and traffic system.

Faced with the above issues, this paper proposes a two-stage joint stochastic service restoration method with the coordination of TPSs and RCs, considering the uncertainty of traffic congestion. A coupling model of DNs and vehicle routing of TPSs and RCs is proposed, where the TPSs serve as emergency power supply sources while the RCs determine the state of the faulted lines. To address the uncertainty of traffic congestion, this paper derives probability distributions of the travel time spent on each road based on the Nesterov user equilibrium (UE) model [

17] and then applies a two-stage stochastic program to determine the optimal routings of TPSs and RCs. Motivated by [18], the model size is reduced by the partitioning of the damaged components. When generating numerous scenarios in the stochastic program, a novel two-phase scenario reduction method is proposed to reduce the computational complexity. Moreover, the adaptive progressive hedging (A-PH) algorithm is utilized to decompose and solve the proposed stochastic mixed-integer linear program (MILP).

The main contributions of this paper can be summarized as follows.

1) A novel two-stage restoration framework is developed with the coordination of transportable resource dispatching and DN operation, where the first stage is to determine the dispatching decisions of TPS and RC while the second stage is to optimize DN operation. In this manner, all the resources from both traffic networks and power networks can be more effectively coordinated for fast and secure restoration.

2) Uncertain traffic congestion during the routing of TPS and RC is modeled based on the Nesterov UE model, which is effective to quantify the impacts of traffic conditions on the system restoration process under contingencies. To the best of the author’s knowledge, it is the first paper to address the uncertainty of traffic congestion in the system restoration problem.

3) A two-phase scenario reduction method and improved decomposition algorithm are proposed to efficiently solve the whole optimization model. Compared with the existing methods, both solution efficiency and accuracy are significantly improved.

The rest of this paper is organized as follows. Section II presents the problem description. Section III presents the mathematical formulation. Section IV details the solution method. Section V presents the case studies. Section VI concludes the paper.

II. Problem Description

A. Joint Stochastic Service Restoration

The entire framework of the proposed joint stochastic service restoration method is shown in Fig. 1, where the mathematical equations will be detailed in the sequel. The main objective of the proposed method is to maximize the amount of the restored loads by optimizing the decisions of routings of TPSs and RCs, topology reconfiguration, and DN power flow. Considering the uncertainty of traffic congestion, this paper designs a two-stage restoration framework based on stochastic optimization. The first stage is to determine the routing decisions of TPSs and RCs ahead of restoration implementation. The uncertain traffic congestion is considered as multiple potential scenarios of travel time with a corresponding probability. The objective function is modeled as the expected value of restored loads in all scenarios. The second stage is implemented after the restoration begins and uncertainty is realized. In the second stage, decisions of topology reconfiguration and DN power flow will be re-optimized to achieve optimal service restoration in real time.

Fig. 1 Framework of proposed joint stochastic service restoration method.

B. Uncertainty Handling

Multiple scenarios of travel time are generated from uncertain traffic congestion. To achieve it, this paper firstly derives a simplified traffic network (TN) by extracting the locations of depots, damaged components, and candidate charging points. Then, based on the Nesterov UE model, the probability distribution of travel time is derived from the traffic demand uncertainty between each OD pair. Finally, the travel time scenarios are generated based on the proposed two-phase scenario reduction method.

III. Mathematical Formulation

The detailed mathematical model of the proposed joint stochastic service restoration method could finally be formulated as a two-stage stochastic program. The objective and detailed constraints of these two stages are shown as follows.

A. Objective Function

Firstly, the objective function of the proposed two-stage stochastic program is given as:

\underset{λ}{m a x} E [\underset{y}{m a x} (\sum_{t} \sum_{i} W_{i} ε_{i}^{t} (ξ) \cdot P D_{i}^{t})]

(1)

The objective of the co-optimization model is to maximize the expected weighted sum of picked-up loads in the designed timeframe. Different load weights can be set based on priority levels [

19].

The item in the $(\cdot)$ is the optimal value of the second-stage problem and $E [\underset{y}{m a x} (\sum_{t} \sum_{i} W_{i} ε_{i}^{t} (ξ) \cdot P D_{i}^{t})]$ is the expected value of the second-stage problem. The inner max is used to gain the optimal value of the second stage with the first-stage decision, and the external max gains the total expected optimal maximum value.

B. Constraints for the First Stage

The routing problem aims to find an optimal routing for RCs and TPSs to travel among depots, damaged components, and candidate charging points. Assume that the routing problem can be defined by two undirected graphs $G_{R C} = (V, E_{R C})$ and $G_{T P S} = (N_{m}, E_{T P S})$ . The dp and tp in these two graphs represent the depot of each graph. Detailed routing constraints are modeled as:

\sum_{\forall n \in V \ \{m\}} α_{m, n, k}^{} = \sum_{\forall n \in V \ \{m\}} α_{n, m, k}^{} \forall k, m

(2)

\sum_{\forall m \in V \ \{d p\}} α_{d p, m, k}^{} - \sum_{\forall m \in V \ \{d p\}} α_{m, d p, k}^{} = 1 \forall k

(3)

\sum_{\forall m \in V \ \{t p\}} α_{m, t p, k}^{} - \sum_{\forall m \in V \ \{t p\}} α_{t p, m, k}^{} = 1 \forall k

(4)

Y_{m, k}^{T P S} = \sum_{\forall n \in N_{m} \ \{m\}} α_{m, n, k}^{} \forall k \in ω_{2}, m \in N_{m}

(5)

Y_{m, k}^{R C} = \sum_{\forall n \in V \ \{m\}} α_{m, n, k}^{} \forall k \in ω_{1}, m \in V

(6)

\sum_{\forall k \in ω_{1}} Y_{m, k}^{R C} \leq 1 \forall m \in V_{1}

(7)

Constraints from (2) to (7) represent the routing constraints for TPSs and RCs. Constraint (2) guarantees that each TPS or RC leaves the candidate charging point and damaged component once it completes the action, referred to as the flow conservation constraint in the vehicle routing problem (VRP) model [

20]. Constraints (3) and (4) define that each fleet starts from and returns to the designed depot after the assigned repair work is completed. Constraint (7) indicates that a maximum of one RC can fix each damaged line. The difference between TPSs and RCs is that the candidate charging points of TPSs can be visited multiple times by different TPS fleets.

C. Modeling for Uncertain Traffic Congestion

The uncertain traffic congestion can be reflected as multiple scenarios of travel time, which is expressed as a probability distribution. However, the traffic condition has a major impact on travel time. The improvement of relevant statistical methods and traffic observation tools allows for the accurate evaluation of road conditions and the monitoring of traffic demand [

21]. Nonetheless, to acquire the travel time of RCs and TPSs to any given destination, the topology of the traffic networks should also be considered. With the above consideration, the Nesterov UE model is used in this paper to determine the probability distribution of travel time based on traffic demand uncertainty.

The traffic demand uncertainty is modeled through a set of discrete scenarios, which are realizations of a uniform distribution with given upper and lower bounds $[(1 - ϱ) \bar{d_{i}}, (1 + ϱ) \bar{d_{i}}]$ [

22], where

ϱ

represents the uncertainty level of the traffic demand; and

\bar{d_{i}}

is the mean value of the demand. Since the traffic demand does not always vary on a short-time basis [23], [24], the road condition of each scenario could reach a state of equilibrium. Therefore, a scenario-based UE model [25] is employed to derive the link travel time of each scenario from the travel demand [26]:

m i n f_{T N} = C^{T N} \sum_{a} t r_{a} \cdot x_{a}

(8)

x_{a} = \sum_{π} \sum_{p} f_{p}^{π} δ_{a, p}^{π} \forall π \in T_{π}, p \in P^{π}

(9)

t r_{a} = t r_{a}^{0} [1 + ϕ {(\frac{x_{a}}{C_{a}})}^{4}] \forall a \in T_{A}

(10)

x_{a} \leq C_{a} \forall a \in T_{A}

(11)

\sum_{p} f_{p}^{π} = q_{π} (ξ) \forall π \in T_{π}, p \in P^{π}

(12)

\{\begin{array}{l} f_{p}^{π} (c_{p}^{π} - υ^{π}) = 0 \\ (c_{p}^{π} - υ^{π}) \geq 0 \end{array} \forall π \in T_{π}, p \in P^{π}

(13)

f_{p}^{π} \geq 0 \forall π \in T_{π}, p \in P^{π}

(14)

c_{p}^{π} = \sum_{a} t r_{a} \cdot δ_{a, p}^{π} \forall a, p

(15)

According to (8), the objective of this optimization process is to minimize the total path (link) travel cost in this equilibrium state. Constraint (9) states the relationship between the link and path flow. Based on the typical bureau of public roads (BPR) function, the link travel time of each scenario can be expressed as (10). Constraint (11) guarantees that the link flow is within the range of the link capacity. Constraint (12) denotes the traffic flow conservation, which represents that in each scenario, the sum of the path flow on each OD pair should meet the travel demand $q_{π}$ . Constraints (13)-(14) determine the equilibrium state. When the flow distribution reaches equilibrium, the equilibrium travel cost between each OD pair $υ^{π}$ is always no greater than that of any path. Constraint (15) explains that $c_{p}^{π}$ is the travel time (cost) of path p between OD $π$ . Since the travel time function (10) and UE condition (13) have several nonlinear items, e.g., the biquadratic item in (10), the piecewise linearization and a big M method are utilized to linearize nonlinear terms to ease the solving bottlenecks [

25].

By elaborating the above Nesterov UE model as the uncertainty revealing process, $t r_{a}$ in each scenario could be obtained. Related parameter settings of TN could be found in [

21], and the TN topology is shown in Appendix A Fig. A1.

D. Constraints for the Second Stage

1) Modeling of RCs

After determining the optimal routing of RCs at the first stage, the state of RCs should be determined considering the uncertainty of traffic congestion, which acts as the connection between the DN and the routing model. The travel time uncertainty in this paper, as different from other kinds of uncertainties like electricity load fluctuation, influences the RC arrival time at any given destination, and further affects the repair states of damages. However, the repair state is typically depicted in the general service restoration model based on the fixed time step, and it is difficult to integrate the travel time directly with the restoration model. Thus, we derive the following arrival time constraints to couple the RC arrival time and the corresponding repair states of damage. Use the binary variable $ζ_{m, t}^{R C}$ to indicate the state transition of the distribution line from the faulted state ( $ζ_{m, t}^{R C} = 0$ ) to the health state ( $ζ_{m, t}^{R C} = 1$ ). Relevant constraints are modeled as follows.

Firstly, arrival time constraints are stated as:

\begin{array}{l} a T_{m, k}^{R C} + r t_{m, k} + t r_{m, n} (ξ) - a T_{n, k}^{R C} \leq (1 - α_{m, n, k}^{}) M \\ \forall k \in ω_{1}, m \in V, n \in V_{1} \end{array}

(16)

\begin{array}{l} - (1 - α_{m, n, k}^{}) M \leq a T_{m, k}^{R C} + r t_{m, k} + t r_{m, n} (ξ) - a T_{n, k}^{R C} \\ \forall k \in ω_{1}, m \in V, n \in V_{1} \end{array}

(17)

For example, the RC $k$ arrives at the damaged component $m$ at time $a T_{m, k}^{R C}$ . Once arrived, $r t_{m, k}$ is spent for $k$ to repair $m$ . Then, after a $t r_{m, n} (ξ)$ travel period from $m$ to $n$ , the RC arrives at component $n$ at time $a T_{n, k}^{R C}$ . The big M method is applied to decouple the arrival time at $m$ and $n$ if RC $k$ does not travel from $m$ to $n$ . Then, to determine the repair completion time of each damaged component, the following formulae are enforced:

\sum_{\forall t} ζ_{m, t}^{R C} \leq 1 \forall m \in V_{1}

(18)

\sum_{\forall t} t ζ_{m, t}^{R C} \geq \sum_{\forall k} (a T_{m, k}^{R C} + r t_{m, k} \cdot Y_{m, k}^{R C}) \forall k \in ω_{1}, m \in V

(19)

\sum_{\forall t} t ζ_{m, t}^{R C} \leq \sum_{\forall k} (a T_{m, k}^{R C} + r t_{m, k} \cdot Y_{m, k}^{R C}) + 1 - ϑ \forall k \in ω_{1}, m \in V

(20)

0 \leq a T_{m, k}^{R C} \leq Y_{m, k}^{R C} M \forall k \in ω_{1}, m \in V

(21)

a T_{d p, k}^{R C} = 0 \forall k \in ω_{1}

(22)

Formula (18) enforces that each component is repaired once. Constraints (19) and (20) define the time when a component is repaired with the corresponding RC arrival time and the required repair time. For example, if RC $k = 1$ arrives at component $m = 1$ at time $t = 1$ , repair time $r t_{1,1} = 2$ is used to repair the component, so the damaged component is repaired at $t ζ_{1,3}^{R C} = 3$ . $ϑ$ is utilized since the schedule time horizon is the integer value. Constraint (21) indicates that $a T_{m, k}^{R C}$ equals 0 if m is not in the route of RC $k$ . Equation (22) enforces that all RCs depart from the depot when $t = 0$ .

Resource availability is guaranteed by constraint (23), which states that each RC’s resource capacity could satisfy the total resource demand of damaged components in its assigned route.

\sum_{m \in V} Y_{m, k}^{R C} \cdot r s_{m} \leq R S^{k} \forall k \in ω_{1}

(23)

At each time step, the repair results of the damaged components will reflect on the connection status of the DN and then affect the DN operation. The interdependent constraints between component repair and DN operation can be expressed as:

u_{i_{m} j_{m}}^{t + 1} \leq \sum_{τ = 1}^{t} ζ_{m, τ}^{R C} \forall m \in V_{1}, \forall t \leq T s - 1

(24)

u_{i j}^{t} \leq u_{i j}^{t + 1} \forall (i, j) \in E, \forall t \leq T s - 1

(25)

u_{i j}^{t} = 1 \forall (i, j) \in E \ L_{1}, \forall t

(26)

u_{i j}^{1} = 0 \forall (i, j) \in L_{1}

(27)

a_{i j}^{t} \leq u_{i j}^{t} \forall (i, j) \in E, \forall t

(28)

Constraint (24) enforces that the damaged lines will be operable once an RC repairs it in the previous time step. Constraints (25)-(27) restrict the line status. Constraint (25) indicates that lines should remain operable once it is repaired. Constraint (26) states that all unbroken lines are operable. Constraint (27) sets the initial status of the damaged line to complement the time horizon of $u_{i j}^{t}$ . Constraint (28) is the independent constraint that states one line could be closed only if it is operable.

2) Modeling of TPSs

The TPSs, including TES and truck-mounted emergency generators (TEGs), should be appropriately dispatched among the candidate charging points to supply critical loads after the outage. Similar to RC, the arrival time constraints are also employed to derive the TPS travel time to the operation time step. The scheduling principle of the TPSs is roughly similar to that of RCs with differences: one damaged component could be repaired only once by one RC, but one candidate charging point could be connected by more than one TPS at the same time. In addition, the repair time for RCs to repair the damaged component is a parameter, but the time one TPS spends at a candidate station is a variable. Detailed constraints are expressed as:

\begin{array}{l} a T_{m, k}^{T P S} + c t_{m, k} + t r_{m, n} (ξ) - a T_{n, k}^{T P S} \leq (1 - α_{m, n, k}^{}) M \\ \forall k \in ω_{2}, m \in N_{m}, n \in N_{m}^{} \ \{m\} \end{array}

(29)

\begin{array}{l} - (1 - α_{m, n, k}^{}) M \leq a T_{m, k}^{T P S} + c t_{m, k} + t r_{m, n} (ξ) - a T_{n, k}^{T P S} \\ \forall k \in ω_{2}, m \in N_{m}, n \in N_{m}^{} \ \{m\} \end{array}

(30)

The arrival time constraints (29) and (30) of TPSs are similar to those of RCs. With travel time variation $t r_{m, n} (ξ)$ , the time one TPS arrives at and leaves a charging point could be influenced. The major difference is that the time one TPS spends at a candidate charging point $c t_{m, k}$ is a variable.

The accessing state of the TPSs depends on the arrival time and the length of stay, which is modeled in detail as:

\sum_{\forall t} ζ_{m, t, k}^{T P S} \leq 1 \forall k \in ω_{2}, m \in N_{m}^{}

(31)

\sum_{\forall t} ζ_{m, t, k}^{T P S} = \sum_{\forall t} l_{m, t, k}^{T P S} \forall k \in ω_{2}, m \in N_{m}^{}

(32)

a T_{m, k}^{T P S} \leq \sum_{\forall t} t ζ_{m, t, k}^{T P S} \leq a T_{m, k}^{T P S} + 1 - ϑ \forall k \in ω_{2}, m \in N_{m}

(33)

\sum_{\forall t} t l_{m, t, k}^{T P S} \geq \sum_{\forall k} (a T_{m, k}^{T P S} + c t_{m, k}) \forall k \in ω_{2}, m \in N_{m}

(34)

\sum_{\forall t} t l_{m, t, k}^{T P S} \leq \sum_{\forall k} (a T_{m, k}^{T P S} + c t_{m, k}) + 1 - ϑ \forall k \in ω_{2}, m \in N_{m}

(35)

0 \leq a T_{m, k}^{T P S} \leq Y_{m, k}^{T P S} M \forall k \in ω_{2}, m \in N_{m}

(36)

0 \leq c t_{m, k} \leq T s \forall k \in ω_{2}, m \in N_{m}

(37)

a T_{d p, k}^{T P S} = 0 \forall k \in ω_{2}

(38)

Constraint (31) states that each TPS only visits one candidate charging point once. Flow conservation constraint is guaranteed by constraint (32), which enforces that once a TPS arrives at a candidate charging point, it must also leave it. Similar to constraints (19) and (20), constraints (33)-(35) determine the time that one TPS arrives at and leaves a charging point with its corresponding arrival time and length of stay. For example, if TPS $k = 1$ arrives at charging point $m = 1$ at time $t = 2$ and discharges for $c t_{1,1} = 1$ hour, $t ζ_{1,2, 1}^{T P S} = 2$ and $t l_{1,3, 1}^{T P S} = 2 + 1 = 3$ . Constraint (36) indicates that the arrival time $a T_{m, k}^{T P S}$ equals 0 if TPS $k$ does not visit candidate charging point m. Constraint (37) restricts the range of the time one TPS spends at a candidate station. Constraint (38) indicates the initial location of all TPSs at $t = 0$ .

Besides the time-related constraints, constraint (39) restricts the number of TPSs that can be connected to one candidate charging point.

\sum_{\forall k} z_{m, t}^{k} \leq c a p_{i} \forall k \in ω_{2}, m \in N_{m}^{}, \forall t

(39)

To better couple the dispatch of TPSs with the DN operation, a binary variable $z_{m, t}^{k}$ is introduced to represent the status of the TPS.

z_{m, t}^{k} = \sum_{τ = 1}^{t - 1} (ζ_{m, t, k}^{T P S} - l_{m, t, k}^{T P S}) \forall k \in ω_{2}, m \in N_{m}^{}

(40)

\sum_{\forall m} z_{m, t}^{k} \leq 1 \forall m \in N_{m}^{}, \forall t, k \in ω_{2}

(41)

Specifically, constraint (40) indicates that $z_{m, t}^{k} = 1$ if and only if the TPS is staying at the candidate charging point; otherwise, it equals 0 if the TPS did not arrive at the candidate charging point or already left the point. Constraint (41) enforces that one TPS could only visit one candidate charging point.

The dispatch and operation of TPSs, e.g., discharging of the TES and TEG, could also affect the DN operation:

P_{i}^{t} = \sum_{k \in M_{1}} z_{i, t}^{k} P_{k, t}^{g e} + \sum_{k \in M_{2}} z_{i, t}^{k} (P_{k, t}^{d} - P_{k, t}^{c}) \forall i \in N_{m}, \forall t

(42)

Q_{i}^{t} = \sum_{k \in M} z_{i, t}^{k} Q_{k, t}^{g e} \forall i \in N_{m}, \forall t

(43)

P_{i}^{t} = Q_{i}^{t} = 0 \forall i \in N \ (N_{m} ⋃ i_{g}), \forall t

(44)

Constraints (42) and (43) indicate the contribution of the TPS output to the DN. The total injection from the TPS is the sum of the output of TESs and TEGs. Constraint (44) restricts the power injection to zero for the buses without connection to TPSs or substations.

It should be noted that there are several bilinear items in (42) and (43). A typical reformulation-linearization method [

27] converts the non-linear model into an MILP.

During the scheduling process of the TPSs, when they are staying at the candidate charging points, power charging or discharging between TPSs and the DN will occur. The operation constraints of the mentioned two types of the TPSs can be modeled as:

0 \leq P_{k, t}^{g e} \leq \sum_{i \in N_{m}} z_{i, t}^{k} {\bar{P}}_{k}^{g e} \forall k \in M_{1}, \forall t

(45)

0 \leq Q_{k, t}^{g e} \leq \sum_{i \in N_{m}} z_{i, t}^{k} {\bar{Q}}_{k}^{g e} \forall k \in M_{1}, \forall t

(46)

0 \leq P_{k, t}^{c} \leq I c_{t}^{k} \cdot {\bar{P}}_{k}^{c} \forall k \in M_{2}, \forall t

(47)

0 \leq P_{k, t}^{d} \leq I d_{t}^{k} \cdot {\bar{P}}_{k}^{d} \forall k \in M_{2}, \forall t

(48)

I c_{t}^{k} + I d_{t}^{k} \leq \sum_{i \in N_{m}} z_{i, t}^{k} \forall k \in M_{2}, \forall t

(49)

0 \leq Q_{k, t}^{g e} \leq (I c_{t}^{k} + I d_{t}^{k}) {\bar{Q}}_{k, t}^{g e} \forall k \in M_{2}, \forall t

(50)

s o c_{t + 1}^{k} = s o c_{t}^{k} + (P_{k, t}^{c} η_{k}^{c} - P_{k, t}^{d} / η_{k}^{d}) Δ t \forall k \in M_{2}, \forall t \leq T s - 1

(51)

{\underset{̲}{s o c}}^{k} \leq s o c_{t}^{k} \leq {\bar{s o c}}^{k} \forall k \in M_{2}, \forall t

(52)

Constraints (45) and (46) are the real and reactive power output limits of TEGs. These two constraints enforce that only if the TEG is connected to a candidate charging point, both $P_{k, t}^{g e}$ and $Q_{k, t}^{g e}$ could have positive values. Charging and discharging power limits of TESs are indicated in (47) and (48). Constraint (49) states the operation modes of the TES. Specifically, only if $\sum_{i \in N_{m}} z_{i, t}^{k} = 1$ , which means that the TES is staying at a candidate charging point, it could operate in charging or discharging mode; otherwise, both $I c_{t}^{k}$ and $I d_{t}^{k}$ will be zero. Similar to (46), (50) restricts the reactive power output of the TES. The state of charge (SOC) variation of the TES is modeled by (51); and (52) limits the upper and lower bounds of SOC for the TES.

3) Operation Constraints of DNs

A linearized Distflow model is utilized in this paper for DN operation analysis [

28], [29]. The detailed formulation can be expressed as:

P G_{i}^{t} - ε_{i}^{t} \cdot P D_{i}^{t} = \sum_{(i, j) \in E} P_{i j}^{t} - \sum_{(k, i) \in E} P_{k i}^{t} \forall i \in N, t \in T s

(53)

Q G_{i}^{t} - ε_{i}^{t} \cdot Q D_{i}^{t} = \sum_{(i, j) \in E} Q_{i j}^{t} - \sum_{(k, i) \in E} Q_{k i}^{t} \forall i \in N, t \in T s

(54)

V_{j}^{t} - V_{i}^{t} \leq M (1 - a_{i j}^{t}) - (P_{i j}^{t} r_{i j} + Q_{i j}^{t} x_{i j}) / V_{0} \forall (i, j), t

(55)

V_{j}^{t} - V_{i}^{t} \geq - M (1 - a_{i j}^{t}) - (P_{i j}^{t} r_{i j} + Q_{i j}^{t} x_{i j}) / V_{0} \forall (i, j), t

(56)

- a_{i j}^{t} S_{i j}^{m a x} \leq P_{i j}^{t} \leq a_{i j}^{t} S_{i j}^{m a x} \forall (i, j) \in E, t \in T s

(57)

- a_{i j}^{t} S_{i j}^{m a x} \leq Q_{i j}^{t} \leq a_{i j}^{t} S_{i j}^{m a x} \forall (i, j) \in E, t \in T s

(58)

- \sqrt[]{2} a_{i j}^{t} S_{i j}^{m a x} \leq P_{i j}^{t} + Q_{i j}^{t} \leq \sqrt[]{2} a_{i j}^{t} S_{i j}^{m a x} \forall (i, j) \in E, t \in T s

(59)

- \sqrt[]{2} a_{i j}^{t} S_{i j}^{m a x} \leq P_{i j}^{t} - Q_{i j}^{t} \leq \sqrt[]{2} a_{i j}^{t} S_{i j}^{m a x} \forall (i, j) \in E, t \in T s

(60)

V_{i}^{m i n} \leq V_{i}^{t} \leq V_{i}^{m a x} \forall i \in N \ D, t \in T s

(61)

ε_{i}^{t} \leq ε_{i}^{t + 1} \forall i, \forall t \leq T s - 1

(62)

Constraints (53) and (54) represent the real and reactive power balance constraints. The line voltage drop constraints are formulated by (55) and (56). The big M method decouples the disconnected buses [

30]. Constraints (57) to (60) indicate the line flow limitations and enforce that the line flow equals zero if the line is open. Constraint (61) sets the upper and lower bounds of the voltage magnitude. Constraint (62) enforces that a load should keep energized during the following period once it is restored.

The DN is dynamically reconfigured in the restoration process to maintain the radial structure by using remotely controlled switches. In this paper, the spanning forest model is used to ensure the radial arrangement of the DN during each step of the restoration process. To strictly guarantee the radial construction, the single-commodity flow (SCF) model is utilized in this paper as the radiality constraints [

31]:

\sum_{(i, j) \in E} b_{i j}^{t} = |N| - |N s| \forall t

(63)

\sum_{(j, i) \in E} c f_{j i}^{t} - \sum_{(i, j) \in E} c f_{i j}^{t} = 1 \forall i \in N \ i_{g}, \forall t

(64)

- b_{i j}^{t} \cdot \bar{c f} \leq c f_{i j}^{t} \leq b_{i j}^{t} \cdot \bar{c f} \forall (i, j) \in E, \forall t

(65)

\bar{c f} = |N| - |N s|

(66)

a_{i j}^{t} \leq b_{i j}^{t} \forall (i, j) \in E, \forall t

(67)

a_{i j}^{t} = 1 \forall (i, j) \in E \ \{L_{1} ⋃ L_{2}\}, \forall t

(68)

As the SCF model states, a DN with N buses and Ns substations could remain a radial topology only if the following two conditions are satisfied.

1) The DN should have $N - N s$ closed lines.

2) All buses of the DN should be connected.

The first condition is satisfied by (63). The single-commodity flow constraints (64) and (65) are used to guarantee the second condition. Constraint (64) states that the substation should feed each bus with one unit of fictitious load. Constraint (65) ensures that there is no fictitious flow on disconnected lines. The capacity value of the fictitious flow is set as $|N| - |N s|$ , as shown in (66).

To ensure the radial structure of the DN in each time step, the binary variable $b_{i j}^{t}$ is set as the connection status of the line $(i, j)$ at time t in fictitious network, and the variable $a_{i j}^{t}$ is the connection status of the line $(i, j)$ at time t in real DN. Therefore, constraint (67) guarantees that at each time step, the closed lines in DN are a subset of lines in a fictitious network so that the feasible solutions $a_{i j}^{t}$ could form a spanning forest at each time t. Except for the damaged lines and lines with remotely-controlled switches, the connection status of all other lines are kept closed by (68).

Finally, the extensive form (EF) of the proposed two-stage stochastic program is formulated by objective (1) and constraints (2)-(7) and (16)-(68) presented above. For better understanding, the constraints are classified in detail as:

1) Constraints for the first stage: routing of RCs and TPSs: (2)-(7).

2) Constraints for the second stage: ① modeling of RCs: (16)-(28); ② modeling of TRSs: (29)-(52); and ③ operation constraints of DNs: (53)-(68).

Combined with the clarification of the general model in (1), the first-stage variables of vector $λ$ contain binary variables $α_{m, n, k}^{}$ , $Y_{m, k}^{R C}$ , $Y_{m, k}^{T P S}$ which determine the routing of TPSs and RCs. Then, the vector $y$ denotes the remaining decision valuables in the second stage, i.e., $ζ_{m, t}^{R C}$ , $ζ_{m, t, k}^{T P S}$ , $l_{m, t, k}^{T P S}$ , $z_{m, t}^{k}$ , $ε_{i}^{t}$ , $a_{i j}^{t}$ , $b_{i j}^{t}$ , $u_{i j}^{t}$ , $I c_{t}^{k}$ , $I d_{t}^{k}$ , $a T_{m, k}^{R C}$ , $a T_{m, k}^{T P S}$ , $c t_{m, k}$ , $P_{k, t}^{g e}$ , $Q_{k, t}^{g e}$ , $P_{k, t}^{c}$ , $P_{k, t}^{d}$ , $s o c_{t}^{k}$ , $P G_{i}^{t}$ , $Q G_{i}^{t}$ , $P_{i j}^{t}$ , $Q_{i j}^{t}$ , $V_{i}^{t}$ , $c f_{i j}^{t}$ .

IV. Solution Method

A. Scenario-based Two-stage Stochastic Optimization

To handle the challenge posed by the uncertainty source, scenario-based two-stage stochastic optimization is utilized. As described in Section III, the stochastic variation of the traffic demand between each OD pair is modeled through a known distribution. The Monte Carlo sampling method is implemented to generate multiple demand realizations. Then, possible scenarios of different travel demands are reduced to a reasonable number through a novel two-phase scenario reduction method, which will be detailedly introduced in the following subsection.

Through the employing of the Nesterov UE model, the link travel time realizations (in terms of scenarios) which reflect the uncertain traffic congestion are derived. After the first stage of decision-making, the uncertain travel time is revealed and is considered in the second stage of decision-making. As shown in (1), two-stage stochastic program aims to find a solution that maximizes the expected restored load over all the simulated scenarios. Thus, all the second-stage variables shall be scenario-related and could be expressed with the added subscript of s (denoted as scenario), i.e., $ζ_{m, t, s}^{R C}$ , $ζ_{m, t, k, s}^{T P S}$ , $l_{m, t, k, s}^{T P S}$ , $z_{m, t, s}^{k}$ , $ε_{i, s}^{t}$ , $a_{i j, s}^{t}$ , $b_{i j, s}^{t}$ , $u_{i j, s}^{t}$ , $I c_{t, s}^{k}$ , $I d_{t, s}^{k}$ , $a T_{m, k, s}^{R C}$ , $a T_{m, k, s}^{T P S}$ , $c t_{m, k, s}$ , $P_{k, t, s}^{g e}$ , $Q_{k, t, s}^{g e}$ , $P_{k, t, s}^{c}$ , $P_{k, t, s}^{d}$ , $s o c_{t, s}^{k}$ , $P G_{i, s}^{t}$ , $Q G_{i, s}^{t}$ , $P_{i j, s}^{t}$ , $Q_{i j, s}^{t}$ , $V_{i, s}^{t}$ , $c f_{i j, s}^{t}$ .

B. Depot-based Partitioning of Damaged Components

After a large disaster, local neighboring utilities can provide assistance by deploying additional RCs to help restore the destroyed network. When there are multiple RC depots in the disaster-affected area, the computational burden can be reduced by first partitioning the damages to different depots [

8]. By assigning damaged components to depots, the entire destroyed network can be divided into several small clusters/regions. Thus, the large-scale NP-hard VRP for the routing problem can be decomposed into several small VRPs. Therefore, the model scale and computational complexity are greatly reduced.

The detailed partitioning model is formulated as:

m i n \sum_{m \in V_{1}} \sum_{n \in D} β_{m, n} d_{m, n}

(69)

\sum_{n \in D} β_{m, n} \leq 1 \forall m \in V_{1}

(70)

u_{i_{m} j_{m}}^{} \leq \sum_{n \in D} β_{m, n} \forall m \in V_{1}

(71)

\sum_{(i, j) \in E} P_{i j}^{} - \sum_{(k, i) \in E} P_{k i}^{} - P D_{i}^{} = 0 \forall i \in N

(72)

\sum_{(i, j) \in E} Q_{i j}^{} - \sum_{(k, i) \in E} Q_{k i}^{} - Q D_{i}^{} = 0 \forall i \in N

(73)

(26), (28), (55)-(61), (63)-(68) (without superscript t)

(74)

The objective of the above small MILP (69) is to minimize the total distance between the damages and their assigned depots. Constraint (70) restricts that one component can be clustered to a maximum of one depot. Motivated by [

18], constraints (71)-(74) guarantee that rather than all the damaged components, the above model is applied to select a minimum number of damages that could fully pick up the DN loads without the support of TPSs. The schematic of the depot-based partitioning method is given in Fig. 2.

Fig. 2 Schematic of depot-based partitioning method.

C. Two-phase Scenario Reduction

1) Phase 1: Clustering for Point Operation

A large number of demand realizations of each OD could cause a relatively huge number of combinations of different OD demands as the total traffic demand matrix scenarios. Since the operation objects in this phase are the points that represent the possible demand values of a specific OD, a clustering technique is applied to select several representative points to form the demand matrix for a further reduction in the second phase.

In this phase, a modified K-means algorithm, K-means++, is implemented for clustering data. The K-means algorithm has the characteristics of low computational complexity and fast speed [

32], which is suitable for the aforementioned point value operation. Furthermore, the K-means++ algorithm determines the cluster centers sequentially according to their distances, which will make the distribution of selected cluster centers more uniform.

2) Phase 2: Kantorovich Distance (KD) for Matrix Scenarios

After clustering in phase 1, a smaller representative number of demand values for each OD is obtained. Then, after the combination, a computable number of traffic demand matrix realizations are selected for the secondary reduction in this phase. Here, an efficient KD-based backward method is used to reduce the scenarios [

33]. By constructing the Kantorovich distance matrix (KDM), scenarios could be eliminated based on their KD with other scenarios and their probability of occurrence.

Finally, by applying the proposed two-phase scenario reduction method shown in Fig. 3, a small representative number of scenarios remain. Thus, the model size is significantly reduced, and the solving process can be sped up greatly.

Fig. 3 Framework of two-phase scenario reduction method.

D. A-PH Algorithm

The two-stage stochastic problem can be solved as a single-stage large linear program with duplicated constraints in each scenario when the random vector has a finite number of scenarios. However, when the number of integer variables is large in each scenario, the computational burden could be increased by directly solving the single-stage large linear program. This section proposes an A-PH algorithm, as shown in Algorithm 1, for decomposing and solving the stochastic MILP, where $𝒬_{s}$ represents the subproblem constraints in each scenario; $P_{s}^{r}$ is the probability of scenario s; $τ_{1}$ and $τ_{2}$ denote at which iteration $ρ$ is adjusted to be larger or smaller; $β_{1}$ and $β_{2}$ denote in which proportion $ρ$ is adjusted to be larger or smaller; and $σ^{ν}$ denotes the solution gap of iteration v.

Firstly, in Steps 1-3, a non-anticipative “initial guess” is obtained by solving the scenario subproblems in the initialization phase. Then, by applying a penalty parameter $ρ$ , the non-anticipative is enforced through the updating of the multiplier $φ_{s}^{ν}$ in Step 5, where the superscript v is the iteration number and the subscript s denotes the scenario counter. The subproblems are then solved by augmenting linear and quadratic proximal items in Step 6. Finally, as shown in Steps 8 and 11, a non-anticipative solution is yielded once all first-stage decisions converge on a common $\bar{λ}$ .

Algorithm 1 : the proposed A-PH algorithm
Step 1: initialization: $ν \leftarrow 0$ , $φ_{s}^{ν} \leftarrow 0$ , $ρ \leftarrow ρ_{0}$ , $i \leftarrow 1$ , $j \leftarrow 1$ , $\forall s \in S$
Step 2: iteration 0: $\forall s \in S$ , $λ_{s}^{ν} = \underset{λ}{a r g m i n} \{m^{T} λ + n_{s}^{T} y_{s} : (λ, y_{s}) \in Q_{s}\}$
Step 3: aggregation: ${\bar{λ}}^{ν} = \sum_{s \in S} P_{s}^{r} λ_{s}^{ν}$
Step 4: iteration update: $ν \leftarrow ν + 1$
Step 5: multiplier update: $φ_{s}^{ν} \leftarrow φ_{s}^{ν - 1} + ρ (λ_{s}^{ν - 1} - {\bar{λ}}^{ν - 1}),$ $\forall s \in S$
Step 6: iteration $ν$ : $\forall s \in S,$ $λ_{s}^{ν} = \underset{λ}{a r g m i n} \{m^{T} λ + n_{s}^{T} y_{s} + φ_{s}^{ν} λ + \frac{ρ}{2} {‖λ - {\bar{λ}}^{ν - 1}‖}^{2} : (λ, y_{s}) \in Q_{s}\}$
Step 7: aggregation: ${\bar{λ}}^{ν} = \sum_{s \in S} P_{s}^{r} λ_{s}^{ν}$
Step 8: $σ^{ν} = \sum_{s \in S} P_{s}^{r} ‖λ - {\bar{λ}}^{ν}‖$
Step 9: threshold value check: if $σ^{ν - 1} - σ^{ν} \leq ψ_{1} σ^{ν - 1}$ , then $i \leftarrow i + 1$ . Else if $σ^{ν - 1} - σ^{ν} \geq ψ_{2} σ^{ν - 1}$ , then $j \leftarrow j + 1$ . Else, $i \leftarrow 1$ , $j \leftarrow 1$ . End if.
Step 10: penalty factor adjustment: if $i = τ_{1}$ , then $ρ = (1 + β_{1}) ρ$ ; $i \leftarrow 1$ . Else if $j = τ_{2}$ , then $ρ = (1 + β_{2}) ρ$ ; $j \leftarrow 1$ . End if.
Step 11: convergence check: if $σ^{ν} < ϵ$ , halt. Otherwise, go to Step 4.

From Algorithm 1, we know that $ρ$ is the parameter that determines the movement range of $λ$ and in turn, affects the PH convergence. If $ρ$ is small, PH convergence will be greatly delayed due to minor movement in $λ$ . If $ρ$ is large, the accuracy of the optimization results cannot be guaranteed. Therefore, the optimal $ρ$ needs not to be a fixed value for each stochastic MILP. However, most existing research applies the fixed $ρ$ value to obtain the PH convergence. Some research works use a fixed position to adjust the $ρ$ value, but that position needs to be changed once the problem parameter changes.

In the proposed A-PH algorithm, the penalty factor is self-adjusting during the solving process according to the change of the gap valve. The change of gap value $σ^{v}$ is checked in Step 9 during each iteration. If the changes remain in a small range after a continuous fixed number of iterations, the penalty factor will increase proportionally as Step 10 shows. In contrast, if the changes remain in a relatively large range after a continuous fixed number of iterations, $ρ$ will decrease proportionally. These actions will repeat the whole solving process until the program reaches convergence. By doing this, the adjustment position of the $ρ$ value could be adaptive and will do self-adjusting according to different problems.

V. Case Studies

This section uses a modified IEEE 33-bus system as the test case for the proposed restoration procedure [

34]. The routing of TPSs and RCs in two comparison cases is shown in Fig. 4. The proposed method is modeled by using MATLAB 2021a, and all subproblems and the depot-based partitioning model are solved by Gurobi (or other off-the-shelf solvers), on an Intel^(R) Xeon^(R) 3.7 GHz PC with 16 GB RAM.

Fig. 4 Routing of TPSs and RCs in two comparison cases. (a) Benchmark case. (b) Stochastic case.

A. Test System

As Fig. 4 shows, we assume six damaged lines, two RCs, and two TPS fleets (one TES and one TEG) in this test system. Both RCs and TPSs are supposed to have identical properties. Two RCs start from one depot, and the initial point of both TPS fleets is set at the same candidate charging point. Assume that the outages happen at midnight, and let the time step equals 0.5 hour. Both required repair time and resources for different damaged components are listed in Table I. For resource constraints, the resource capacity of the RC is set to be 10. Detailed parameters of TPSs are shown in Table II.

TABLE I Required Repair Time and Resources for Different Damaged Components

Damaged component	Required repair time (time step)		Required resources (units) of all RCs
Damaged component	RC1	RC2	Required resources (units) of all RCs
D1	1	2	3
D2	2	1	3
D3	1	2	1
D4	1	1	3
D5	2	1	2
D6	1	2	1

TABLE II Detailed Parameters of TPSs

Type	Parameter	Value
TEG	Initial position (DN bus number)	8.0
	Real power (MW)	0.5
	Reactive power (Mvar)	0.3
TES	Initial position (DN bus number)	8.0
	Charging/discharging power (MW)	0.2
	Energy capacity (MWh)	0.2
	Initial $\underset{̲}{s o c}$	0.9
	$s o c_{m a x}, s o c_{m i n}$	0.9, 0.1
	Charging/discharging efficiency	0.95

For the DN test feeder, the upper and lower bounds of the voltage value are set to be $\pm 5 %$ of the nominal level, which is 1.0 p.u.. According to constraint (58), the commodity flow capacity is set to be 32. In the optimal case, each line could be equipped with a remotely controlled switch. To reduce the number of switching actions and related costs [

35], seven switches are allocated (lines 9-10, 23-24, 28-29, 8-21, 12-22, 18-33, 25-29). The buses without switches can be expressed as “branch blocks”, which could be used to simplify the topology structure of the test system, and then reduce the computational complexity [36].

Meanwhile, the number of damaged components that need to be repaired could also influence the calculation complication of the proposed method. In the aforementioned test case with one RC depot, the assigning algorithm in Section IV is implemented to pre-assign the minimum set of repair tasks to the depots. Four lines are selected by processing this assigning algorithm before the main optimization, and the problem size is further reduced.

B. Simulation Results and Analysis

1) Comparison Study

The co-optimization problem under the deterministic situation is set as the benchmark for comparison. The conventional deterministic case is tested without considering traffic uncertainty, and the traffic demand is set as their expected values.

To show the importance of the uncertainty consideration, the deterministic dispatch solution of RCs and TPSs is tested using 30 random scenarios of traffic demands, as shown in Fig. 5, where all the colored lines represent random scenarios. The load of around 30% of scenarios cannot be thoroughly picked up under the deterministic benchmark case during the scheduled period.

Fig. 5 Load restoration of deterministic benchmark case in random scenarios.

Therefore, based on the main aim of load restoration, it is necessary to consider the traffic demand uncertainty in the proposed service restoration method. One thousand random scenarios for each OD are generated using the Monte Carlo sampling technique to describe the stochastic nature of the traffic demand. The uncertainty level of the travel demand is set to be 0.4, and the value of each OD demand is given in Table III. Then, the proposed two-phase scenario reduction method is applied. After the first phase of clustering, $5^{4} = 625$ groups of different OD demand combinations remain. Then, after the second-phase reduction with KD, ten traffic demand matrixes (scenarios) are selected as the most representative ones for stochastic program.

TABLE III Values of OD Demand for 10 Representative Scenarios

OD pairs	Expected value	Value of OD demand
OD pairs	Expected value	1	2	3	4	5	6	7	8	9	10
π₄(T₁-T₁₂)	70	41.680	41.680	70.310	97.920	55.660	55.660	55.660	84.150	84.150	84.150
π₃(T₁-T₁₁)	60	60.360	72.260	72.260	72.260	49.280	60.360	72.260	49.280	49.280	60.360
π₁(T₁-T₆)	50	69.150	38.470	38.470	38.470	38.470	69.150	38.470	38.470	38.470	69.150
π₂(T₁-T₁₀)	70	53.340	96.110	80.050	53.340	66.000	80.050	53.340	53.340	96.110	53.340
Probability		0.072	0.064	0.096	0.094	0.123	0.093	0.107	0.114	0.110	0.126

It can be observed from Fig. 4 that in the stochastic case, one more line is assigned to RC1 for repair. Since traffic demand varies and traffic congestion may occur in this case, increasing repair work could ensure the restoration plan with sufficient capability to fully restore the system load. Besides, the repair order of RC2 changes. Compared with the benchmark case, though microgrids are connected to the main grid, loads like 27 still need extra connections like line 32-33 to satisfy operational constraints for picking up. Since the deterministic travel time from depot to line 32-33 is also farther than that of line 19-20, repairing line 32-33 firstly makes the restoration plan more robust when facing traffic demand fluctuation. In addition to the RCs, the routing of TES1 also changes. Since traffic congestion may occur on its travel paths, to be more robust against the possible congestion, TES1 is assigned to reduce the number of position changes and stay at suitable points for supporting power. Meanwhile, during the restoration process, lines 25-29 and 18-33 are permanently assigned to be closed; and bus 29 is a good injection point to support the maximum loads based on the topology. Therefore, TES1 changes its routes and is assigned to be connected to bus 29 directly.

Regarding the scheduling solution of the second stage, we use one possible realization (scenario) for detailed analysis in the stochastic case. The dispatching solution of RCs in two comparison cases is shown in Table IV. We take TES1 as the example for analyzing the scheduling solution. Buses 8, 15, 29, and 25 of the DN are represented by charging points 1-4. The fluctuation of the traffic demand causes congestion on the paths of TPSs. Thus, in this scenario, TPSs take more time to arrive at their assigned locations. Therefore, from Fig. 6, we can observe that TES1 starts injecting power from $t = 2$ in the benchmark case, but in the stochastic case, the discharging of TES1 starts from $t = 3$ due to the traffic congestion. During the following period, in the benchmark case, the amount of TES1 power injection remains high until line 32-33 is repaired. TES1 is needed to satisfy the operational constraints of buses at the end of the topology. However, in this test case, since the lines 32-33 and 9-15 are repaired at $t = 5$ and $t = 6$ firstly, TES1 could be disconnected after $t = 5$ . Finally, as soon as line 19-20 is repaired, all TPS fleets could be disconnected, and network loads could be fully restored.

TABLE IV Dispatching Solution of RCs in Two Comparison Cases

Case	Fleet	Time step
Case	Fleet	0	1	2	3	4	5	6	7	8
Benchmark	RC1	Depot	Travel	Line 2-3	Line 2-3	Return to depot	Return to depot	Return to depot	Return to depot	Return to depot
Benchmark	RC2	Depot	Travel	Line 19-20	Line 19-20	Travel	Line 32-33	Return to depot	Return to depot	Return to depot
Stochastic	RC1	Depot	Travel	Line 2-3	Line 2-3	Travel	Lines 9-15	Return to depot	Return to depot	Return to depot
Stochastic	RC2	Depot	Travel	Travel	Line 32-33	Line 32-33	Travel	Line 19-20	Return to depot	Return to depot

Fig. 6 Scheduling of TES1 in two comparison cases. (a) Benchmark case.(b) Selected scenario in stochastic case.

Unlike the deterministic case, Fig. 7 shows that by solving the two-stage stochastic model, loads in all scenarios could be fully restored in the designed planning period. Table V shows the statistical results of two comparison cases. Through modeling with the traffic uncertainty, the mean value of the objective has been improved and the variance of objective values in stochastic case is smaller than that in the benchmark case. Therefore, it means that the proposed joint stochastic service restoration method could generate a more robust recovery plan and could better deal with the uncertain traffic condition for DN load restoration.

Fig. 7 Load restoration for 10 representative scenarios by stochastic case.

TABLE V Statistical Results of Two Comparison Cases

Case	Optimal value (MWh)	Mean value (MWh)	Variance (MWh)
Benchmark	5770.9	5556.7	518540.6
Stochastic	5560.9	5578.2	494590.3

2) Performance of Solution

Firstly, the proposed joint stochastic service restoration method on the aforementioned test case is utilized to demonstrate the performance of the proposed methods. The performance comparison is shown in Table VI.

TABLE VI Performance Comparison Among Different Algorithms

Algorithm	Computation time (min)	Objective value (MWh)	Whether reaches optimality
EF	120	5390.1	No ( $% g a p = 3.2 % > 0.01 %$ )
PH	120	5472.6	No ( $σ = 2.3 > 0.01$ )
A-PH	15.3	5560.9	Yes ( $σ = 1.6 \times 10^{- 14} < 0.01$ )

In Table VI, the %gap is the computational gap obtained from directly solving the EF using Gurobi, and the last column is the relative gap $σ$ we set in PH and A-PH algorithms. We can observe that the A-PH algorithm converges in 15.3 min, the solution of which is the highest among the three algorithms, and the relative gap is below the given threshold $ϵ = 0.01$ . The %gap shows that solving EF cannot gain an optimal solution in 2 hours. The high complexity of EF also results in the lowest objective value in the designed timeframe. According to the $σ$ value, the PH also does not achieve an optimal solution with a constant $ρ$ . In test cases, the same initial penalty factor value is applied in both PH and A-PH algorithms. If the penalty factor is fixed during the solving process, PH may not be able to converge to a smaller gap when compared with A-PH algorithm in reasonable timeframes.

Moreover, the advantage of using PH over EF is also shown from the results. Though both solutions are not optimal, PH gains a higher objective value than EF.

Therefore, the computational burden could be significantly reduced by creating the decomposed stochastic program using the A-PH algorithm while a reasonable accuracy is maintained. Thus, more representative scenarios could be used, and solutions with more robustness will be obtained.

To make the proposed joint stochastic service restoration method more realistic and applicable against the unexpected contingency, the performance analysis regarding the damage partitioning with multiple depots is conducted, as shown in Table VII.

TABLE VII Performance Analysis Regarding Damage Partitioning with Multiple Depots

Case	Objective value (MWh)	Computation time (s)
Case 1	5560.9	918.0
Case 2	7062.3	402.0
Case 3	7062.3	117.5

Case 1 is described in Section V-A with only one depot which contains two RCs. Case 2 is the one with two depots (each with one RC) and without partitioning of damages. Case 3 has two depots (each with one RC), and the damages are divided into two regions through the method in Section IV. As shown in Table VII, by comparing the performance between Case 1 and Case 2, the total served loads increase, and the computation time is reduced when the number of depots increases from 1 to 2.

When the total served loads are basically the same, it can be observed from the comparison between Case 2 and Case 3 that by dividing the damaged network into smaller clusters/regions, the computation time of Case 3 is greatly reduced to 117.5 s.

Thus, it can be concluded that through the partitioning of damages with multiple depots, the computational complexity can be greatly reduced, and a fast response can be realized after the contingency.

VI. Conclusion

This paper proposes a joint stochastic service restoration method to co-optimize the DN repair and restoration. The dispatching of TPSs and RCs and dynamic DN restoration strategies are coordinated together with traffic uncertainty. The mathematical model of the proposed method is formulated as a two-stage stochastic program. In the first stage, the optimal routings of TPSs and RCs are decided. The second stage derives the scheduling sequences of the discharging of TPSs and the corresponding DN reconfiguration. The proposed method is finally linearized and transformed into an MILP. An improved decomposition algorithm, A-PH algorithm, is applied to solve the proposed stochastic MILP. The numerical results demonstrate the effectiveness and robustness of the joint stochastic service restoration method compared with the deterministic case. The results also show that the improved decomposition method could significantly reduce the computational complexity while maintaining the necessary accuracy. Moreover, to ensure a fast response regarding post-disaster recovery, the advantages of employing the depot-based partitioning technique have also been verified.

Nomenclature

Symbol	——	Definition
A.	——	Indices and Sets
$ω_{1}, ω_{2}$	——	Sets for repair crews (RCs) and transportable power sources (TPSs)
$E$	——	Set of distribution network (DN) lines, indexed by line $(i, j)$
$i_{g}$	——	Set of DN substations
$i_{m}, j_{m}$	——	Indices of damaged line m
$k$	——	Index of TPSs and RCs
$L_{1}, L_{2}$	——	Sets of damaged lines and lines with switches
$M_{1}$	——	Set of truck-mounted emergency generators (TEGs)
$M_{2}$	——	Set of transportable energy storage (TESs)
m, n	——	Indices of traffic nodes for damaged components and charging points, and depots
$N$	——	Set of DN buses, indexed by line (i, j)
$N_{s}$	——	Set of DN substations
$N_{m}$	——	Set of TPS charging points
$P^{π}$	——	Set of paths belonging to each origin-destination (OD) pair $π$
$T_{A}$	——	Link set of road
$T_{π}$	——	Set of OD pairs
$V$	——	Set of damaged components and depots
$V_{1}$	——	Set of damaged components
B.	——	Parameters
$ϑ$	——	Error factor
$C_{a}$	——	Capacity of road link a
$c a p_{i}$	——	Allowed number of TPSs connected to point i
$d_{m, n}$	——	Distance between component m and depot n
${\bar{P}}_{k}^{g e}, {\bar{Q}}_{k}^{g e}$	——	The maximum active and reactive power outputs of TEG k
${\bar{P}}_{k}^{c}, {\bar{P}}_{k}^{d}$	——	The maximum charging and discharging power of TES k
$P D_{i}^{t}, Q D_{i}^{t}$	——	Active and reactive loads at bus i at time t
$q_{π}$	——	Traffic demand of OD pair $π$ in each scenario
$R S^{k}$	——	Resource capacity of RC $k$
$r_{i j}, x_{i j}$	——	Resistance and reactance of line (i, j)
$r s_{m}$	——	Required resources to repair component m
$r t_{m, k}$	——	Time for RC $k$ to repair component m
$t r_{m, n}$	——	Travel time from component m to depot n
$T s$	——	Number of time steps
$W_{i}$	——	Weight of load at bus i
C.	——	Variables
$α_{m, n, k}^{}$	——	Binary variable indicating whether RC or TPS $k$ travels from component m to depot n
$β_{m, n}$	——	Binary variable indicating whether component m is assigned to depot n
$λ$ , $y$	——	Decision vectors in the first and second stages
$ξ$	——	Random vector
$δ_{a, p}^{π}$	——	Binary variable, which equals 1 when path p goes through link a
$ζ_{m, t}^{R C}$	——	Binary variable, which equals 1 if component $m$ is repaired at time t
$ζ_{m, t, k}^{T P S}$	——	Binary variable equals 1 if charging point m is visited by TPS k at time t
$ε_{i}^{t}$	——	Binary variable indicating operation status of load at bus i, and $ε_{i}^{t} = 1$ means the load is restored at time t in scenario s, $ε_{i}^{t} = 0$ , otherwise
$a_{i j}^{t}$	——	Connection status of line (i, j) at time t in real DN
$a T_{m, k}^{R C}$	——	Time when RC $k$ arrives at component m
$a T_{}^{T P S}$	——	Time when TPS k arrives at charging point
$b_{i j}^{t}$	——	Connection status of line (i, j) at time t in fictitious network
$C_{}^{T N}$	——	Unit monetary value of traffic time
$c f_{j i}^{t}$	——	Fictitious flow among line (i, j) at time t
$c_{p}^{π}$	——	Travel time (cost) on path p between OD pair $π$ in each scenario
$c t_{m, k}$	——	Time TPS k spends at charging point m
$f_{p}^{π}$	——	Traffic flow on path p between OD pair $π$ in each scenario
$I c_{t}^{k}, I d_{t}^{k}$	——	Charging and discharging statuses of TPS k at time t
$l_{m, t, k}^{T P S}$	——	Binary variable, which equals 1 if TPS k leaves charging point m at time t
$P_{k, t}^{g e}, Q_{k, t}^{g e}$	——	Active and reactive power outputs of TEG k at time t
$P_{k, t}^{c}, P_{k, t}^{d}$	——	Charging and discharging power of TES k at time t
$P G_{i}^{t}, Q G_{i}^{t}$	——	Active and reactive power generated at bus i at time t
$P_{i j}^{t}, Q_{i j}^{t}$	——	Active and reactive power flows on line (i, j) at time t
$t_{a}^{0}$	——	Free travel time on link a
$t r_{a}$	——	Travel time on link a in each scenario
$u_{i j}^{t}$	——	Operation status of DN line (i, j) at time t
$V_{i}^{t}$	——	Voltage magnitude of bus i at time t
$x_{a}$	——	Traffic flow on link a in each scenario
$Y_{m, k}^{R C}$	——	Binary variable, which equals 1 if component m is fixed by RC $k$
$Y_{m, k}^{T P S}$	——	Binary variable, which equals 1 if TPS k visits charging point m
$z_{m, t}^{k}$	——	Binary variable, which equals 1 if TPS k is staying at charging point m at time t

Appendix

Appendix A

Fig. A1 Traffic network topology with coupled components.

References

C. Chen and B. Chen, “Modernizing distribution system restoration to achieve resiliency against extreme weather events,” in Proceedings of 2018 IEEE Global Conference on Signal and Information Processing, Anaheim, USA, Nov. 2018, pp. 895-896. [Baidu Scholar]

H. Gao, Y. Chen, S. Mei et al., “Resilience-oriented pre-hurricane resource allocation in distribution systems considering electric buses,” Proceedings of the IEEE, vol. 105, no. 7, pp. 1214-1233, Jul. 2017. [Baidu Scholar]

S. Yao, P. Wang, and T. Zhao, “Transportable energy storage for more resilient distribution systems with multiple microgrids,” IEEE Transactions on Smart Grid, vol. 10, no. 3, pp. 3331-3341, May 2019. [Baidu Scholar]

S. Yao, P. Wang, X. Liu et al., “Rolling optimization of mobile energy storage fleets for resilient service restoration,” IEEE Transactions on Smart Grid, vol. 11, no. 2, pp. 1030-1043, Mar. 2020. [Baidu Scholar]

P. Van Hentenryck, C. Coffrin, and R. Bent, “Vehicle routing for the last mile of power system restoration,” in Proceedings of 17th Power Systems Computation Conference, Stockholm, Sweden, Nov. 2011, pp. 836-843. [Baidu Scholar]

G. Zhang, F. Zhang, X. Zhang et al., “Sequential disaster recovery model for distribution systems with co-optimization of maintenance and restoration crew dispatch,” IEEE Transactions on Smart Grid, vol. 11, no. 6, pp. 4700-4713, Nov. 2020. [Baidu Scholar]

S. Lei, C. Chen, H. Zhou et al., “Routing and scheduling of mobile power sources for distribution system resilience enhancement,” IEEE Transactions on Smart Grid, vol. 10, no. 5, pp. 5650-5662, Sept. 2019. [Baidu Scholar]

A. Arif, Z. Wang, J. Wang et al., “Power distribution system outage management with co-optimization of repairs, reconfiguration, and DG dispatch,” IEEE Transactions on Smart Grid, vol. 9, no. 5, pp. 4109-4118, Sept. 2018. [Baidu Scholar]

A. Golshani, W. Sun, Q. Zhou et al., “Incorporating wind energy in power system restoration planning,” IEEE Transactions on Smart Grid, vol. 10, no. 1, pp. 16-28, Jan. 2019. [Baidu Scholar]

A. Inanlouganji, G. Pedrielli, T. A. Reddy et al., “A computational approach for real-time stochastic recovery of electric power networks during a disaster,” Transportation Research Part E: Logistics and Transportation Review, vol. 163, p. 102752, Jul. 2022. [Baidu Scholar]

A. Arif, S. Ma, Z. Wang et al., “Optimizing service restoration in distribution systems with uncertain repair time and demand,” IEEE Transactions on Power Systems, vol. 33, no. 6, pp. 6828-6838, Nov. 2018. [Baidu Scholar]

M. Panteli and P. Mancarella, “Influence of extreme weather and climate change on the resilience of power systems: impacts and possible mitigation strategies,” Electric Power Systems Research, vol. 127, pp. 259-270, Oct. 2015. [Baidu Scholar]

H. Chen, J. Wang, J. Zhu et al., “A two-stage stochastic mixed-integer programming model for resilience enhancement of active distribution networks,” Journal of Modern Power Systems and Clean Energy,vol. 11, no. 1, pp. 94-106, Jan. 2023. [Baidu Scholar]

S. Clark and D. Watling, “Modelling network travel time reliability under stochastic demand,” Transportation Research Part B: Methodological, vol. 39, no. 2, pp. 119-140, Feb. 2005. [Baidu Scholar]

H. Gan and Y. Bai, “The effect of travel time variability on route choice decision: a generalized linear mixed model based analysis,” Transportation, vol. 41, no. 2, pp. 339-350, Mar. 2014. [Baidu Scholar]

B. A. Alkhaleel, H. Liao, and K. M. Sullivan, “Risk and resilience-based optimal post-disruption restoration for critical infrastructures under uncertainty,” European Journal of Operational Research, vol. 296, no. 1, pp. 174-202, Jan. 2022. [Baidu Scholar]

J. G. Wardrop, “Road paper. Some theoretical aspects of road traffic research,” Proceedings of the Institution of Civil Engineers, vol. 1, no. 3, pp. 325-362, May 1952. [Baidu Scholar]

S. Lei, C. Chen, Y. Li et al., “Resilient disaster recovery logistics of distribution systems: co-optimize service restoration with repair crew and mobile power source dispatch,” IEEE Transactions on Smart Grid, vol. 10, no. 6, pp. 6187-6202, Nov. 2019. [Baidu Scholar]

W. Li, Risk Assessment of Power Systems: Models, Methods, and Applications. New York: IEEE Press, Jan. 2005. [Baidu Scholar]

D. C. Paraskevopoulos, G. Laporte, P. P. Repoussis et al., “Resource constrained routing and scheduling: review and research prospects,” European Journal of Operational Research, vol. 263, no. 3, pp. 737-754, Dec. 2017. [Baidu Scholar]

W. Wei, L. Wu, J. Wang et al., “Network equilibrium of coupled transportation and power distribution systems,” IEEE Transactions on Smart Grid, vol. 9, no. 6, pp. 6764-6779, Nov. 2018. [Baidu Scholar]

Q. Tian, Y. Lin, and D. Wang, “Autonomous and conventional bus fleet optimization for fixed-route operations considering demand uncertainty,” Transportation, vol. 48, no. 5, pp. 2735-2763, Oct. 2021. [Baidu Scholar]

H. Shao, W. H. K. Lam, and M. L. Tam, “A reliability-based stochastic traffic assignment model for network with multiple user classes under uncertainty in demand,” Networks and Spatial Economics, vol. 6, pp. 173-204, Aug. 2006. [Baidu Scholar]

S. Nakayama and D. Watling, “Consistent formulation of network equilibrium with stochastic flows,” Transportation Research Part B, vol. 66, pp. 50-69, Aug. 2014. [Baidu Scholar]

S. Xie, Z. Hu, and J. Wang, “Scenario-based comprehensive expansion planning model for a coupled transportation and active distribution system,” Applied Energy, vol. 255, p. 113782, Dec. 2019. [Baidu Scholar]

J. Wang, M. Ehrgott, and A. Chen, “A bi-objective user equilibrium model of travel time reliability in a road network,” Transportation Research Part B: Methodological, vol. 66, pp. 4-15, Aug. 2014. [Baidu Scholar]

H. D. Sherali and W. P. Adams, “A hierarchy of relaxations and convex hull characterizations for mixed-integer zero-one programming problems,” Discrete Applied Mathematics, vol. 52, no. 1, pp. 83-106, Jul. 1994. [Baidu Scholar]

T. Ding, Y. Lin, G. Li et al., “A new model for resilient distribution systems by microgrids formation,” IEEE Transactions on Power Systems, vol. 32, no. 5, pp. 4145-4147, Sept. 2017. [Baidu Scholar]

Z. Li, Y. Xu, X. Feng et al., “Optimal stochastic deployment of heterogeneous energy storage in a residential multienergy microgrid with demand-side management,” IEEE Transactions on Industrial Informatics, vol. 17, no. 2, pp. 991-1004, Feb. 2021. [Baidu Scholar]

T. Ding, Y. Lin, Z. Bie et al., “A resilient microgrid formation strategy for load restoration considering master-slave distributed generators and topology reconfiguration,” Applied Energy, vol. 199, pp. 205-216, Aug. 2017. [Baidu Scholar]

M. Lavorato, J. F. Franco, M. J. Rider et al., “Imposing radiality constraints in distribution system optimization problems,” IEEE Transactions on Power Systems, vol. 28, no. 1, pp. 568-569, Feb. 2012. [Baidu Scholar]

N. Shi, X. Liu, and Y. Guan, “Research on k-means clustering algorithm: an improved k-means clustering algorithm,” in Proceedings of 3rd International Symposium on Intelligent Information Technology and Security Informatics, Jian, China, Apr. 2010, pp. 63-67. [Baidu Scholar]

L. Wu, M. Shahidehpour, and T. Li, “Stochastic security-constrained unit commitment,” IEEE Transactions on Power Systems, vol. 22, no. 2, pp. 800-811, May 2007. [Baidu Scholar]

M. E. Baran and F. F. Wu, “Network reconfiguration in distribution systems for loss reduction and load balancing,” IEEE Power Engineering Review, vol. 9, no. 4, pp. 101-102, Apr. 1989. [Baidu Scholar]

S. Lei, J. Wang, and Y. Hou, “Remote-controlled switch allocation enabling prompt restoration of distribution systems,” IEEE Transactions on Power Systems, vol. 33, no. 3, pp. 3129-3142, May 2018. [Baidu Scholar]

B. Chen, C. Chen, J. Wang et al., “Sequential service restoration for unbalanced distribution systems and microgrids,” IEEE Transactions on Power Systems, vol. 33, no. 2, pp. 1507-1520, Mar. 2018. [Baidu Scholar]

Zhao Shi received the B.E. degree from Shandong University, Ji’nan, China, in 2017, and the M.S. degree from the University of Bath, Bath, UK, in 2018. He is currently pursuing the Ph.D. degree with the School of Electrical and Electronic Engineering, Nanyang Technological University, Singapore. His research interests include power system optimization, mobile power source integration, and interdependent network resilience. [Baidu Scholar]

Yan Xu received the B.E. and M.E degrees from South China University of Technology, Guangzhou, China, in 2008 and 2011, respectively, and the Ph.D. degree from The University of Newcastle, Newcastle, Australia, in 2013. He conducted postdoctoral research with the University of Sydney Postdoctoral Fellowship, Sydney, Australia, and then joined Nanyang Technological University (NTU), Singapore, with The Nanyang Assistant Professorship. He is now an Associate Professor at School of Electrical and Electronic Engineering and a Cluster Director at Energy Research Institute @ NTU (ERI@N), Singapore. He is serving as the Chairman for IEEE Power & Energy Society Singapore Chapter. His research interests include power system stability and control, microgrid, and data analytics for smart grid applications. [Baidu Scholar]

Dunjian Xie received both B.Eng and M.Sc degrees in electrical engineering from Zhejiang University, Hangzhou, China, in 2016 and 2019, respectively. He is currently persuing a Ph.D. degree at Nanyang Technological University, Singapore. His research interests include modeling and optimization of the demand side resources and their further application to enhance resilience of the power grid. [Baidu Scholar]

Shiwei Xie received the Ph.D. degree in electrical engineering from Wuhan University, Wuhan, China, in 2021. From 2019 to 2020, he was a Research Assistant with the School of Electrical and Electronic Engineering (EEE), Nanyang Technological University, Singapore. He is currently a Lecturer (Associated Researcher) with the School of electrical engineering and automation, Fuzhou University, Fuzhou, China. His current research interests include variational inequality theory, distributed optimization, robust optimization, and their applications in power and transportation systems. [Baidu Scholar]

Amer M. Y. M. Ghias received the B.Sc. degree in electrical engineering from Saint Cloud State University, St Cloud, USA, in 2001, the M.Eng. degree in telecommunications from the University of Limerick, Limerick, Ireland, in 2006, and the Ph.D. degree in electrical engineering from the University of New South Wales (UNSW), Sydney, Australia, in 2014. From February 2002 to July 2009, he held various positions such as Electrical Engineer, Project Engineer, and Project Manager, while working with the top companies in Kuwait. He was with UNSW, during 2014-2015, and the University of Sharjah, United Arab Emirates, during 2015-2018. In 2018, he joined the Nanyang Technological University, sigapore, as an Assistant Professor. He is also a Cluster Director (Power Electronics and Energy Management) for Energy Research Institute @ NTU (ERI@N), Singapore. His research interests include model predictive control, hybrid energy storage, renewable energy sources, multiphase drives, new multilevel converters, and advanced modulations for multilevel converter. [Baidu Scholar]

Address:No.19 Chengxin Avenue, Jiangning District, Nanjing 211106, China

E-mail: mpce@alljournals.cn

Tel:86-25-81093060

Fax:86-25-81093040

Home

Introduction

Editorial Board

For Author

Call For Papers

APC

Sponsor & Publisher