{'fraction_non_alphanumeric': 0.035629453681710214, 'fraction_numerical': 0.012183875925667179, 'mean_word_length': 4.163924963924964, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 9, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 9, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Repair Crew Routing for Power Distribution Network Restoration

As extreme weather events become more frequent and disruptive, service restoration is increasingly important for many infrastructures, e.g., power grids and communication networks. In many studies on service restoration, the logistics issue of traveling over the road network, however, is often overlooked due to the complexity of considering multiple networks simultaneously, resulting in prolonged disruption time. In this work, we address such a problem arising in power systems, where technical crew and utility trucks travel to a number of sites to repair damaged equipment, with the goal of minimizing the total service disruption time within the service region. We call this problem the Power Restoration Traveling Repairman Problem (PRTRP). What makes it significantly more challenging than a typical routing problem is that the service disruption time in a location depends on the interaction of the routing sequence with both networks, i.e., the road network and the power grid. To solve the problem, we develop an exact method based on bi-directional dynamic programming. We then improve the method by reducing the search space with solution upper and lower bounds, and threshold rules derived from the precedence relations in the power grid. We also propose efficient heuristic variants of the method. We present computational results and compare our method with benchmark heuristics.


Introduction
In recent years, extreme weather events and natural hazards, such as hurricanes, wildfires, and storms, have been causing a growing number of significant disruptions to various critical infrastructures. With severe physical damages, service restoration could take days, weeks, or even months, causing an enormous economic loss. Recent examples include weeks-long power outage in Florida in 2017 by Hurricane Irma [20] and months-long power outage in Puerto Rico in 2017 caused by restore the disrupted parts on the network to ensure the power service is restored as quickly as possible. Our study mainly focuses on response strategies.
The literature on disaster recovery for power systems is relatively young, with most of the studies published after 2010 [8,28]. These studies can be divided into two groups based on the part of the power network considered: distribution network restoration and transmission network restoration. Although the problems in both networks are similar due to the underlying precedence relations, the high voltage on the transmission network requires a precise model of the power flow problem. [5] is one of the early works on addressing the routing issue of repair crew on distribution networks. They propose a two-stage model, where the first stage allocates damaged locations to clusters and the second stage determines the route in each cluster. The number of disruptions in their test instances goes up to 17. [29] follows a different approach. In case of a disruption, they first solve the reconfiguration problem to ensure all nodes are supplied, and then solve the restoration problem to determine the order of repairs without considering the routing decisions. [28] emphasizes the difference between disruption restoration and blackouts restoration. Blackouts generally occur due to an instability in the power system, and the service restoration involves determining the sequence of switching actions to re-energize the demand nodes. On the other hand, disruptions caused by external factors such as natural disasters need to be attended by a repair crew. The authors study the distribution network repair and restoration problem, where travel time is ignored and a repair schedule is found with heuristics based on linear programming relaxation. [8] emphasizes the routing aspect of the distribution network restoration. They consider two groups of crews with different capabilities and propose an integer programming model to determine routes for the repair crew. To our knowledge, [4] is the only study considering the uncertainty in distribution network restoration. They assume that disruption locations are known, but travel times and demand are uncertain. In their two-stage model, they first determine the routes for the repair crew and then the reconfiguration scheme to re-energize the demand nodes.
In transmission network restoration, [32] and [26] focus on the last-mile disaster recovery for power restoration. They consider a disrupted transmission network to determine the schedules and routes for the repair crew to restore the power as quickly as possible. They define this problem as Power Restoration Vehicle Routing Problem. The authors solve the problem in three steps. In the first step, they determine the minimum set of disrupted items to repair. In the second step, they determine the order of restoration without considering routing. In the third step, they determine the routes for the repair crew considering the partial ordering imposed by the first step. [9] studies the same problem with a focus on comparing the performances of two different power flow models, namely direct current and linear approximation of alternating current (LPAC). [31] considers the transmission network restoration problem using the LPAC model, with a similar formulation and solution methodology to that in [9]. [10] is the first study on the last-mile restoration problem considering the inter-dependency between the power and gas networks. They solve the problem in three steps and use a randomized adaptive decomposition method. They consider disruption scenarios with 10 to 120 damaged locations. [2] determines the optimal repair schedule, unit commitment, and system reconfiguration to restore disrupted power networks. They develop a MIP formulation and use Benders decomposition. Their test instances are based on IEEE-118 instance with seven damaged nodes. [11] proposes a model to determine the schedule to repair the disrupted transmission network. They ignore the routing decision of the repair crew, and test their model on small instances with six disrupted nodes.
Our underlying problem is a routing problem. Considering the cumulative nature of our objective function, the Cumulative Capacitated Vehicle Routing Problem (CCVRP) and the TRP are the two most closely related problems. [17] provides the first exact solution method for the CCVRP based on a branch-and-cut-and-price algorithm. [24], [19], and [25] present metaheuristics for the CCVRP based on large neighborhood search, a memetic algorithm, and multi-start local search, respectively. [33] proposes an exact solution method based on dynamic programming and branch and bound. They find solutions for instances with up to 23 nodes. For the TRP, [16] proposes an exact solution approach for the multiple TRP based on a branch-and-price-and-cut method. They provide solutions for instances up to 50 nodes. [3] develops a 7.18-approximation for the TRP.
One study that is most closely related to ours is [18], in which the joint routing and scheduling problem is considered for restoring an infrastructure network. They propose two MIP models, one for binary restoration and one for proportional restoration, and develop heuristics to solve the models. [30] considers a disrupted road network and develops a solution methodology based on approximate dynamic programming. While this paper has similarities to ours in terms of solution methodology, we consider a power network in parallel to a road network. [12] also emphasizes the logistics aspects of the response strategies, but they consider logistics as a part of a higher level planning with three geographical areas.
The current literature on crew routing for repairing infrastructure networks is still very limited in terms of understanding the impact of precedence relations. [27] addresses this gap in the literature in a recent review. In this study, we develop a new solution approach for the single crew routing problem to minimize the negative impact of infrastructure network disruptions. In the following section, we describe our problem in more detail and present the formulation.

Problem Definition, Formulation, and Complexity
Before presenting the problem, we first provide background information on power distribution networks as an example of infrastructure networks with precedence relations. A typical power grid consists of power stations, substations, the transmission network, and the distribution network. The distribution network is considered to be the last mile of the power distribution. It delivers the low-voltage power to final customers. There are different configurations of power network topology such as radial, multi-ring, and mesh structure, offering different resilience, costs, and complexities to manage. The radial network configuration, where there are no close loops, is the most common configuration for distribution networks [21]. While it is the simplest and cheapest topology, if a node on the network is disrupted, all downstream nodes lose power service. Figure 1 shows the IEEE 13 bus feeder instance as an example of a power distribution network. On this network, node 50 is the source. If there is a disruption on node 84 only, in addition to 84, nodes 52 and 911 lose power as well. We study how to efficiently restore a power distribution network where faults occur in some locations and result in outage for customers over the network. In this paper, we assume that faults happen only at the vertices on the power distribution network and their locations are known. The utility company sends out a repair crew to visit these fault locations sequentially and make repairs. The goal is to determine a sequence of fault locations to visit for the repair crew to minimize the total service disruption time for all customers.

Problem definition
There are two underlying networks in our problem: the road network on which the repair crew travels, and the power distribution network through which electricity is sent to customers. In our network model, we use a vertex to model both the intersection of two road segments in the road network and a fault location in the power network. We first introduce some notation for the two networks.
Let vertex 0 denote the depot where the repair crew is located before the repair starts. Let V c = {1, 2 . . . , n} be the set of all fault locations (power substation or a customer location). All locations in V c receive power from the substation through distribution lines in the power distribution network, which is a directed acyclic graph. We consider the case in which there is only one substation in the network. For i, j ∈ V c , an arc (i, j) denotes that vertex i is the immediate predecessor of vertex j on the unique path from the power source to vertex j in the power distribution network. The tree T = (V c , A c ) characterizes the precedence relations among all fault locations in the power network where the arc set A c is the collection of precedence relations on the power distribution network. On the other hand, we model the road network as a complete directed graph G = (V, A) where V = V c ∪ {0} and each arc (i, j) ∈ A represents a shortest path from location i ∈ V to location j ∈ V . The travel time over arc ( The main complexity of our problem compared to a typical routing problem comes from how the service disruption time at a fault location is computed. Let vertex 1 ∈ V c denote the substation (the unique power source in the power network). For each vertex i ∈ V c , there exists a unique path P 1i in T from the substation to vertex i. The power service at vertex i can only be recovered when all faults on the path P 1i are removed. Therefore, the service disruption time at vertex i is the duration between the time the repair crew leaves the depot (vertex 0) and the time all faults on path P 1i (including the fault at vertex i) are removed. Let r i denote the time it takes until the recovery of power service at vertex i ∈ V c . Note that r i depends on both the topology of the power network T as well as the sequence in which fault locations are visited in the road network G. Our goal is to determine a Hamiltonian cycle over vertices V starting at vertex 0 to minimize the total service disruption time i∈Vc r i .

Formulation
Before introducing the formulation of PRTRP, we assume that the repair duration p i at each vertex i ∈ V c is 0. To see why this assumption is without loss of generality, we can transform the original problem into an equivalent PRTRP with zero repair duration at each vertex, by changing two inputs: set the travel time d ji in the new problem to be d ji + p i for each incoming arc (j, i) ∈ A in the road network and p i = 0. This transformation does not change the objective value of any repair sequence and therefore creates an equivalent PRTRP with no repair duration. We also introduce some notation to assist the calculation of the service disruption time at a vertex. Given two vertices i and j in V c , let i j denote that vertex i is a predecessor of vertex j in the power network T , i.e., vertex j is reachable from vertex i through a directed path on T . If vertex i does not receive power service, neither does vertex j. Let i j denote that either i j or i = j. Given integers i, j with i ≤ j, we use [i : j] to denote the set of integers {i, i + 1, . . . , j}.
We now introduce the formulation of the problem. Let the binary decision variables x ij denote whether the repair crew travels from vertex i to vertex j for each arc (i, j) ∈ A, the continuous decision variables t i denote the arrival time at vertex i for each i ∈ V , and the continuous decision variables r i denote the service disruption time of vertex i for each i ∈ V c . Then the PRTRP can be formulated as a MIP as follows.
j∈V \{i} The objective function in (1) calculates the total service disruption time at all fault locations. Constraints (2) enforce that each faulty vertex on the distribution network should be visited exactly once. Constraints (3) ensure that vertex j will be visited after vertex i if the arc (i, j) is used. Constraints (4) are the unique constraints compared to other routing problems. They ensure that vertex j will only receive the power service after all its predecessors and itself in the power grid are visited and repaired. We need to set the value of M in (3) large enough for the formulation to be valid, for example setting M = (k,l)∈A d kl . It is possible to choose smaller values for M . However, it is well known that this type of big-M formulations provide weak linear programming relaxations and take a long time to solve for large-size instances.
Complexity: It is not surprising that the described problem is NP-hard.
Proof. We will show a polynomial-time reduction from the NP-hard traveling repairman problem to the PRTRP. Consider an instance of the traveling repairman problem over a graph N = (V , A ), where V = {1, . . . , n} and the travel time over arc (i, j) ∈ A is d ij . Construct an instance of PRTRP in the following way. For the road network, set V c = V , put the depot 0 at the same location as vertex 1, and let the travel time from vertex i to vertex j be d ij for (i, j) ∈ V c . For the power network, make vertex 1 the power source and each vertex i ∈ V c \{1} the immediate successor of vertex 1; the power network is essentially a star with n − 1 leaves. Finding an optimal solution to the traveling repairman problem over graph N is equivalent to finding an optimal solution of the PRTRP over the constructed road and power networks, and the transformation takes polynomial time.

Solution Methodology
To solve the PRTRP, we develop an algorithm based on bi-directional dynamic programming (DP). Our algorithm performs backward and forward DP at the same time and uses the combined information to efficiently reduce the search space. It is also strengthened by pre-processing steps that are developed based on the precedence relations. Before presenting the details of our algorithm, we first introduce some definitions and new notation.
We call a path P a return path if it starts at some vertex i ∈ V c and ends at vertex 0.
Definition 2. We call a path P an outgoing path if it starts at vertex 0 and ends at some vertex Given a path P , let V (P ) denote the set of vertices on path P . Given two paths P 1 and P 2 , let P 1 ⊕ P 2 denote a path obtained by concatenating P 2 after P 1 . For example, if P 1 = (0, 3, 5, 4) and P 2 = (1, 2), then P 1 ⊕ P 2 = (0, 3, 5, 4, 1, 2). Assume that P is a return path starting at vertex i ∈ V c . Let v B (P ) denote the service disruption time for all customers in V c between the time the vehicle arrives at vertex i and the time the vehicle returns to vertex 0 following path P . Given an outgoing path P ending at vertex i ∈ V c , let u F (P ) denote the service disruption time for all customers between the time the vehicle leaves the depot and the time vehicle arrives at vertex i following path P . We introduce some key concepts and results that will be used in DP.
Definition 3. Two return paths P 1 and P 2 have the same configuration if P 1 and P 2 start at the same vertex and V (P 1 ) = V (P 2 ). Similarly, two outgoing paths P 1 and P 2 have the same configuration if P 1 and P 2 end at the same vertex and V (P 1 ) = V (P 2 ).
Proposition 2. If two return paths P 1 and P 2 have the same configuration and v B (P 1 ) < v B (P 2 ), then no optimal Hamiltonian cycle of the PRTRP ends with P 2 .
Proof. Suppose that there exists an optimal Hamiltonian cycle P that ends with P 2 . Assume that P = P 3 ⊕ P 2 . Then P = P 3 ⊕ P 1 is a Hamiltonian cycle leading to a smaller total service disruption time than that of P , which contradicts the optimality of path P .
Proposition 3. If two outgoing paths P 1 and P 2 have the same configuration and u F (P 1 ) < u F (P 2 ), then no optimal Hamiltonian cycle of the PRTRP starts with P 2 .
Proof. Suppose that there exists an optimal Hamiltonian cycle P that starts with P 2 . Assume that P = P 2 ⊕ P 3 . Then P = P 1 ⊕ P 3 is a Hamiltonian cycle leading to a smaller total service disruption time than that of P , which contradicts the optimality of P .
In the BiDP algorithm, we construct outgoing paths and return paths simultaneously, and use additional bounding techniques to reduce the search space.

Backward dynamic programming
We denote the state vector of the backward DP by (k, Y ), where vertex k is the current vertex, i.e., the current location of the repair crew, and Y ⊆ V c is the set of vertices that are not visited (repaired) yet. Note that k / ∈ Y . Let v(k, Y ) denote the minimum total service disruption time of all customers between the time the crew leaves vertex k and the time the crew finishes removing all faults. Then we have the following Bellman equations.
where W : 2 V → R >0 and W (S) computes the number of vertices whose power services are disrupted when the set of vertices in S are repaired. Equation (8)

Forward dynamic programming
We denote the state vector of the forward DP by (X, k), where X ⊆ V is the set of vertices that have been visited so far and vertex k is the current location of the repair crew. Note that k ∈ X.
Let u(X, k) denote the minimum total service disruption time of all customers between the vehicle leaves the depot, visits all customers in X, and arrives at vertex k. We can write the Bellman equations for the forward DP as follows.
Let u * be the minimum service disruption time among all Hamiltonian cycles. Then u * = min i∈Vc u(V, i).

Pre-processing and the maximum position of a vertex in optimal solutions
In the pre-processing step, we determine the maximum position of a vertex in any optimal Hamiltonian cycle and use that to reduce the search space during BiDP.
We first introduce some notation that will be used in such bounds. Let The proof of Proposition 4 is given in Appendix B. We can use the lower bound associated with vertex i in Proposition 4 to bound the maximum position of the vertex in an optimal solution.
Proposition 5. Let U be an upper bound for the optimal objective value of the PRTRP. Given a vertex i ∈ V c and some position k > n − |S i |, if the lower bound L k i defined in (14) satisfies that L k i > U , then any optimal Hamiltonian cycle for the PRTRP will not visit vertex i at the k th position or later.
Proof. For the lower bound L k i defined in (14), if L k i > U , then it is not optimal to visit vertex i in the k th position of a Hamiltonian cycle. To prove the proposition, it suffices to show that for a given vertex i, L k i is monotonically non-decreasing in k. According to (14) , Therefore, if it is suboptimal to visit vertex i on the k th position, it is also suboptimal to visit it on later positions.
The monotonicity of the lower bound L k i in Equation (14) is due to the condition that k > n − |S i |. Since |S i |> n − k, some of the vertices visited before the k th vertex are successors of vertex i. These vertices cannot receive power before i is recovered, so delaying the repair of i delays their service recovery, thus increasing the associated lower bound. Proposition 5 implies that for vertex i ∈ V c , if we find a position k that satisfies L k i > U , then any solution where vertex i is visited on k or later could be removed from the search space of the dynamic programming.

Lower bounds of solutions extended from a given path
We derive lower bounds for the objective value of Hamiltonian cycles that can be extended from an outgoing or return path. These lower bounds will be used to reduce the search space during the iterations of DP. Recall that s p denotes the length of the p th shortest arc in the road network G. Proposition 6. Given an outgoing path P that contains k vertices in V c , let w P be the number of vertices without power service after the crew travels through P . Then the total service disruption time of any Hamiltonian cycle extended from P is lower bounded by Proposition 7. Given a return path P that contains k vertices in V c , the total service disruption time of any Hamiltonian cycle extended from P is lower bounded by The proofs of the propositions above are in Appendices C and D respectively.

Upper bounds from heuristics
To create initial solutions and upper bounds for optimal solutions, we use two greedy heuristics. The heuristics use different criteria based on the distance and the priority of vertices on the power distribution network. The heuristics work as follow.
• Greedy in distance (GiD): Always go to the nearest unvisited location.
• Greedy in priority-weighted distance (GiPD): Given the current location i, go the the unvisited vertex j with the smallest d ij /|S j | value.

Bidirectional dynamic programming
We present our bi-directional DP-based solution method in Algorithm 1. It starts with calculating initial feasible solutions and upper bounds, and the maximum position for each vertex in a Hamiltonian cycle. Then we perform the backward and forward DP simultaneously. For each outgoing path or return path, we calculate a lower bound following the ideas in Propositions 6 and 7 to determine whether to accept this path or not. We also update the upper bound periodically by completing the outgoing paths greedily. for each return path P in Γ B l do 5: Generate a path P = {i} ⊕ P .

8:
Compute a lower bound LB(P ) for the return path P according to Proposition 7.

9:
if U ≥ LB(P ) then 10: for each outgoing path P in Γ F l do 15: for each vertex i ∈ V c \ V (P ) such that l ≤ β i do

18:
Compute a lower bound LB(P ) for the outgoing path P according to Proposition 6.

19:
if U ≥ LB(P ) then 20: if P ⊕ P is a Hamiltonian cycle and u F (P ) + v B (P ) < U then 29: U ← u F (P ) + v B (P ) and store P ⊕ P . The subroutines for BackwardUpdate(S, P ) and ForwardUpdate(S, P ) work as follows.

Algorithm 2 The subroutine BackwardUpdate(S, P )
Input: A set S of return paths and a return path P . Output: An updated set S. 1: if there is no return path P ∈ S with the same configuration as P then 2: Add P to S 3: else 4: Let path P ∈ S be the return path with the same configuration as P 5: if v B (P ) < v B (P ) then 6: Delete P from S and add P into S 7: end if 8: end if 9: Return S.

Algorithm 3 The subroutine ForwardUpdate(S, P )
Input: A set S of outgoing paths and an outgoing path P . Output: An updated set S. 1: if there is no outgoing path P ∈ S with the same configuration as P then 2: Add P to S 3: else 4: Let path P ∈ S be the outgoing path with the same configuration as P

5:
if u F (P ) < u F (P ) then 6: Delete P from S and add P into S

Computational Experiments
In this section, we present computational results of our solution method and compare it with benchmarks. The test instances are created from the Electric Power Research Institute (EPRI) instances, which are based on a U.S. power distribution network [14]. We use the ckt5 instance as shown in Figure 2. The instance ckt5 represents the power distribution network, and we assume the road network is a complete graph. To create test instances over a power network of a smaller size, we take subtrees from the original power network. In our naming convention, ckt5 x represents the subtree of the original power network rooted at node x, i.e., node x becomes the power source of the new instance. We assume that a unit distance can be traveled in a second, and this can be adjusted for the speed of travel.
We compare the optimal results obtained by our BiDP with the benchmark heuristics GiD and GiPD introduced in Section 4.4.2. We omit the results obtained by the commercial optimization solver CPLEX since only very small-sized instances can be solved to optimality within the given time limit. All algorithms are implemented in the C programming language, and the experiments are performed on a computer that has a 2.8 GHz Intel i7 CPU with 4 cores and 16 GB of RAM and runs the macOS operating system. In the computational results, we report the resulting total service disruption time (z) and the computation time (t) in seconds for each algorithm. We omit the computation time for the benchmark heuristics since they are very quick to find a solution (in less than a second).

Exact solutions by BiDP
We present solutions found by our BiDP algorithm and benchmark heuristics in Table 1. Our BiDP algorithm finds optimal solutions for instances with fewer than 20 vertices in less than 165 seconds. GiD and GiPD deviate from the optimal objective by 27% and 22% on average, while the maximum deviation goes up to 50% and 62% respectively. The BiDP shows a clear advantage over the two heuristics in terms of solution quality. This could be weighed against the computation time in practice. We also observe that the computation times for instances with the same number of vertices can be significantly different. Consider instances ckt5 361 and ckt5 805. They both have 21 vertices. The optimal solution for ckt5 361 is found in 199.69 seconds while it takes 3459.41 seconds for ckt5 805. Similarly, ckt5 144 and ckt5 820 both have 18 vertices, and computation times are 18.18 and 164.15 seconds respectively. The main factor behind these differences is the structure of the power distribution network in the instances, as shown in Figures 3 and 4. As more vertices branch out (i.e., having multiple immediate successors) and the vertices that branch out are closer to the power source, the computation time of our algorithm significantly increases. This is mainly due to the fact that the upper bounds in instances with fewer vertices branching out are tighter, which keeps the size of the search tree of BiDP small.

Heuristic Implementations
Motivated by the previous observation, we also implement our BiDP algorithm heuristically. In the first setting, we modify the use of the upper bounds. Recall that we initialize the algorithm with an upper bound that is generated by the greedy heuristics, and update the upper bound periodically as we progress with the bi-directional search. The quality of the upper bound improves over time. Therefore, we adjust the upper bound with coefficients θ and δ. Accordingly, we allow the outgoing and return paths into our search tree if the corresponding lower bound is less than or equal to (θ + lδ)U , where l is the number of vertices that have been added to the outgoing or return path. Therefore, we start with an aggressive use of the upper bound and relax it as we progress, creating a dynamic heuristic upper bound. For example, for θ = 0.80 and δ = 0.01, the effective upper bound when the first vertex is added to the outgoing or return path is 0.80U , and it is 0.87U when adding the 8 th vertex. We report the computational results for (θ = 0.80, δ = 0.01), (θ = 0.90, δ = 0.01), and (θ = 1.00, δ = 0.00) in Table 2. Note that the last one corresponds to the exact solution method BiDP. We also report the average, minimum and maximum optimality gap (i.e., percentage deviation from the optimal solution), and the average, minimum and maximum percentage reduction in computation time compared to t(opt) for each (θ, δ) pair. We highlight the cases where we find optimal solutions with the heuristic upper bounds. In the computational results, the average optimality gap is 0.00% for each (θ, δ) pair. Using (θ = 0.80, δ = 0.01), we find optimal solutions in 19 instances out of 21. In the remaining instances, the deviations from optimal solutions are less than or equal to 1%. The average percentage reduction in computation time is 68%, and the maximum reduction is 90%. Using (θ = 0.90, δ = 0.01), we find the optimal solution in all but one instance with 11% average reduction in computation time. In the next heuristic implementation, in addition to the modified upper bounds, we also use a modified maximum position for the power source. After finding greedy solutions with GiP and GiPD, we check the positions of the power source in both solutions, and use the greater position to determine the maximum position for the power source during BiDP. Table 3 presents the corresponding solutions. We repeat the optimal solutions in the table for completeness. Optimal solutions found by the heuristic approaches are highlighted. Using this heuristic implementation, we can still find optimal solutions in almost all instances with 0.00% average optimality gap. The results corresponding to (θ = 1.00, δ = 0.00) show the impact of the heuristic maximum position only. We observe that by just using the heuristic maximum position, computation times to find optimal solutions decrease by 22% on average and can be reduced up to 46%. When heuristic upper bounds for (θ = 0.80, δ = 0.01) are coupled with heuristic maximum position, the computation times to find optimal solutions decrease by 72% on average and can be reduced up to 93%.

Conclusion
In this paper, we study a routing problem arising in power service restoration, where the repair crew travels to a number of locations to repair power equipment with the goal of minimizing the total power service disruption time. The service disruption time of a location depends on both the travel time on the road network and the topology of the power network, which brings a new layer of complexity to this routing problem. We call this problem the PRTRP.
The PRTRP is computationally challenging to solve even for small-sized instances. We propose a bi-directional dynamic programming (BiDP) method to find optimal solutions. We develop several techniques to further reduce the search space of the dynamic programming, including precomputed bounds of the optimal position of a vertex, dynamic lower bounds of optimal objectives of extending a path to Hamiltonian cycle, and upper bounds by two primal heuristics. In addition, we propose heuristic implementations of our BiDP based on dynamic heuristic upper bounds and heuristic bounds for optimal vertex positions. We test our methods on instances created from EPRI instances. Our BiDP method is able to find optimal solutions of instances that cannot be solved by off-the-shelf optimization software. Our heuristic implementations further reduce the computation time significantly while finding optimal solutions for most instances.
There are several directions to extend our model and methods for future research. Firstly, more computationally efficient algorithms will need to be investigated to support quick decisions under emergency response. Secondly, it's worth considering routing multiple repair crew as power utility companies usually dispatch a fleet of utility trucks during restoration. An added complexity is that the total service disruption time over a route could depend on other routes, while in traditional vehicle routing problems the total travel time or distance of a route only depends on the route itself. Thirdly, service restoration often faces uncertainty in practice, e.g., equipment status and road availability. Taking these uncertainties into account in the model will help generate more robust restoration solutions. Lastly, the problem structure studied in this paper is common in service restoration of other networks, such as communication network or gas network. It would be interesting to see our results to be applied to other infrastructure networks.
denote the total service disruption time if the repair crew follows H. Let c p (H) be the number of vertices with no power right before the repair crew arrives at vertex j p following H, for p ∈ [1 : n].
Then we have the following result.
Equation (16) simply states that the total disruption time following H is to sum up the total service disruption time between visiting two consecutive vertices j p−1 and j p on H, which is the product of the number of vertices with no power service c p and the travel time d j p−1 jp .
Let C(H) = (c 1 (H), c 2 (H), . . . , c n (H)) and D(H) = (d 0j 1 , d j 1 j 2 , . . . , d j n−1 jn ). We now introduce the lower-bounding lemma. where the first equality follows from Equation (16), the first inequality follows from the Rearrangement Inequality (15), and the second equality follows from the fact that the components of C are already sorted in a non-increasing order.

B Proof of Proposition 4
Proof of Proposition 4. According to Lemma 1, for any Hamiltonian cycle H, we have T (H) ≥ C(H) · g(D(H)).
We will find vectorsC andD such that C(H) ≥C and g(D(H)) ≥D respectively.
We first construct a lower bound of the vector C(H). Consider any Hamiltonian cycle where vertex i is visited at position k with k > n − |S i |, i.e., (j 0 = 0, j 1 , . . . , j k−1 , i, j k+1 , . . . , j n ). Since k > n − |S i |, the positions on this Hamiltonian cycle can be divided into three groups as follows: Proof of Proposition 7. The proof is similar to the Proof of Proposition 6. Since path P has k vertices, we can denote any Hamiltonian cycle H that is extended from P by H = Q ⊕ P , where Q is an outgoing path that ends at the starting vertex of path P and contains n − k vertices in V c . Let Q = (j 0 = 0, j 1 , . . . , j n−k ). Let T (H) denote the total service disruption time of H. Let c p be the number of vertices with no power right before the repair crew arrives at vertex j p for p ∈ [1 : n]. Let C = (c 1 , . . . , c n−k+1 ) and D = (d j 0 j 1 , . . . , d j n−k j n−k+1 ). Now we try to find vectorsC andD such that C ≥C and g(D) ≥D respectively. Note that c p ≥ n − p + 1 since there are at least n − p + 1 vertices not visited before vertex j p is visited following H, so C ≥C = (n, n − 1, . . . , k). On the other hand, recall that s p is the length of the p th shortest arc in A. Then g(D) ≥D = (s 1 , s 2 , . . . , s n−k ). Therefore, we have T (H) ≥C ·D + v B (P ) = ns 1 + (n − 1)s 2 + . . . + ks n−k+1 + v B (P ).