A Rolling Horizon Restoration Framework for Post-disaster Restoration of Electrical Distribution Networks

Abstract

Severe weather events such as floods, hurricanes, earthquakes, and large wind or ice storms can cause extensive damage to electrical distribution networks, requiring a multi-day restoration effort. Complicating the recovery process is the lack of complete and accurate information regarding the extent and locations of damages, at least during the initial part of the recovery process. These factors make workforce planning challenging. In this paper, we adopt a rolling horizon restoration framework whereby repairs are planned for adjustable finite length restoration windows. Considering both repair times as well as travel times, we show that the optimal scheduling problem with multiple crews, each with their own time budget, can be recast in terms of a cost constrained reward maximizing $m$TSP (traveling salesman problem) on doubly weighted graphs, where the objective is to maximize the aggregate reward earned during the upcoming restoration window, provided no crew violates its time budget and certain electrical continuity constraints are met. We propose a mixed integer linear programming (MILP) model for solving the above problem which is validated on standard IEEE PES test feeder networks.

I Introduction

Disasters, such as typhoons, earthquakes, and other natural calamities, can cause severe damage to distributed power networks. Worldwide, severe weather events have been reported to be the leading cause of major electrical outages [1]. For instance, an ice storm in 2023 severely damaged power lines in Texas, resulting in a blackout that lasted several days and affected more than 400,000 customers [2, 3]. It has been reported that the U.S. economy loses an estimated $20 to $55 billion each year due to weather-related outages [4, 5, 6]. Post-disaster recovery planning can be extremely challenging due to factors such as: (i) the large number of equipment that must be repaired or replaced, (ii) accessibility to job sites may be limited, and (iii) insufficient repair workforce and/or equipment [7]. Consequently, effective management of both equipment and personnel resources necessitates careful coordination and strategic planning. This involves prioritizing essential repairs and optimizing the travel schedule to reduce the unused work time.

A key impediment toward optimizing the restoration process is incomplete or uncertain information regarding: (i) repair locations, (ii) road infrastructure conditions, and (iii) availability of repair personnel. Since postponing the recovery process until better information is available is a luxury that planners can ill afford, the process usually begins with incomplete or uncertain information, and repair schedules/horizons are adjusted along the way as better and more accurate information starts to become available.

In this paper, we adopt a rolling horizon restoration framework where the recovery timeline is split into adjustable finite horizon windows. Repair activities for a particular restoration window are scheduled based on the information available until the end of the prior restoration window. This information could be some combination of newly acquired damage information, refinement of prior damage information, knowledge of completed repairs, and availability of repair equipment/personnel. Updating the damage and resource information before scheduling the repairs for the next restoration window helps adapt the recovery process to prevailing conditions on the ground.

We frame the optimal repair of the distribution network as a version of the traveling salesman problem (TSP) on a doubly weighted graph, i.e., a graph with node weights (modeling the repair times) and edge weights (modeling the travel times between repair sites). We assume that a certain number of repair crews ($m$) is available for the next restoration window, and that each crew has its own time budget, which need not be equal to the length of the window. To each node (repair job), we associate a reward which can be earned if the job is completed during the next window and certain electrical continuity constraints (explained later) are met. However, instead of a traditional TSP (or $m$TSP), we solve a model which closely resembles an instance of a cost constrained reward maximizing $m$TSP on doubly weighted graphs (CCRM-$m$TSP-DW), where the objective is to maximize the aggregate reward earned during the upcoming restoration window, provided no crew violates its time budget and certain additional electrical continuity constraints are met. Additionally, we show that if the power network is of star topology, the electrical continuity constraints are redundant and the restoration problem is exactly equivalent to an instance of the CCRM-$m$TSP-DW (see Section V-C). The reward maximization objective can be changed to a profit maximization objective by assigning a penalty parameter to each node, to be incurred if a job is not selected for repair during the upcoming window. Qualitatively speaking, a properly designed time-adaptive penalty parameter can introduce an element of ‘fairness’ to the solutions, so that repair jobs offering low reward (i.e., jobs viewed as low priority) are not automatically relegated to the back end of the recovery process. To the best of our knowledge, these variants of the TSP have not been studied before and may be of independent interest to researchers in other areas such as computer science and operations research.

II Literature Review

In addressing the complex problem of repairing electrical networks after disasters, researchers have explored a variety of approaches. Studies [8, 9, 10] built mathematical models under two simplifying assumptions: ignoring the travel time between repair nodes and assuming only a single crew was available. While these assumptions helped simplify the problem and speed up the optimization process, they did not fully meet the practical demands of real-world scenarios. Tan et al. in [10] touches upon multi-crews, which was thoroughly considered in [11, 12, 13].

As the field evolved, Ma et al. [14] and Arif et al. [15] both specifically pointed out the equal importance of accounting for both travel and repair times. To tackle the added complexity this brought, a number of studies [16, 17, 18, 19, 20] adopted co-optimization strategies, solving the problem by breaking it down into multiple stages or optimizing multiple objectives concurrently. Additionally, Markov decision process (MDP) [21], machine learning [22], and reinforcement learning algorithms [23] based on the idea of decision trees [24] were investigated.

After extreme weather disasters, the damage to the electrical network updates progressively over time. Accounting for this dynamic damage information is another critical aspect of the repair scheduling problem. In [25, 26, 17, 21, 27], the authors introduced rolling window techniques to adapt to continuous updates in network damage. This innovation reduces the complexity of the model, shortens the solving times, decreases the computational costs, and lets the scheduling start immediately.

Scholarly work has focused on the practical aspects of the issue, including network power load [17, 19, 20], cost-efficiency [25], and the average outage duration for consumers [21]. MILP models have been shown to be effective in optimizing repair routes and schedules. Arif et al. [15] discussed how to simultaneously allocate repair personnel and equipment resources, and Lei et al. [18] investigated how to dispatch the repair crews and mobile power sources at the same time.

However, prior research has seldom considered the different remaining work time budgets for each repair crew, although Xu et al. [13] conducted a preliminary study on incorporating this constraint. Therefore, this research seeks to enhance the scheduling and allocation of repair tasks by integrating a broader spectrum of practical factors, aiming to maximize the completion of repairs within a constrained rolling time frame. We concurrently consider the travel time, multiple crews, and the disparities in the remaining work time budgets among crews within a consistent rolling window. The developed MILP approach ensures a more efficient and effective allocation of resources, contributing significantly to the optimization of repair task scheduling in dynamic environments.

III Problem Definition

We begin with an explanation of the rolling horizon restoration framework we consider in this paper. Referring to Figure 1, let $t=0$ denote the time corresponding to the onset of a disaster event. Repair and restoration work commences at time $T_{0}$, after adequate damage information has been collected and initial repair crews mobilized. The time between $T_{k}$ and $T_{k-1}$ denotes the $k^{th}$ restoration window (planning horizon). We assume that damage information collection is a streaming process and commencement of repair and restoration efforts is not contingent on the availability of complete damage information. The set of all possible repair jobs during the $k^{th}$ window, $L_{k}^{D}$, includes previously known but unscheduled or unfinished jobs from prior windows, as well as new jobs learnt during the $(k-1)^{th}$ window. For notational simplicity, we will drop the ‘window subscript’ and simply use $L^{D}$ to denote the complete set of repair jobs for any restoration window. Associated with each damaged line (in this paper, we will use the terms ‘edge’, ‘line’, and ‘job’ interchangeably) is a reward, $r_{l}$, which we assume is known, determined by the utility companies. While it is certainly possible to assign identical rewards to all lines, we will discuss some alternate approaches later in the paper.

Figure 1: Timeline of the restoration process, which commences at time $T_{0}$ after a disaster at $t=0$. The set of all possible line/edge repairs (jobs) during the $k^{th}$ restoration window, $[T_{k-1},T_{k}]$, comprises previously known but unscheduled jobs from prior windows and new jobs learnt during the $(k-1)^{th}$ restoration window (denoted by $J_{k-1}$). Unfinished jobs from previous windows, possibly due to time overruns, can also be included in the set.

Let $m$ denote the number of repair crews available to work during the $k^{th}$ window. The planning objective for the $k^{th}$ window is to determine the repair schedules for the $m$ crews, choosing jobs from the set $L^{D}$, such that the total reward that can be earned at the end of the current window is maximized. The cost of any schedule is the sum of the repair times of the jobs on the schedule and the travel times necessary for moving from one job site to the next one on the schedule. We assume that the travel times between any pair of jobs in $L^{D}$ are symmetric and correspond to the shortest path (derived from a separate transportation map) between the two job sites. Further, we assume that any job site is reachable from any other job site, failing which a job site which is unreachable, possibly due to loss of transportation infrastructure, should be removed from the set $L^{D}$. Obviously, the cost (length) of any schedule must not exceed the duration of the $k^{th}$ window, $T_{k}-T_{k-1}$. Any line, $l$, repaired during the $k^{th}$ window earns the reward $r_{l}$ only if it can be energized from the previously energized portion of the distribution network, at the end of the $k^{th}$ restoration window. This feasibility depends on whether a fully repaired path exists from the energized portion of the distribution network at $T_{k}$ to $l$ (electrical continuity constraint).

We do not assume that the repair crews are deployed to their first job sites from a common staging location. In practice, this would require provisioning of adequate slack times between successive restoration windows in the restoration timeline. The travel times to the first job sites on the $m$ schedules can then be factored into these slack times and the $k^{th}$ restoration window can be deemed to start when all crews are at their initial job sites. Further, since the repair and travel times are estimated quantities and susceptible to delays, the slack times can be used to provide some cushion against unanticipated time overruns. In this paper, however, we ignore the presence of these slack times.

Mathematically, the broad objective is to construct $m$ repair schedules in order to:

$$\text{maximize}\ \sum_{l\in L^{D}}r_{l}\,I\left(\mathcal{P}(l,S_{k-1})\right),$$

where $S_{k-1}$ denotes the energized portion of the distribution network at time $T_{k-1}$ (the ‘source’) and $I\left(\mathcal{P}(l,S_{k-1})\right)$ is an indicator function which is equal to $1$ if a completely repaired path exists from $l$ to $S_{k-1}$ at time $T_{k}$ and is equal to $0$ otherwise. Henceforth, we will simply use ‘path’ to mean a ‘completely repaired path’ to the source. If the $r_{l}$’s are equal for all $l$, the objective function reduces to maximizing the number of repairs which can earn an immediate reward at the end of the current restoration window.

We conclude this section with a comparison of the objective in this paper vis-a-vis a similar objective considered by Tan et al in [12]. Broadly speaking, the core objectives in this paper and in [12] are similar, namely, restoring an electrical distribution network in the most efficient manner after a disaster. Specifically, in [12], the objective is to $\text{minimize}\,\sum_{l\in L^{D}}w_{l}T_{l}$, where $w_{l}$ is a weight on line $l$ representing its ‘importance’ (possibly related to the number of customers served by its downstream node) and $T_{l}$ is the energization time of line $l$. The objective function can therefore be interpreted as a measure of aggregate inconvenience, which is minimized. Inherent in this framework is the assumption that line $l$ will be energized after repairs as soon as it is feasible to do so, i.e., as soon as a path is established between $l$ and the source. The interpretation of $T_{l}$ is therefore it is the earliest possible timeline $l$ can be energized along the restoration timeline. In our framework, however, even if it is feasible to energize a line in the middle of the current restoration window after completion of repairs, we wait till the end of the current restoration windows, at which time, all lines which can be energized are energized, after necessary operational/safety checks. Restoration is therefore a staged process in our framework and broad swathes of the network are energized at certain planned restoration horizons, as opposed to a ‘restore whenever feasible’ framework in [12].

IV Problem Formulation

Figure 2: (a) Graph of an example radial distribution network, $G=(\mathcal{N},\mathcal{E})$, where $\mathcal{N}$ is the set of nodes and $\mathcal{E}$ is the set of edges. Edges $(1,2)$ and $(1,3)$ have been repaired and power has (possibly) been restored to nodes $2$ and $3$ at the start of the $k^{th}$ restoration window. The set of damaged edges yet to be repaired is therefore: $L^{D}=\{(2,4),(2,5),(3,6),(3,7),(4,8),(7,9)\}$. Each damaged edge represents a ‘job/repair site’. The estimated repair time and reward of the edge $(i,j)$ are $w_{ij}$ and $r_{ij}$ respectively.
(b) Graph $G$ modified by contracting the repaired/energized portion in $G$ into a supernode (root), $r$, and adding an auxiliary node, $0$, with an edge connecting $0$ and $r$. We denote this graph by $G_{m}$. The repair time and reward of the dummy edge $(0,r)$, $w_{0r}$ and $r_{0r}$, are both set to $0$.
(c) Line graph of $G_{m}$, which we denote by $L(G_{m})$. In this figure, we adopt a two-index notation for the node labels in the line graph simply for clarity. Repair times and rewards on edges in $G_{m}$ now appear as weights and rewards on nodes in $L(G_{m})$.
(d) The graph $L(G_{m})$ converted to a complete symmetric directed graph, which we denote by $L_{K}(G_{m})$. The weight of the edge $(i,j)\in L_{K}(G_{m})$, $c_{ij}$, represents the shortest travel time between job sites $i$ and $j$, $c_{ij}=c_{ji}$. Travel times between the node $(0,r)$ and any other node are set to $0$. We refer to the node $(0,r)$ as the root of $L_{K}(G_{m})$.

Consider the radial electrical distribution network, $G$, shown in Figure 2(a). The dotted portion of the network is assumed to have been repaired and energized at restoration horizon $T_{k-1}$ and the rest of the edges (job sites) in the network constitute the set $L^{D}$. Let $r_{ij}$ denote the reward associated with line $(i,j)$ and $w_{ij}$ denote its estimated mean repair time. If no repair is needed on line $(i,j)$, both the reward and repair time parameters are set equal to zero. In a practical setting, time overruns are possible due to unforeseen eventualities and it could well be the case that a repair which started in a previous restoration window and was scheduled for completion within that window may not actually have been finished. To accommodate such time overruns, we interpret $w_{ij}$ as the estimated mean residual repair time of line $(i,j)$. In this paper, we concentrate on single source radial networks (more generally, acyclic graphs), but it is possible to extend our approach to any network with an arbitrary number of sources. We assume that a transportation graph is also available (not shown in Figure 2), which allows us to determine the shortest travel times between all pairs of job sites. Additionally, as mentioned in Section III, we assume that the travel times are symmetric and it is possible to reach any site from any other site. In order to embed the repair times (associated with the power distribution network) and the travel times (associated with the transportation network) within the same graph, we proceed as follows:

• Step-$1$: Contract the already energized portion of the network into a supernode (root), add an auxiliary node (node $0$), and add an undirected edge between the auxiliary node and the root. Assign to this edge the parameter vector $\left(w_{0r}=0,r_{0r}=0\right)$, where $w_{ij}$ and $r_{ij}$ are the estimated repair time and reward of the edge $(i,j)$. The modified graph is denoted by $G_{m}$, as shown in Figure 2(b).

• Step-$2$: Construct the line graph of $G_{m}$, which we denote by $L(G_{m})$. Note that the job sites appear as nodes in $L(G_{m})$. Furthermore, each node in $L(G_{m})$ is endowed with two parameters, a reward parameter and a repair time parameter, as shown in Figure 2(c). Both these parameters are zero if a node does not need repair.

• Step-$3$: Convert the line graph $L(G_{m})$ to a complete graph, which we denote by $L_{K}(G_{m})$; $L(G_{m})\subseteq L_{K}(G_{m})$. This step is shown in Figure 2(d). We will refer to the node labeled $(0,r)$, interpreted as a dummy job, as the root of $L_{K}(G_{m})$. The repair time and reward parameters associated with each node in $L(G_{m})$ carry over to $L_{K}(G_{m})$. We will now assign weights to each undirected edge in $L_{K}(G_{m})$. Let $c_{ij}=c_{ji}$, representing the shortest travel time between nodes $i$ and $j$, denote the weight of the edge $(i,j)\in L_{K}(G_{m})$. For example, the weight of the undirected edge between the nodes $(4,8)$ and $(7,9)$ in Figure 2(d) represents the travel time for a crew moving from the line/job $(4,8)$ in Figure 2(a) to the line/job $(7,9)$ or vice versa. Travel times between the node $(0,r)$ and any other node are set to $0$, the reason for which will be apparent shortly.

A couple of observations are in order here. First, expansion of $L(G_{m})$ to $L_{K}(G_{m})$ in Step-$3$ allows the repair times (appears as node weights in Figure 2(d)) and the travel times (appears as edge weights in Figure 2(d)) to be embedded within the same undirected doubly weighted graph. If the number of damaged edges in Figure 2(a) is $|L^{D}|$, we end up with a doubly weighted complete graph in Figure 2(d) on $|L^{D}|+1$ nodes. Second, while the presence of each edge in $L(G_{m})$ reflects the physical connectivity of the underlying power distribution network $G$ (e.g., an edge exists between nodes $(r,4)$ and $(4,8)$ in Figure 2(c) since the undirected edges $(r,4)$ and $(4,8)$ share a node in Figure 2(b)), the interpretation of each edge in $L_{K}(G_{m})$ is the time needed to travel the shortest path between the end nodes of that edge. This observation will be important when we discuss the nature of feasible solutions next.

It so turns out that the optimal repair schedules can be interpreted in the context of a new variant of the traveling salesman problem. Assume that the number of repair crews is $m=1$ (one ‘repair salesman’). Operating within a certain time budget, the repair salesman starts from the node/home city $(0,r)$ in Figure 2(d) and selectively chooses to visit certain cities (repair sites) such that the reward collected at the end of the tour is maximized. If the salesman decides to visit a certain city, he must spend a prescribed amount of time at that city, equal to the repair time of that job site, before moving to the next city following the shortest path between the two cities. Since the travel time from the home city to any other node is $0$, our repair salesman does not incur any cost while traveling to the first repair site, which is consistent with our assumption (see Section III) that any restoration window commences when all repair crew are at their first job sites. The total cost incurred by the salesman on his tour is the sum of the repair times (node weights in $L_{K}(G_{m})$) and travel times (edge weights in $L_{K}(G_{m})$) and the aggregate reward is the sum of the individual rewards collected at the cities visited. Since every non-root node in $L_{K}(G_{m})$ is connected to the root by an edge with weight (travel time) $0$, the salesman can make his way back to his home city if visiting an additional city would result in a violation of the time budget, $\text{min}\left(\mathcal{T},\Delta T_{k}\right)$, where $\mathcal{T}$ represents the crew’s time budget and $\Delta T_{k}:=T_{k}-T_{k-1}$ denotes the length of the current restoration window.

Let us now take a look at some example ‘tours’ on $L_{K}(G_{m})$ (Figure 2(d)) which will help illustrate the nature of solutions we are looking for, as well as the similarities and differences between the traveling salesman framework outlined in the previous paragraph and the repair scheduling problem. For simplicity, we will assume that every node/city in $L_{K}(G_{m})$, other than the root, corresponds to a repair site and offers an unit reward.

• Suppose the repair salesman’s solution is $(0,r)\to(4,8)\to(r,4)\to(r,6)\to(0,r)$. This is a valid tour from a TSP standpoint since there are no subtours in the solution. After removing the first and last cities from the tour (doing so does not alter the cost of the tour since the weight of any edge connecting to the home city $(0,r)$ is zero), we obtain the repair sequence, $(4,8)\to(2,4)\to(3,6)$, in the context of Figure 2(a). This repair schedule is also valid from an electrical standpoint since a connected path exists between the source node $(0,r)$ and the repaired nodes $(4,8)$, $(r,4)$, and $(r,6)$ in $L(G_{m})$ (Figure 2(c)). Note that the same conclusion could have been reached from either Figure 2(a) or Figure 2(b). Since all three distribution lines $(2,4)$, $(4,8)$, and $(3,6)$ in Figure 2(a) can be energized at the end of the current restoration window, the salesman can claim a reward of 3 units.

• Suppose the repair salesman’s solution is $(0,r)\to(r,4)\to(4,8)\to(7,9)\to(0,r)$, or equivalently, the repair schedule $(2,4)\to(4,8)\to(7,9)$ in the context of Figure 2(a). Although this is a valid solution from a TSP perspective, this is not an electrically valid solution since the set of repaired nodes, $(0,r)$, $(r,4)$, $(4,8)$, and $(7,9)$, do not induce a connected subgraph in $L(G_{m})$ (Figure 2(c)). Alternately, from Figure 2(a), we can observe that the distribution line $(7,9)$ cannot be energized unless line $(3,7)$ is repaired. We say that this solution satisfies the path continuity condition (no subtours) but violates the electrical continuity condition.

• Suppose the repair salesman’s solution is $(0,r)\to(r,4)\to(4,8)\to(0,r),(r,5)\to(r,6)\to(r,7)\to(r,5)$. Clearly, this is not a valid solution from a TSP perspective since it has subtours, although it is electrically valid since the set of repaired nodes, $(0,r)$, $(r,4)$, $(4,8)$, $(r,5)$, $(r,6)$, and $(r,7)$, induce a connected subgraph in $L(G_{m})$. Therefore, this solution satisfies the electrical continuity condition but violates the path continuity condition.

• Suppose there are two repair salesmen (both departing from the home city, operating under their individual time constraints) and their solutions are $(0,r)\to(4,8)\to(7,9)\to(0,r)$ and $(0,r)\to(r,4)\to(r,7)\to(0,r)$. Both solutions are free from subtours and therefore valid from a TSP perspective. These solutions are also electrically valid since the aggregate set of repaired nodes induce a connected subgraph in $L(G_{m})$ (although, the first salesman’s solution, by itself, is not electrically valid). We say that the two solutions satisfy path continuity constraints (individually) as well as the electrical continuity constraint (aggregated). The total number of rewards claimed by the two salesmen is therefore four.

The previous examples illustrate the two types of continuity constraints, path continuity (alternately, subtour elimination constraint) and electrical continuity, which we require for a repair schedule to be valid. Next, we will discuss how electrical continuity constraints may be imposed. Modeling such constraints is dependent on the specific topology of the power distribution network. We show that for single-source radial distribution networks, such constraints translate to precedence relations between the set of repaired nodes.

Figure 3: (a) Graph of the radial distribution network shown in Figure 2. Arrow directions on the edges denote the direction of power flow. Since this is a directed graph, we denote it by $G^{d}$. (b) Graph $G^{d}$ modified by contracting the repaired/energized portion in $G^{d}$ into a supernode, $r$, and adding an auxiliary node, $0$, connected to node $r$. We denote this directed graph by $G_{m}^{d}$. (c) Line graph of the directed graph shown in panel (b), $L(G_{m}^{d})$.

Consider the power distribution network shown in Figure 3(a), which is the same network shown in Figure 2(a), except arrows have been added indicating the direction of power flow. We denote this directed graph by $G^{d}$. Next, we contract the repaired/energized portion in $G^{d}$ into a supernode, $r$, and add an auxiliary node, $0$, connected to node $r$, as shown in Figure 3(b). We denote this modified graph by $G_{m}^{d}$. Finally, we convert $G_{m}^{d}$ to its line graph representation, $L(G_{m}^{d})$, as shown in Figure 3(c). This directed line graph captures all the precedence relations we will require our solution to satisfy. For example, if node $(7,9)$ is chosen to be repaired during the current window, we will require that node $(r,7)$ be repaired first. Similarly, if node $(4,8)$ is chosen to be repaired during the current window, we will require that node $(r,4)$ be repaired first. No precedence relations are required for nodes $(r,4)$, $(r,5)$, $(r,6)$, and $(r,7)$. It is easy to verify that if these precedence relations are enforced, the solution $(0,r)\to(r,4)\to(4,8)\to(7,9)\to(0,r)$ (which violates the electrical continuity constraint) will be rendered infeasible.

IV-A Electrical continuity constraints for partially damaged networks

In the previous paragraph, we assumed that all nodes in Figure 3(c), other than the root node $(0,r)$, are damaged. For partially damaged scenarios, we adopt one more pre-processing step whereby the line graph $L(G_{m}^{d})$ is “collapsed” to eliminate all non-damaged nodes. In order to avoid additional notation, we will continue to refer to the collapsed graph as $L(G_{m}^{d})$, i.e., $L(G_{m}^{d})\to\text{collapse procedure}\to L(G_{m}^{d})$. Precedence relations are extracted as part of this procedure.

Consider the directed graph $L(G_{m}^{d})$ shown in Figure 4(a), where we have adopted single-index node labels for simplicity, unlike the node labels in Figure 3. Figure 4(b) shows the modified $L(G_{m}^{d})$, which is obtained as follows:

• Step-$1$: First, all undamaged leaf nodes in $L(G_{m}^{d})$ are deleted, along with the directed edges incident on these nodes. Deletion of an undamaged leaf node may create an additional undamaged leaf node, which is also deleted. At the end of this step, there must not be any undamaged leaf node in $L(G_{m}^{d})$. Executing this step on Figure 4(a) leads to the deletion of nodes labeled 10, 11, 6, and 1. Note that the last three nodes are deleted sequentially (after node 11 is deleted, node 6 becomes a leaf node, which is deleted next, followed by node 1).

• Step-$2$: Next, for every damaged node $j\in L(G_{m}^{d})$, we search for its nearest ancestor which is also damaged. If no such node is found, we record its nearest ancestor to be the root node (which, by definition, is undamaged). The precedence relation for node $j$ is therefore of the form $i\to j$, where $i$ is the nearest ancestor of $j$ which is also damaged. If $i=\emptyset$, we assign $i:=\text{root node of $L(G_{m}^{d})$}$. Executing this step on Figure 4(a) leads to the following precedence relations: (i) $0\to 7$ (node 2 deleted), (ii) $0\to 3$, (iii) $3\to 12$ (node 8 deleted), (iv) $0\to 4$, (v) $4\to 14$ (nodes 9 and 13 deleted), and (ii) $0\to 5$. These precedence relations are embodied in Figure 4(b).

The previous discussion leads to an interesting question: ‘for a partially damaged network, how should one assign the rewards for a damaged node?’ For example, referring to Figure 4(a), we see that if node 4 is repaired and energized, it may be possible to energize nodes 9 and 13 simultaneously, provided there is a switch between nodes 13 and 14. In this case, the true reward for repairing node 4 is the sum of the individual rewards associated with nodes 4, 9, and 13. Accounting for node rewards in this manner, aided by knowledge of switch locations, may lead to improved decision making. A similar argument in favor of switch location guided repair scheduling can also be made for a fully damaged network. We plan on investigating this issue in future.

Remark IV.1

In order to keep the set of nodes in $L_{K}(G_{m})$ (see Figure 2(d)) and $L(G_{m}^{d})$ identical, nodes which are deleted as part of Step-$2$ above should also be deleted (along with incident edges) from $L_{K}(G_{m})$.

Figure 4: Illustrating the “graph collapse” procedure for a partially damaged network. (a) $L(G_{m}^{d})$ as obtained after application of the procedure illustrated in Figure 3. (b) $L(G_{m}^{d})$ after the undamaged nodes (other than the root) are deleted.

We are now in a position to state the MILP formulation for the $m$-crew repair scheduling problem. Although we have cycled through different graph representations, our model will be defined on the undirected complete graph $L_{K}(G_{m})$ (see Figure 2(d)) and the directed graph $L(G_{m}^{d})$ (see Figure 3(c)). We will use $L_{K}(G_{m})$ to set up the repair scheduling problem as a traveling salesman-like problem and $L(G_{m}^{d})$ to impose the additional electrical continuity constraints. For the rest of the paper, we will refer to these two graphs as working graphs and denote them by $G_{w}$ and $G_{wd}$ respectively, the subscript $w$ indicating ‘working’ and the subscript $d$ indicating ‘directed’. Variables defined on undirected edges will bear the subscript $(ij)$; e.g., $x_{(ij)c}=1$ would mean that the undirected edge between $i$ and $j$ is traveled by crew $c$. Variables defined on directed edges will bear the subscript $ij$; e.g., $x_{ijc}=1$ would mean that the directed edge from $i$ to $j$ is traveled by crew $c$. Table I summarizes the major notation we will use in the rest of the paper. We also summarize all our assumptions below:

1. The power distribution network is radial with a single source. Our assumption of a radial distribution network isn’t too restrictive since, even if such networks are topologically non-radial (e.g., meshed), they are operated radially. If operational decisions are made prior to repair planning/scheduling, allowing the non-radial network to be pruned and pre-configured to a radial structure, the MILP model discussed in this paper will be applicable.

2. Every repair site is reachable from every other repair site.

3. The travel time between any two repair sites, following the shortest path between the two sites, is symmetric.

4. The reward parameters (possibly determined by utilities or other planning agencies) are known.

5. Time accounting for every crew starts when they are at their first repair site. That is, we do not assume that crews are dispatched from a central staging location, in which case additional travel time needs to be accounted for while commuting from the staging location to the first repair site.

Notation	Description
$G_{w}$	undirected complete working graph, $L_{K}(G_{m})$ (see Figure 2(d)).
$\mathcal{E}_{w}$	set of undirected edges in $G_{w}$.
$G_{wd}$	directed working graph, $L(G_{m}^{d})$ (see Figure 3(c)).
$\mathcal{E}_{wd}$	set of directed edges in $G_{wd}$.
$\mathcal{N}_{w}$	set of all nodes in $G_{w}$ or $G_{wd}$; $\mathcal{N}_{w}=\{0,1,2,\ldots n-1\}$, where node $0$ is the root/home city.
$n$	number of nodes in $G_{w}$ or $G_{wd}$
$\mathcal{N}_{w}^{D}$	subset of nodes in $\mathcal{N}_{w}$ which are damaged and require repair; note that node 0 is not in $\mathcal{N}_{w}^{D}$.
$r_{i}$ $(\geq 0)$	reward of node $i\in\mathcal{N}_{w}$; $r_{i}=0$, $\forall i\in\left\{\mathcal{N}_{w}\setminus\mathcal{N}_{w}^{D}\right\}$.
$w_{i}$	repair time of node $i\in\mathcal{N}_{w}$; $w_{i}=0$, $\forall i\in\left\{\mathcal{N}_{w}\setminus\mathcal{N}_{w}^{D}\right\}$.
$c_{(ij)}$	travel time between nodes $i,j\in\mathcal{N}_{w}$, following the shortest path between $i$ and $j$; i.e., $c_{(ij)}=c_{ij}=c_{ji}$. If either $i=0$ or $j=0$, then $c_{(ij)}=0$
$\mathcal{C}$	set of crew indices; $\mathcal{C}=\{1,2,\ldots m\}$
$m$	number of repair crew (repair salesmen)
$\mathcal{T}_{c}$	time budget of crew $c$
$\Delta T_{k}$	time duration of the current ($k^{th}$) restoration window.
$y_{ic}$	binary variable which takes the value $1$ if crew $c$ visits/repairs node $i\in\mathcal{N}_{w}^{D}$ or node $i$ does not need repair, i.e., $i\in\left\{\mathcal{N}_{w}\setminus\left\{\mathcal{N}_{w}^{D}\cup 0\right\}\right\}$, $0$ otherwise.
$x_{(ij)c}$	integer variable which takes the value $0$ if crew $c$ does not travel between nodes $i,j\in\mathcal{N}_{w}$, $i\neq j$. If crew $c$ travels between nodes $i$ and $j$, then $x_{(ij)c}\in\{1,2\}$, provided the undirected edge $(i,j)$ is between node $0$ and some other node. If $(i,j)$ is not incident on the root, $x_{(ij)c}=1$.
$(i,j)$	undirected edge between nodes $i$ and $j$.
$i\to j$	directed edge from node $i$ to $j$.
$\{i,j\}$	set of nodes $i$ and $j$.

TABLE I: Table of notations used in MILP modeling of the repair scheduling problem.

V MILP formulation

In this section, we discuss a MILP model for the repair scheduling problem. Given a distribution network $G$, we first perform a series of pre-processing steps to create the two working graphs, $G_{w}$ and $G_{wd}$ (see Table I), summarized below:

$$G\,\to\, \,G_{m}\,\to\, L(G_{m})\,\to\,\overbrace{L_{K}(G_{m})}^{\text{Figure~{}\ref{fig:Fig1_2}(d)}} :=G_{w}=\left(\mathcal{N}_{w},\mathcal{E}_{w}\right)$$

$$\downarrow$$

$$G^{d}\,\to\, \,G_{m}^{d}\to \underbrace{L(G_{m}^{d})}_{\text{Figure~{}\ref{fig:Fig5}(c)}} :=G_{wd}=\left(\mathcal{N}_{w},\mathcal{E}_{wd}\right)$$

Without any loss of generality, we assume that all nodes in $G_{w}$ and $G_{wd}$, other than the root node, are damaged. For a partially damaged scenario, this requires the “graph collapse” procedure explained in Section IV-A to be implemented (see also Remark IV.1). As mentioned previously, we will use $G_{w}$ to set up the repair scheduling problem as a traveling salesman-like problem and $G_{wd}$ to impose the additional electrical continuity constraints. Note that $\left|\mathcal{E}_{w}\right|=n(n-1)/2$ since $G_{w}$ is an undirected graph. Accordingly, we will choose the set of edges in $G_{w}$ to be $\mathcal{E}_{w}=\{(i,j):0\leq i\leq n-2,\,j>i\}$. Our decision variables are $\left\{y_{ic}:i\in\mathcal{N}_{w},\,c\in\mathcal{C}\right\}$ and $\left\{x_{(ij)c}:(i,j)\in\mathcal{E}_{w},\,c\in\mathcal{C}\right\}$. The interpretations of these variables are discussed in Table I.

The MILP formulation for $m$-crew repair scheduling is listed in Figure 5. We dispatch $m$ ‘repair salesmen’, each salesman starting and ending his tour at the root node $0$. Additionally, each salesman operates under his own time budget and each node which needs repair can be visited by at most one repair salesman. Our objective function (eqn. (1a)) is maximization of the aggregate reward that can be collected by all crew at the end of the current restoration window. If $r_{i}=1$, for all $i\in\mathcal{N}_{w}^{D}$, the objective function reduces to maximization of the number of nodes that can be repaired within the current restoration window. The roles of the constraint equations are discussed below.

• Eqn. (1b) forces each crew to visit the root (node $0$).

• Eqn. (1c) imposes binary conditions on all nodes in $\mathcal{N}_{w}^{D}$ (set of all nodes which need repair), on a per crew basis.

• Eqns. (1d) and (1e) impose integrality conditions on the triple indexed $x_{(ij)c}$ variables. The set of values these variables can take depends on whether the undirected edge $(i,j)$ is incident on the root node or not. If $(i,j)$ is not incident on the root, $x_{(ij)c}$ is binary valued (eqn. (1e)). In contrast, if $(i,j)$ is incident on the root, $x_{(ij)c}$ takes a value from the set $\{0,1,2\}$ (eqn. (1d)). We allow $x_{(ij)c}$ to take the value $2$ in the latter case since we accept cycles of length $2$ in the optimal solution. For example, referring to Figure 2(d), the optimal solution for $m=2$ could be $(0,r)\to(r,4)\to(0,r)$ and $(0,r)\to(r,5)\to(r,6)\to(0,r)$. In this case, the undirected edge between $(0,r)$ and $(r,4)$ is used twice by the first crew, but all edges traveled by the second crew are used once.

• Eqn. (1f) stipulates that all nodes needing repair be visited/repaired by at most one crew.

• Eqn. (1g) relates the $x$ variables to the $y$ variables, on a per crew basis. If crew $c$ visits node $i$, where $i$ is either the root node or belongs to the set of nodes needing repair, then $y_{ic}=1$, which implies $\sum_{j:(ij)\,\in\,\mathcal{E}_{w}}x_{(ij)c}=2$. This leaves open two possibilities: (i) some undirected edge incident on node $i$ is used twice by the crew, which happens only when the edge connects to the root and the optimal solution for that crew is a cycle of length $2$ (i.e., the tour is of the form $0\to i\to 0$), or (ii) two different edges incident on $i$ are used, once each, which happens when the optimal solution for that crew is a cycle of length greater than $2$, in which case one edge is used to enter the node and another to exit (i.e., the tour is of the form $0\to i\to q\cdots\to 0$ or $0\to\cdots\to p\to i\to q\cdots\to 0$). The role of this constraint, when applied to the home node, i.e., node $0$, in determining model feasibility under limited time budgets is discussed in Section V-A.

• Eqn. (1h) imposes a time constraint on every tour, on a per crew basis. For crew index $c$, the maximum tour cost (in units of time) is the minimum of $\mathcal{T}_{c}$ and $\Delta T_{k}$, where $\mathcal{T}_{c}$ is the time budget for crew $c$ and $\Delta T_{k}$ is the time duration of the $k^{th}$ restoration window. Since we do not allow fractional repairs (and therefore, prorated rewards), the l.h.s of eqn. (1h) may not equal the r.h.s in the optimal solution. Note that both repair times and travel times are accounted for in the l.h.s of eqn. (1h).

• Eqn. (1i) enforces the electrical continuity constraint by ensuring that for every node which is visited/repaired by some crew, a directed fully repaired path exists from the source to that node. This is achieved by imposing the complete set of precedence relations dictated by the working graph $G_{wd}$ (i.e., $L(G_{m}^{d})$; see Figure 3(c), Figure 4 for a partially damaged network, and Remark IV.1).

  In essence, the electrical continuity constraints yield a pattern of energized swathes which expand radially outward from the source over time, which makes intuitive sense for single-source radial distribution networks. For other topologies with either single or multiple sources, these constraints may need to be modeled differently.

• Eqn. (1j) is a placeholder for subtour elimination constraints (SECs) which we address below.

 

max $$\sum_{i\,\in\,\mathcal{N}_{w}^{D}}r_{i}\,\sum_{c\,\in\,\mathcal{C}}y_{ic}$$ (1a)

s.t

$$y_{ic} =1; \forall c\in\mathcal{C},\ i=0$$ (1b)

$$y_{ic} \in\{0,1\}; \forall c\in\mathcal{C},\ \forall i\in\mathcal{N}_{w}^{D}$$ (1c)

$$x_{(ij)c} \in\{0,1,2\};\ \forall c\in\mathcal{C},\ \forall(i,j)\,\in\,\mathcal{E}_{w},\ i=0$$ (1d)

$$x_{(ij)c} \in\{0,1\}; \forall c\in\mathcal{C},\ \forall(i,j)\,\in\,\mathcal{E}_{w},\ i\neq 0$$ (1e)

$$\sum_{c\,\in\,\mathcal{C}} y_{ic} \leq 1;\quad \forall i\in\mathcal{N}_{w}^{D}\ (\text{note: }0\not\in\mathcal{N}_{w}^{D})$$ (1f)

$$\sum_{j:(ij)\,\in\,\mathcal{E}_{w}} x_{(ij)c}-2\,y_{ic} =0; \forall c\in\mathcal{C},\ \forall i\in\left\{\mathcal{N}_{w}^{D}\cup 0\right\}\ \ \text{/* see eqns.~{}(\ref{eq:const1}), (\ref{eq:const2}) */}$$ (1g)

$$\sum_{i\,\in\,\mathcal{N}_{w}^{D}} w_{i}\,y_{ic}+\sum_{(i,j)\,\in\,\mathcal{E}_{w}}c_{(ij)}\,x_{(ij)c} \leq\text{min}\left(\mathcal{T}_{c},\Delta T_{k}\right);\quad \forall c\in\mathcal{C}$$ (1h)

$$\sum_{c\,\in\,\mathcal{C}} y_{ic}-\sum_{c\,\in\,\mathcal{C}}y_{jc} \geq 0;\quad \forall(i\to j)\in\mathcal{E}_{wd}\ \ \mbox{/* electrical continuity constraints */}$$ (1i)

subtour elimination constraints (1j)

 

Figure 5: MILP formulation for $m$-crew repair scheduling during the $k^{th}$ restoration window. Note that the electrical continuity constraints, eqn. (1i), are valid only for single source, radial power distribution networks.

V-A Discussion on problem infeasibility due to limited time budgets

In this section, we discuss why the MILP model can produce counterintuitive results with respect to the number of crews $m$, in the sense that the problem can be feasible for some values of $m$ and $m+2$, but infeasible for $m+1$. To illustrate why and when this might happen, consider the distribution network $1\rightarrow 2\rightarrow 3\rightarrow 4$, where node $1$ is the source (undamaged) and the directed arrows represent the direction of power flow. For simplicity, assume all travel times are zero, and the repair times of nodes 2, 3, and 4 are 30 minutes, 30 minutes, and 20 minutes respectively. Suppose we can choose to deploy one, two, or three crews, and the time budgets for the first and third crews are 30 minutes each, while that for the second crew is 20 minutes. When only the first crew is available, it can be deployed to node 2. Similarly, when all three crews are available, the first and third crews can be deployed to nodes 2 and 3, while the second crew can be deployed to node 4. However, if only the first two crews are available, the first crew can be deployed to node 2, but no job can be assigned to the second crew (electrical continuity constraints prevent deployment of the second crew to node 4 until node 3 is repaired), rendering the model infeasible. Whenever we encountered such infeasibilities during simulations, we discarded those problem instances.

Although not implemented, there exist at least two approaches by which model infeasibility due to inadequate time budgets can be avoided. Before we discuss possible solutions, let us first understand why eqn. 1g:

$$\sum_{j:(ij)\,\in\,\mathcal{E}_{w}}\ x_{(ij)c}-2\,y_{ic}=0;\ \forall c\in\mathcal{C},\ \forall i\in\left\{\mathcal{N}_{w}^{D}\cup 0\right\},$$

is responsible for inducing such infeasibilities. Since $y_{0c}=1$ for all $c$, $\sum_{j:(0j)\,\in\,\mathcal{E}_{w}}\ x_{(0j)c}=2$, which forces every crew to venture out of the home node, even if its time budget isn’t enough to finish a single job meeting the electrical continuity constraints. This explains why the 2-crew instance in the previous example is infeasible. The first approach to circumventing this infeasibility issue involves introducing a dummy node, say node 0, so that the new network is $0-1\rightarrow 2\rightarrow 3\rightarrow 4$. The reason for showing the edge between nodes 0 and 1 with an undirected line is that, while node 0 plays no role as far as power flow (electrical continuity) is concerned, the edge $0-1$ is available for travel between nodes 0 and 1. Now, if we set a zero reward for node 0 and zero travel time on the undirected edge $0-1$, the second crew in the above example is free to make the tour $1\to 0\to 1$, for zero reward and cost, rendering a 2-crew solution feasible. The other approach involves splitting eqn. 1g into two:

i) $$\sum_{j:(ij)\,\in\,\mathcal{E}_{w}}\ x_{(ij)c}-2\,y_{ic} =0;\ \forall c\in\mathcal{C},\ \forall i\in\mathcal{N}_{w}^{D}$$

ii) $$\sum_{j:(ij)\,\in\,\mathcal{E}_{w}}\ x_{(ij)c}-2\,y_{ic} \leq 0;\ \forall c\in\mathcal{C},\ i=0$$

Notice that the first equation (equality) is applicable for the damaged nodes only, while the second equation (inequality) is for the home node (undamaged). Since the equality is loosened to an inequality for the home node, $\sum_{j:(0j)\,\in\,\mathcal{E}_{w}}\ x_{(0j)c}\leq 2$, and the second crew in our above example can choose not to travel at all from the home node. However, if the crew does travel from the home node to some damaged node $i$, the equality constraint in the first equation will force it to either travel to another damaged node $j$ (the tour will be of the form $0\to i\to j\to\ldots\to 0$) or come back to the home node (the tour will be of the form $0\to i\to 0$), depending on the remaining time.

V-B Subtour elimination constraints

There exists a rich body of literature on different approaches for subtour elimination in the traveling salesman problem (TSP). Some of these approaches utilize an exponential number of constraints, while others are polynomial formulations. It is not our intent to review relevant existing literature in this paper. In our work, we have adopted and modified the single commodity flow (SCF) based subtour elimination procedure proposed by Gavish and Graves [28]. In the context of the TSP, the SCF-based subtour elimination procedure requires that the home city be able to send a unit of flow to each of the cities visited, under the constraint that a positive flow can be activated on the edge $i\to j$ only if salesman has traveled from $i$ to $j$. Clearly, for the flow constraints to be satisfied, there must not be any subtour in the optimal solution.

To accommodate a flow model, we view each undirected edge $(i,j)\in\mathcal{E}_{w}$ as a bidirected edge with corresponding (continuous-valued and non-negative) flow variables $f_{ij}$ and $f_{ji}$, where $f_{ij}$ is interpreted as a flow from $i$ to $j$. In other words, each undirected edge $(i,j)\in\mathcal{E}_{w}$ maps to two flow variables, $f_{ij}$ and $f_{ji}$. We define the ‘flow differential’ at any node as the aggregate flow out of the node minus the aggregate flow into the node. Our SCF based subtour elimination formulation is shown in eqn. (3), which requires an additional $O(n^{2})$ continuous variables and $O(n^{2})$ constraints. Notice that all inequalities are coupling constraints, in the sense that every equality/inequality depends on the collective job completions over all $m$ crew.

Eqn. (3a) ensures that the flow differential at the root is equal to the total number of jobs completed by all crew. Eqn. (3b) ensures that the flow differential at any node which needs repair is equal to $-1$ if it has been visited/repaired by any crew, and $0$ otherwise.

Eqn. (3c) is a flow activation condition which ensures that flows on $i\to j$ and $j\to i$ can be activated only if some crew has traveled the undirected edge $(i,j)$, i.e., $\sum_{c\,\in\,\mathcal{C}}x_{(ij)c}\geq 1$. If flows are activated, it also imposes a capacity constraint which limits the maximum flow to either $2mn$ or $2n^{2}$. To see why, consider some node $i\in\left\{\mathcal{N}_{w}^{D}\cup 0\right\}$. Summing eqn. (1g) over $c$, we find:

$$\sum_{c\,\in\,\mathcal{C}}\sum_{j:(ij)\,\in\,\mathcal{E}_{w}}x_{(ij)c}\leq 2\sum_{c\,\in\,\mathcal{C}}y_{ic},\ \forall i\in\left\{\mathcal{N}_{w}^{D}\cup 0\right\}$$ (2)

If $i=0$, the r.h.s of eqn. (2) evaluates to $2m$, using eqn. (1b). Since:

$$\sum_{c\,\in\,\mathcal{C}}\sum_{j:(ij)\,\in\,\mathcal{E}_{w}}x_{(ij)c}\leq 2m\ \Rightarrow\ \sum_{c\,\in\,\mathcal{C}}x_{(ij)c}\leq 2m,$$

eqn. (3c) reduces to:

$$f_{ij},\,f_{ji}\leq n\sum_{c\,\in\,\mathcal{C}}x_{(ij)c}=2mn$$

However, if $i\in\mathcal{N}_{w}^{D}$, the r.h.s of eqn. (2) evaluates to $2n$, using eqn. (1f), and eqn. (3c) reduces to:

$$f_{ij},\,f_{ji}\leq n\sum_{c\,\in\,\mathcal{C}}x_{(ij)c}=2n^{2}$$

Eqn. (3d) significantly tightens the flow upper bound established by eqn. (3c) by restricting the maximum flow on any directed edge to the total number of nodes visited/repaired by all crew, which can be at most equal to the cardinality of $\mathcal{N}_{w}^{D}$ (i.e., number of nodes needing repair). Although a tighter upper bound would have been $\mbox{max}_{c\,\in\,\mathcal{C}}\left(\sum_{i\,\in\,\mathcal{N}_{w}^{D}}y_{ic}\right)$, i.e., the maximum number of nodes visited/repaired by any crew, we did not implement this in order to avoid introducing additional auxiliary binary variables.

 

$$\sum_{j=1}^{n}f_{0j}-\sum_{j=1}^{n}f_{j0} =\sum_{c\,\in\,\mathcal{C}}\,\sum_{i\,\in\,\mathcal{N}_{w}^{D}}y_{ic};$$ (3a)

$$\sum_{\begin{subarray}{c}j=0\\ j\neq i\end{subarray}}^{n}f_{ij}-\sum_{\begin{subarray}{c}j=0\\ j\neq i\end{subarray}}^{n}f_{ji} =-\sum_{c=1}^{m}y_{ic};\ \forall i\in\mathcal{N}_{w}^{D}$$ (3b)

$$f_{ij},\,f_{ji} \leq n\sum_{c\,\in\,\mathcal{C}}x_{(ij)c};\ \forall(i,j)\in\mathcal{E}_{w}$$ (3c)

$$f_{ij},\,f_{ji} \leq\sum_{c\,\in\,\mathcal{C}}\sum_{i\,\in\,\mathcal{N}_{w}^{D}}y_{ic};\ \forall(i,j)\in\mathcal{E}_{w}$$ (3d)

$$f_{ij},\,f_{ji} \geq 0;\, \forall(i,j)\in\mathcal{E}_{w}$$ (3e)

 

Figure 6: Single commodity flow based subtour elimination constraints to be used in conjunction with the MILP model in Figure 5. Each undirected edge $(i,j)$ in $\mathcal{E}_{w}$ gives rise to two directed flow variables, $f_{ij}$ and $f_{ji}$.

V-C A special case: when $\left\{\mathcal{N}_{w}\setminus\left\{\mathcal{N}_{w}^{D}\cup 0\right\}\right\}=\emptyset$ and the power distribution network is a star graph: introducing a new variant of the TSP

An even more special case arises when the power distribution network is a star graph and $\left\{\mathcal{N}_{w}\setminus\left\{\mathcal{N}_{w}^{D}\cup 0\right\}\right\}=\emptyset$. Note that all star graphs are radial, but the converse is not true. Since the line graph of a star graph on $n$ nodes is a complete graph on $n-1$ nodes, the graph $L(G_{m})$ (see Figure 2(c)) itself would be a complete graph, i.e., $L(G_{m})=L_{K}(G_{m})$ (in general, $L(G_{m})\subseteq L_{K}(G_{m})$). Since the set of edges in $L(G_{m})$ and $L_{K}(G_{m})$ are identical, $L(G_{m})$ itself is a complete descriptor of the physical connectivity of the power network as well as the connectivity of the transportation graph. Another way of stating this is that, energization of any node after repair is not contingent on the repair status of any other node. In fact, the precedence relations embodied by $\mathcal{E}_{wd}$ would state:

$$\sum_{c\,\in\,\mathcal{C}}y_{0c}\geq\sum_{c\,\in\,\mathcal{C}}y_{ic};\quad\forall i\in\left\{\mathcal{N}_{w}\setminus 0\right\},$$

which will always be true since the l.h.s of the inequality evaluates to the number of crew, $m$ ($\geq 1$, see eqn. (1b)), and the r.h.s is $\leq 1$ (see eqn. (1f)). This proves that the electrical continuity constraint in eqn. (1i) can be dropped and the pre-processing steps outlined in Figure 3 can be omitted. See Figure 7 for an illustration.

The $m$-crew repair scheduling problem therefore reduces to: “given a complete doubly weighted (node costs as well as edge costs) graph on $n$ nodes, where each node offers a reward and each crew has a time budget, find optimal tours for $m$ crews such that the aggregate reward collected by all crews is maximized.”

Essentially, what we have arrived at is, to the best of our knowledge, a new variant of the multiple traveling salesman problem which has not been studied before, which we call the “cost constrained reward maximizing mTSP on doubly weighted graphs”, or CCRM-$m$TSP-DW. The variant of the TSP which is closest to CCRM-$m$TSP-DW is the “selective TSP (S-TSP)” which was introduced by Laporte and Martello [29]. The objective in the S-TSP is to visit a subset of $n$ cities such that the aggregate reward collected by the salesman at the end of the tour is maximized, under the constraint that the total cost incurred by the salesman be no more than some pre-determined amount. Therefore, the CCRM-$m$TSP-DW can be viewed as a generalization of S-TSP along two fronts: (i) multiple salesmen instead of one, and (ii) doubly weighted input graph with node weights and edge weights, instead of an edge-weighted graph. From this perspective, CCRM-$m$TSP-DW $\equiv$ S-$m$TSP-DW. Based on our preceding discussion, it follows that the MILP model in Figure 5, along with appropriately defined subtour elimination constraints (e.g., using the equations in Figure 6), can also be used to solve an instance of the CCRM-$m$TSP-DW problem, after setting the r.h.s of eqn. (1h) equal to $\mathcal{T}_{c}$ since $\Delta T_{k}$ plays no role in a TSP framework.

Figure 7: (a) Graph of a star connected distribution network, $G$. (b) Graph $G$ modified by contracting the repaired/energized portion in $G$ into a supernode (root), $r$, and adding an auxiliary node, $0$, with an edge connecting $0$ and $r$. We denote this graph by $G_{m}$. (c) Line graph of $G_{m}$, which we denote by $L(G_{m})$. Since $L(G_{m})$ is itself a complete graph, this is our undirected working graph $G_{w}$. (d) Directed version of $G$, $G^{d}$, showing the direction of power flow. (e) Graph $G^{d}$ modified by contracting the repaired/energized portion in $G^{d}$ into a supernode, $r$, and adding an auxiliary node, $0$, connected to node $r$. We denote this directed graph by $G_{m}^{d}$. (f) Line graph of the directed graph shown in panel (e), $L(G_{m}^{d})$. This is our directed working graph $G_{wd}$.

V-D Valid inequalities

The MILP model in Figure 5 can be fortified through addition of valid inequalities. One set of such inequalities can be: “if undirected edges $(i,j)$ and $(j,k)$ are traveled on by any crew, then $(i,k)$ must not be traveled on by that crew, where $i\neq j\neq k\neq 0,\,\left\{i,j,k\right\}\in\mathcal{N}_{w}$”.

$$\hskip 12.91663ptx_{(ij)c}+x_{(jk)c}+x_{(ik)c}\,\leq\,y_{ic}+y_{jc};\ \forall c\in\mathcal{C}$$ (4a)

$$\hskip 12.91663ptx_{(ij)c}+x_{(jk)c}+x_{(ik)c}\,\leq\,y_{jc}+y_{kc};\ \forall c\in\mathcal{C}$$ (4b)

$$\hskip 12.91663ptx_{(ij)c}+x_{(jk)c}+x_{(ik)c}\,\leq\,y_{ic}+y_{kc};\ \forall c\in\mathcal{C}$$ (4c)

The above system of inequalities ensures that:

• if $i$, $j$ and $k$ are all repaired by crew $c$, then $x_{(ij)c}+x_{(jk)c}+x_{(ik)c}\leq 2$.

• if any two from $i$, $j$ and $k$ are repaired, then $x_{(ij)c}+x_{(jk)c}+x_{(ik)c}\leq 1$, and

• if any one from $i$, $j$ and $k$ or none are repaired, then $x_{(ij)c}+x_{(jk)c}+x_{(ik)c}\leq 0$.

In essence, eqn. (4), which introduces $3\,\binom{n}{3}$ constraints, is equivalent to:

$$x_{(ij)c}+x_{(jk)c}+x_{(ik)c}\leq\text{max}\left(0,y_{ic}+y_{jc}+y_{kc}-1\right).$$

V-E Alternate objective function: profit maximization instead of reward maximization

Since the objective function in eqn. (1a) will always prefer nodes with higher rewards, repairs of lower priority nodes (low reward) will tend to be deferred to later restoration windows. This phenomenon will be most acute if there are no high priority nodes downstream of the low priority node along a power flow path from the root (see Figure 3). In order to mitigate this situation, we can introduce penalty parameters, $\{p_{i}:i\in\mathcal{N}_{w}^{D}\}$, which could be time dependent so that the penalty associated with a node, if not chosen for repair during the current restoration window, depends on the time elapsed since arrival of damage information for that node. Eqn. (1a) can then be modified to as follows:

$$\text{max }\left(\sum_{i\,\in\,\mathcal{N}_{w}^{D}}r_{i}\,\sum_{c\,\in\,\mathcal{C}}y_{ic}-\sum_{i\,\in\,\mathcal{N}_{w}^{D}}p_{i}\,\sum_{c\,\in\,\mathcal{C}}\left(1-y_{ic}\right)\right),$$ (5)

where the first term is the aggregate reward and the second term is the aggregate penalty. If $\left\{\mathcal{N}_{w}\setminus\left\{\mathcal{N}_{w}^{D}\cup 0\right\}\right\}=\emptyset$ and the power distribution network is a star graph, the resulting TSP variant can be dubbed the “cost constrained profit maximizing mTSP on doubly weighted graphs”, or CCPM-$m$TSP-DW.

VI Computational Results

In this section, we evaluate our base model shown in Figure 5 (i.e, without the valid inequalities discussed in Section V-D) on three IEEE test feeder networks, shown in Figures A1, A2 and A3 in the Appendix. Some modifications were made to the standard IEEE networks, which are discussed below:

• Link distances in the standard IEEE test networks [30], which are specified as 0 ft are set equal to 1 ft. For example, in the case of the 13-node network, the distance of link $(633,634)$ was set to 1 ft.

• Since our model is designed exclusively for radial networks, certain links were deleted from the standard IEEE networks so that the modified network is radial. For example, the links $(151,300)$ and $(54,94)$ were deleted from the standard IEEE-123 node network (see Figure A3 in the Appendix). Additionally, for the 123-node network, the switches between nodes 151 and 300 and between nodes 54 and 94 were deemed to be in an open state.

Our simulations were conducted using the Gurobi [31] solver with the following parameters:

• the rewards for all jobs are equal to $1$, i.e., $r_{i}=1,\,\forall i$,

• the penalties for all unfinished jobs are equal to $0$, i.e., $p_{i}=0,\,\forall i$ (therefore, we solve a reward maximization problem (eqn. 1a) instead of profit maximization (eqn. 5)),

• job repair times were drawn from a Weibull distribution with varying means, subject to a lower bound of $30$ minutes (consequently, all repair times are 30 minutes or longer),

• travel distance between job sites $i$ and $j$ was set equal to the shortest path distance between $i$ and $j$, which in turn was derived from IEEE specifications, subject to the modification discussed above (in other words, we assume that the road infrastructure network is identical to the power network),

• travel speed was adjusted to yield varying travel time matrices (the travel time between jobs $i$ and $j$ is equal to the shortest path distance between $i$ and $j$, divided by the travel speed) matrices,

• all nodes other than the source are damaged and require repair,

• the time budgets for all crew are identical, and

• the ‘home nodes’ (node 0) are node 650, 800, and 150 (see Figures A1, A2 and A3 in the Appendix) for the 13, 34, and 123-node networks respectively.

Simulations on the 13-node network were carried out on an M1Pro laptop with 16 GB RAM. Simulations on the 34 and 123 node networks were conducted on University of Washington’s Hyak supercomputing system [32]. For each simulation on the supercomputing system, we allocated 20 tasks per node and 80 GB RAM (total). Since model complexity depends on the network size (as discussed in Section V), networks of varying sizes require different computational times. All problem instances on the IEEE 13-node network were solved to optimality. However, for the 34 and 123-node networks, we specified a termination criterion of 13% duality gap or 12 hours solver time, whichever occurs first. The mean and worst-case gap (defined as UB - LB, where UB denotes the best feasible upper bound and LB denotes the best feasible lower bound) we observed after program termination, over all 34-node trials and all 123-node trials, are: (i) 34-node: mean 1.25, worst 3, (ii) 123-node: mean 2.05, worst 7.

In the subsequent subsections, we evaluate our model on three performance metrics: (i) normalized aggregate reward per crew, defined as the aggregate reward (i.e, the value of the objective function upon program termination) divided by the product of number of nodes ($n-1$) and number of crew ($m$), as a function of mean repair time and time budget per crew, (ii) aggregate reward as a function of number of crew and repair to travel time ratio, and (iii) normalized unused work time, defined as the ratio of total unused work time over all crew divided by the total work time over all crew, as a function of repair to travel time ratio and number of crew. Note that since fractional job completions are not allowed in our model (rewards earned on job completion, provided electrical continuity constraints are met, are therefore all-or-nothing), the total work time (i.e., the total time spent on repair and travel) for any crew need not be exactly equal to the time budget allocated to that crew. We refer to this residual time as the unused work time for that crew. There is another reason why the unused work time need not be zero. Depending on the number of crews deployed and the time budget for each crew, it may be possible for one or more crews to complete their repair schedule strictly before their respective time budgets. In this case too, the unused work time will be non-zero.

VI-A Discussion on normalized aggregate reward per crew vs. mean repair time and time budget per crew

Figure 8 shows the normalized aggregate reward (NAR) per crew as a function of mean repair time and time budget per crew for all three network sizes. Each data point in Figure 8 is obtained by averaging over multiple trials: $100$ for the IEEE-13 node network and $10$ for the 34 and 123-node networks. Since deploying an identical number of crews for all three network sizes is unrealistic, we elected to plot the NAR/crew instead of the NAR. Simulations for the 13-node network were conducted using 4 crews ($m=4$), while those for the 34 and 123-node networks were conducted using $m=8$. Observe that the maximum NAR per crew that can be achieved for the 13-node network with $m=4$ is $1/4=0.25$, while for the 34 and 123-node networks, it is equal to $1/8=0.125$ when $m=8$. Additionally, the travel speed was set such that the mean repair to travel time ratio is 3:1 for all network sizes when the mean repair time is 47.725 minutes. This required using a speed of 141 ft/min ($\approx 1.6$ mph) for the 13-node network, 4740 ft/min ($\approx 54$ mph) for the 34-node network, and 225 ft/min ($\approx 2.6$ mph) for the 123-node network. Although the speeds for 13 and 123-node networks are unrealistic, using these values in conjunction with IEEE-specified edge lengths resulted in realistic travel times, in the range of $[0.0071,36.17]$ minutes for the 13-node network, $[0.00021,40.825]$ minutes for the 34-node network, and $[0.025,216.325]$ minutes for the 123-node network. Keeping the travel speeds (or equivalently, travel times) fixed at the above values, we then varied the mean repair time and the time budget/crew to generate the figures shown in Figure 8.

Figure 8(a) shows the NAR/crew for the 13-node network. When the time budget/crew is very high (e.g., 240, 300, or 360 minutes) relative to repair and travel times, we observe that the NAR/crew achieves its highest possible value (0.25) when the mean repair time is swept between $47-66$ minutes (approximately). Conversely, for very low time budgets (60 minutes), the NAR/crew remains constant at approximately $1/12\approx 0.0833$, which corresponds to 4 completed jobs out of possible 12. For moderate time budgets, we expect the NAR/crew to be monotonically decreasing as the mean repair time is increased (keeping the travel times fixed), and this is indeed what we observe when the time budget per crew is 120 or 180 minutes. Similar inferences can be made for the 34 and 123-node networks (see Figures 8(b) and 8(c)). However, since instances of these networks weren’t solved to optimality, we observe local deviations from the expected non-decreasing behavior (e.g., in the neighborhood of mean repair time $=52$ minutes when the time budget is 240 or 300 minutes in the case of the 34-node network). Additionally, when the time budget is very small (60 minutes) and the mean repair time exceeds 56 minutes (approximately), no jobs could be completed by any crew, ensuring electrical continuity constraints, which resulted in a NAR/crew of zero. This explains why the $x$-axis is limited to approximately 52 minutes when the time budget/crew is 60 minutes.

VI-B Discussion on aggregate reward vs. number of crew and mean repair to travel time ratio

Figure 9 shows the aggregate reward as a function of number of crew and repair to travel time ratio for all three network sizes. For convenience, the $y$-axes are dual scaled to show the aggregate reward as well as the normalized aggregate reward (NAR). Note that per crew normalization is not performed (unlike Figure 8) since our intent is to understand the marginal benefit due to deployment of an additional crew. For this analysis, we fixed the repair time vectors (approximately 88 minutes for 13-node network and 79 minutes for 34 and 123-node networks) for each of the three network sizes, and varied the travel speed to obtain different travel time matrices/distributions (see Figures A4, A5, and A6 in Appendix for representative distributions). This resulted in different mean repair time to mean travel time ratios. Since the travel speed and the repair time vectors are assumed to be deterministic, each data point in Figure 9 is derived from a single trial.

Intuitively, for a fixed repair to travel time ratio, we expect the aggregate reward to be non-decreasing, which is confirmed in Figure 9. Additionally, for a fixed repair time vector, a larger repair to travel time ratio (e.g., 12:1) implies smaller commute times which allows for more jobs to be completed (higher aggregate reward), while a smaller repair to travel time ratio (e.g., 1:12) implies larger commute time which results in fewer job completions (smaller aggregate reward). A visual examination of the color coded sequence of plots in Figure 9 confirms this behavior. We also observe that the trend for the 13-node network shows a saturation behavior as $m$ increases, which is unlike what we see for the other two network sizes. In general, this saturation phenomenon can occur for any network size if all jobs can be completed by a certain number of crew, resulting in zero marginal aggregate reward (defined as $\text{AR}(m+1)-\text{AR}(m)$, where $\text{AR}(m)$ denotes the aggregate reward from $m$ crews, or alternately, the slope of the aggregate reward curve between successive values of $m$) due to deployment of an additional crew. See Figure A7 in the Appendix for plots of marginal aggregate reward corresponding to Figure 9.

Since an analysis of marginal aggregate reward (MAR) w.r.t number of crew can provide useful insights for personnel planning purposes, an important question we wanted to address in this section is, ‘are there critical values of $m$ at which the MAR increases sharply?’ While we did not observe such behavior over $m\in[5,10]$, an examination of Figure A7(c) in the Appendix (for 123-node network) does provide some interesting insights. Recall that the repair times are held constant. When the mean repair to travel ratio is relatively small (e.g., 1:12, 1:9, 1:6, and 1:3), we observe that the MAR hovers between 1 and 2. For example, when the ratio is 1:12, (i) MAR $=1$ when the number of crews increases from 5 to 6, (ii) MAR $=2$ when the number of crews increases from 6 to 7, (iii) MAR $=1$ when the number of crews increases from 7 to 8, and (iv) MAR $=2$ when the number of crews increases from 8 to 9, or from 9 to 10. In contrast, when the mean repair to travel ratio is 1 or higher, the MAR is either 2 or 3. For example, when the ratio is 12:1, (i) MAR $=2$ when the number of crews increases from 5 to 6, (ii) MAR $=3$ when the number of crews increases from 6 to 7, (iii) MAR $=2$ when the number of crews increases from 7 to 8, and (iv) MAR $=3$ when the number of crews increases from 8 to 9, or from 9 to 10. The preceding discussion should not be construed to imply that critical values of $m$, as alluded to in the beginning of this paragraph, do not exist. In general, the network topology is an important factor which can affect the marginal aggregate reward w.r.t the number of crews deployed and the mean repair to travel time ratio. For example, in a rural environment, the distribution network, despite being radial, may have locally dense clusters of repair jobs with large travel distances/times for any inter-cluster travel. For such network topologies, it is conceivable that deployment of an additional crew deployment may result in a sharp increase in the marginal aggregate reward.

VI-C Discussion on normalized unused work time vs. repair to travel time ratio and number of crew

Figure 10 shows the normalized unused work time (NUWT, can be viewed as a measure of crew utilization) as a function of mean repair to travel time ratio and number of crew. The repair and travel time parameters for Figure 10 are identical to those used for Figure 9. As in Figure 9, each data point in Figure 10 is derived from a single trial.

Figure 10(c) shows the normalized unused work time for the 123-node network. For a fixed repair time vector, when the ratio of the number of available repair jobs to number of crew is high (recall that we assume that all but the source node is damaged) and the time budgets are reasonable, we can expect that most of the crews will be ’maximally utilized’, i.e., be engaged in either repair/travel for a substantial fraction of their time budgets, particularly when travel times are small (or alternately, large mean repair to travel time ratios). This should result in small unused work time. For large travel times, it may not be possible for crews to make multiple commutes, possibly resulting in a larger unused work time. Figure 10(c) confirms this behavior. We also observe that the unity repair to travel ratio is an elbow point of the family of NUWT curves.

The 34-node network (see Figure 10(b)) behaves similarly to the 123-node network, except there is a steep transition between the 1:3 ratio and the unity ratio. Although we are not entirely sure what causes this sharp transitory behavior for the 34-node network, we posit that it is due to the relative scarcity of the number of repair jobs (33 for the 34-node network compared to 122 for the 123-node network). When travel times are large, it is plausible that the repair crews are faced with relatively few choices for selecting additional repairs without violating their time budget, leading to higher (compared to the 123-node network) unused work time when the mean repair to travel time ratio is less than unity. When the ratio is unity or higher, the travel times become comparable to or smaller than the repair times in an average sense, thereby opening up more choices for the crews to travel and complete additional repair jobs.

Finally, we observe that the 13-node network (see Figure 10(a)) behaves somewhat differently compared to the 34 and 123-node networks. With only 12 repair jobs, deploying 10 crews will naturally lead to one or two crews completing multiple jobs while the rest of the crews complete only one job (no travel time is involved for these crews since we assume that time accounting starts when a crew is stationed at the first repair site), leading to a high unused work time. This is manifested in the relatively flat characteristic throughout the range of repair to travel ratios when $m=10$. As the number of crews decreases, we start to see a behavior which resembles the behavior in Figures 10(b) and 10(c), with an ‘elbow effect’ appearing for fewer crews at a ratio of 3:1.

VII Conclusion

In this paper, we address the problem of optimal restoration of electrical distribution networks after a severe weather event. Since the collection of damage information is a streaming process and prone to uncertainty or incompleteness, particularly during the initial hours/days after the event, we adopt a rolling horizon restoration framework. Incorporating both repair times and travel times, we propose a mixed integer linear programming model for solving the problem optimally on a restoration window basis. A key feature of this model is the inclusion of electrical connectivity constraints which ensures that all nodes scheduled for repair in the upcoming restoration window can be energized at the end of the window. We show that when the power network is of star topology, the optimal recovery problem reduces to an instance of the cost constrained reward maximizing $m$TSP on doubly weighted graphs, a novel variant of the traveling salesman problem, which, to the best of our knowledge has not been studied before. We validate our model on multiple IEEE test networks and evaluate the objective function w.r.t number of crews deployed, time budgets, mean repair time, and mean travel time.

(a)

(b)

(c)

Figure 8: Plots of normalized aggregate reward (NAR) per crew as a function of mean repair time and time budget per crew. (a) 13-node and $m=4$, where $m$ denotes the number of crew, (b) 34-node and $m=8$, and (c) 123-node and $m=8$.

(a)

(b)

(c)

Figure 9: Plots of aggregate reward (AR) and normalized aggregate reward (NAR) as a function of number of crew and mean repair to travel time ratio. (a) 13-node (b) 34-node, and (c) 123-node. The time budget for each crew is 180 minutes.

(a)

(b)

(c)

Figure 10: Plots of normalized unused work time (NUWT) as a function of mean repair to travel time ratio and number of crew $m$. (a) 13-node (b) 34-node, and (c) 123-node. The time budget for each crew is 180 minutes.

Figure A1: IEEE 13-node network.

Figure A2: IEEE 34-node network.

Figure A3: IEEE 123-node network.

Figure A4: Histogram of travel times for IEEE-13 node network. The mean/median travel times are 26.22/29.42 minutes respectively. A speed of 85 ft/min ($\approx 1$ mph) was used to compute the travel time matrix from the IEEE specified distance matrix.

Figure A5: Histogram of travel times for IEEE-34 node network. The mean/median travel times are 24.82/24.18 minutes respectively. A speed of 2875 ft/min ($\approx 32.7$ mph) was used to compute the travel time matrix from the IEEE specified distance matrix. Approximately 20% of travel times are less than 5 minutes. This distribution is significantly different from the travel time distributions of the 13-node (Figure A4) and 123-node (Figure A6) networks, even though the mean/median travel time for all three network sizes are similar.

Figure A6: Histogram of travel times for IEEE-123 node network. The mean/median travel times are 26.26/25.20 minutes respectively. A speed of 137 ft/min ($\approx 1.6$ mph) was used to compute the travel time matrix from the IEEE specified distance matrix.

(a)

(b)

(c)

Figure A7: Plots of marginal aggregate reward (MAR) as a function of mean repair to travel time ratio and number of crew $m$. (a) 13-node (b) 34-node, and (c) 123-node. The time budget for each crew is 180 minutes.