Collaborative Optimization of Post-Disaster Damage Repair and Power System Operation

: After disasters, enhancing the resilience of power systems and restoring power systems rapidly can effectively reduce the economy damage and bad social impacts. Reasonable post-disaster restoration strategies are the most critical part of power system restoration work. This paper co-optimizes post-disaster damage repair and power system operation together to formulate the optimal repair route, the unit output and transmission switching plan. The power outage loss will be minimized, with possible small expense of damage repair and power system operation cost. The co-optimization model is formulated as a mixed integer second order cone program (MISOCP), while the AC-power-ﬂow model, the complex power system restoration constraints and the changing processes of component available states are synthetically considered to make the model more realistic. Lagrange relaxation (LR) decomposes the model into the damage repair routing sub problem and the power system operation sub problem, which can be solved iteratively. An acceleration strategy is used to improve the solving efﬁciency. The proposed model and algorithm are validated by the IEEE 57-bus test system and the results indicate that the proposed model can realize the enhancement of resilience and the economic restoration of post-disaster power systems.


Introduction
Power systems have become an indispensable guarantee for the prosperity and development of modern society. Traditionally, power systems are designed to be reliable enough under normal conditions and to sustain single outage contingency ("N−1" security criterion). However, on one hand, due to the current climate change, power systems are more likely to suffer disturbances from extreme natural disasters, which possibly lead to cascading outages and blackouts. In 2008, a severe ice storm hit southern China [1], resulting in the failure of 7541 lines; in 2011, more than 4 million household power supply was affected by Japan's earthquake [2]. These disasters include earthquake, strong wind, ice storm, lightning, extreme temperature, rainstorm and flood and so forth. On the other hand, increasing man-made and terrorist attacks have also greatly harmed the safe and stable operation of power systems. In 2013, the PG&E substation in Metcalf of California was shot, causing the breakdowns of 17 transformers and a serious power outage [3]. To make power systems stronger and smarter to resist disasters, the concept of resilience of power systems which has aroused extensive attention is proposed.
According to [4], power system resilience is "the ability of the power system to anticipate high-impact low-probability events, rapidly recover from these disruptive contingencies and absorb lessons for adapting its operation and structure for preventing or mitigating the impact of similar events in future." According to the evolving nature of general disasters, reference [5] divided the research of power system resilience into three parts: pre-disaster system toughness, during-disaster system resistance and post-disaster system restoration ability, as shown in Figure 1. As the core part of resilience research, the enhancement strategies of resilience can also be implemented from the three aspects. Pre-disaster improvement mainly belongs to the scope of hardening measures, while duringdisaster and post-disaster improvement are dominated by smart operational measures. Before the occurrence of disasters, it is essentially to optimize power system resilience-oriented planning [6][7][8], identify and strengthen the weak links [9,10], as well as allocate the disaster repair resources optimally [11,12]. During disasters, effective and flexible control measures [13,14] should be taken to avoid cascading failures and restrict the further expansion of failures. After disasters, how to rapidly repair the damaged components and restore curtailed load with limited resources [15][16][17] is the most urgent and significant task for power system restoration. In [18], an analysis of lifelines in case of seismic disaster under the analysis FMECA was presented. Reference [19] reviewed the challenges of maintaining the linear assets and provided a conceptual framework for the use of autonomous inspection and maintenance practices. This paper emphasizes the enhancement of resilience in the post-disaster stage.  In the previous power system post-disaster resilience improvement studies, concentrations are focused more on the restoration after large-scale blackouts, while the researches are primarily carried out from three aspects: black-start [20][21][22], network reconfiguration [23][24][25] and load restoration [26,27]. However, the interconnection of regional large power systems, reasonable protection and automatic devices, along with advanced alarm processing and fault diagnosis technology, can restrict most electrical faults and accidents in local areas [5]. Therefore, the happening probability of extreme weather events leading to large-scale blackouts is tiny. Compared with the system-wide blackout, local outages are more likely to occur. Nevertheless, the studies on the restoration of local power systems are relatively scarce and immature. Thus, it is more meaningful and practical to study the restoration of local power systems especially considering the repair of damaged components.

System Function
Currently, many studies are trying to take the damage repair work into consideration when restoring the post-disaster system. For example, reference [28] discussed the last-mile disaster recovery for power systems, that is, how to schedule and route repair crews to restore the power network as fast as possible. Further work was developed in [29], where a randomized adaptive vehicle decomposition technique was used. However, these methods neglect the coordination of damage repair and power system restoration, which is regarded as a challenging problem [30]. Reference [30] modeled the transportation of repair crews as a vehicle routing problem (VRP) in detail and a two-stage method for the repair and restoration of distribution networks was proposed to In the previous power system post-disaster resilience improvement studies, concentrations are focused more on the restoration after large-scale blackouts, while the researches are primarily carried out from three aspects: black-start [20][21][22], network reconfiguration [23][24][25] and load restoration [26,27]. However, the interconnection of regional large power systems, reasonable protection and automatic devices, along with advanced alarm processing and fault diagnosis technology, can restrict most electrical faults and accidents in local areas [5]. Therefore, the happening probability of extreme weather events leading to large-scale blackouts is tiny. Compared with the system-wide blackout, local outages are more likely to occur. Nevertheless, the studies on the restoration of local power systems are relatively scarce and immature. Thus, it is more meaningful and practical to study the restoration of local power systems especially considering the repair of damaged components.
Currently, many studies are trying to take the damage repair work into consideration when restoring the post-disaster system. For example, reference [28] discussed the last-mile disaster recovery for power systems, that is, how to schedule and route repair crews to restore the power network as fast as possible. Further work was developed in [29], where a randomized adaptive vehicle decomposition technique was used. However, these methods neglect the coordination of damage repair and power system restoration, which is regarded as a challenging problem [30]. Reference [30] modeled the transportation of repair crews as a vehicle routing problem (VRP) in detail and a two-stage method for the repair and restoration of distribution networks was proposed to minimize the sizes and durations of outages. But the model is for distribution networks and not suitable for the generation and transmission systems. Besides, existing models cannot consider many complex post-disaster system restoration constraints precisely and synthetically or depict the dynamic changing processes of component available states. Ignorance of these constraints may make the obtained restoration scheme infeasible in reality. In addition, the model of power system restoration considering damage repair is a non-convex and non-linear optimization problem, which cannot be solved easily.
Based on the above studies, this paper proposes a novel co-optimization model to coordinate the damage repair work with system operation in post-disaster restoration for generation and transmission system. Once the damage information is available for utilities, the maintenance department can make the best damage repair route. Simultaneously, according to the anticipated available states of fault components, the dispatching department will change the system operation state. Consequently, the economic loss of outage will be minimized, with reasonable repair and system operation cost. The AC-power-flow model, the post-disaster power system restoration complex constraints and the changing processes of component states are involved systematically and integrally, which makes up the shortcomings of existing studies. Moreover, to improve the efficiency of the algorithm, Lagrange relaxation is used to decompose the model into the upper sub problem of routing repair crews and the lower sub problem of power system operation optimization. Further, the acceleration strategy is adopted to speed up the convergence of the MISOCP model. The proposed model is tested on the IEEE 57-bus test system.
The remainder of this paper is organized as follows. Section 2 describes the proposed model. Section 3 presents the mathematical formulation of the model. The LR-based acceleration strategy is discussed in Section 4. Numerical studies are provided in Section 5 to exhibit the effectiveness and application values of the proposed model. Finally, the conclusions are drawn in Section 6.

Model Description
The electrical equipment in the power system, especially buses and transmission lines, are often exposed outdoors, accounting for why they are vulnerable to the external disturbances. The damaged components will have adverse effects on the save and stable operation of power systems, leading to the occurrence of generator tripping, load shedding, system splitting or even the breakdown of power systems. Once the disaster strikes, it is critical to route repair crew to repair the damaged components, restore outage load as soon as possible and provide reliable power supply for power users.
In our co-optimization model, we propose a decision-making method for utilities to restore post-disaster power systems. After a disaster, damage assessment will be conducted to acquire the states of components, locate faults and estimate the expected time to repair (TTR) as well as the required resource number of damaged components. Afterwards the best repair scheme is made and repair teams will drive to damaged points and then fix them up based on the optimized task and route. Meanwhile, power system dispatchers will obtain the available states of system equipment and formulate the system operation mode through adjusting the output of units and switching transmission lines. The whole repair-dispatching procedure will maximally reduce the customer interruption cost, while the system operation and repair expense will be minimized.
It is worth noticing that compared to existing models, this research incorporates the second-order-based AC power flow constraints and takes the important constraints during the restoration of power systems into account, which makes the model more accurate and practical to describe the actual restoration process of post-disaster power systems.
After the acquisition of the repair and the system restoration scheme, the anticipated restoration effects can be assessed by the method proposed in [5], which helps to reflect the effectiveness and economic benefits of the strategy and compare the pros and cons of different restoration scenarios.

Mathematical Model
The main aim of the model is to restore curtailed load through repairing damaged components, starting up outage units, switching transmission lines and adjusting unit output. During the restoration process, to maintain the stability and improve the control ability of the system, generators are kept in running state. When all the faults are cleared and the power system returns to the normal state, the optimal unit commitment scheme is implemented by adjusting the start-stop status and output of units, in order to realize the economic dispatch, which is beyond the scope of this study.
It is necessary to explain the assumptions first before presenting the model.

•
All repair teams are capable of repairing any type of damage. Once a damaged component is repaired, repair team will leave for the next one immediately, until all the tasks are completed.

•
The repair time and resource to fix a fault is known and certain; the vehicle speed of repair teams is fixed.

•
During the process of damage repair and system restoration, no extra new equipment damage occurs.

•
The repair expense of fault points and the outage unit restoration cost are fixed, which have nothing to do with the repair moment or repair teams. • Generators are not damaged by disasters because they are often located indoors.

Objective Function
The objective is to minimize the power outage loss, with possible small expense of damage repair and power generation cost, as shown in (1) where C ope indicates the generation costs of all generators [31], C rep represents the damage repair cost, including labor and transportation cost and C loss is the total economic loss due to outage during the restoration horizon. α 1 , α 2 and α 3 are weight factors. As the restoration of lost load is the most urgent mission when there is power outage, the value of α 3 is correspondingly larger.
Besides, we need to point out that because the repair cost and unit restoration cost are fixed once the damaged components are identified, they are not involved in the objective function.

•
Constraints of damage repair routing Energies 2018, 11, 2611 The damage repair routing constraints [28] are displayed in this part, helping to find the optimal routes and tasks for each team. In the repair route of team c, if c leaves damaged component x for y, M x,y,c equals to 1. Constraint (2) ensures that team c will leave damaged component y once y is repaired. Constraints (3) and (4) ensure that each team starts from the repair center and constraints (5)-(7) ensure that each team finally returns to the repair center after all the missions are completed. If damaged component x is repaired by team c, N x,c equals 1. Constraint (8) ensures that each damage point will be repaired by one and only one team, avoiding the overlap of tasks. Constraint (9) builds the relationship between N x,c and M x,y,c . If team c travels from x to y, then N x,c equals to 1.

Constraints of repair resources
Constraints in this part help to reflect the resource abundance and availability, which is an important factor to be considered in the damage repair work. Constraint (10) indicates that the repair resources needed by team c to finish its assignment should not exceed its capacity limit, while constraint (11) ensures that the repair center has enough resources to serve all repair teams.

•
Constraints of damaged component states The time it takes for x to return to the normal state is composed of the waiting time, the route time of repair teams and the repair time. Constraint (12) helps to find the arrival time of team c at fault component y. Team c arrives at x at AT x,c and time RT x,c is spent to repair the damage. Then it takes d x,y time to drive to y. The big M method ensures that the constraint does not work if M x,y,c is 0. If damaged component x is not repaired by team c, N x,c equals to 0 and consequently AT x,c equals to 0, or AT x,c is not affected by (13). When FT x,t equals to 1, it means that damaged component x is repaired at t. Constraint (14) ensures that each damage is repaired once, which is essential for the model. Constraints (15) and (16) help to calculate the restoration moment of damaged component x, which is represented by the result of T ∑ t=1 tFT x,t . As mentioned above, N x,c and AT x,c are 0 if x is not repaired by team c, having no influence on (15) and (16). ε is a very small positive number. Constraint (17) helps to obtain the available states of damaged components, which is important information for system dispatchers.

Constraints of Power System Operation
• Constraints of power system safe and stable operation The power system must operate under some constraints to keep safe and stable when restored, that is, power flow, load balance and operation limits of electrical equipment and so forth. Besides, the reactive power balance and voltage security should also be emphasized on, which are the magnitude guarantee of the stable post-disaster restoration. Hence, the power system operation model is formulated based on AC power flow, which helps to allow for incorporating reactive power and voltage security. However, the AC-power-flow-based model cannot be solved efficiently yet, so in this paper, the power flow is constrained based on the second order cone formulation according to [32], which is an approximation but can produce more accurate calculation result than DC-power-flow-based model and bring more computational efficiency.
Energies 2018, 11, 2611 7 of 21 The power system operation constraints are presented. Constraints (18) and (19) indicate the active and reactive power balance at each bus and time t. Constraints (20)- (23) are the limits of bus voltage magnitude/phase and actual active/reactive power supply, which should not exceed their upper and lower bounds. The operation state BS i,t of bus i at time t helps to ensure these variables are limited to 0 if bus i is damaged or under repair.
Constraints (24) and (25) indicate the relationship between bus voltage magnitude and power flow. These two constraints are linear since V i,t 2 is regarded as a variable instead of V i,t . The big M method is used here to ensure that the voltage magnitudes of bus i and j will be irrelevant if the line i→j is out of service. Similarly, the relationship between bus voltage phase and line flow are constrained by (26) and (27) using the big M formulation. The voltage magnitude V i (c) and V j (c) are constant in (26) and (27). They can be either 1.0 or historical average values. h ij is the square of current on line i→j and constraint (28) is relaxed from equality constraint to inequality constraint for convexity of the model. The capacity limit of apparent power flow on transmission lines is ensured by (29). Moreover, the current on lines should be within the thermal current limit, as shown in (30). The operation state LS i,t of line i→j at time t helps to ensure these variables are limited to 0 if line i→j is out of service.
Constraints (31) and (33) define the active and reactive power output limits for unit g, if unit g does not malfunction or locate on damaged buses. The ramp speed of generators is limited by (32).

•
The restoration characteristic of generators The damage of bus will result in the outage of generators on it and these generators will not be started up until the bus is repaired. If SS g,tsg equals to 1, unit g is started up at tsg, as shown in Figure 2 [33]. Once the generator is started up, it has to absorb power from the power system (without considering generators with black start ability). Then the generator has the capability to transmit power to the system after its auxiliary power is self-sufficient. The simplified restoration process of generators is shown in (34)-(39).
If the generator is started up at tsg, in the interval 0~tsg, the output of the generator is 0, as shown in (34); constraint (35) indicates that the unit will absorb P g start power to start for time tdg; then the power absorbed by the unit is decreased gradually and it begins to transmit power to the system, as shown in (36). The ramp speed of generators is limited by (37). Constraints (38) and (39) define the active and reactive power output limits. If the generator is started up at tsg, in the interval 0~tsg, the output of the generator is 0, as shown in (34); constraint (35) indicates that the unit will absorb Pg start power to start for time tdg; then the power absorbed by the unit is decreased gradually and it begins to transmit power to the system, as shown in (36). The ramp speed of generators is limited by (37). Constraints (38) and (39) define the active and reactive power output limits. •

Constraints of component operation states
The operation state variable BSi,t of bus i at time t equals to 1 if the bus is neither damaged nor under repair, as shown in (40). Constraint (41) indicates that the available state variable ALSij,t of line i→j at time t equals to 1 if the line and its linked buses are normal. Constraint (42) presents that the available state ALSij,t of the undamaged line i→j is determining by BSi,t and BSj,t. As for the operation state of line i→j, LSij,t is 0 when ALSij,t equals to 0 and 1 or 0, according to the decision of system dispatchers when ALSij,t equals to 1, as expressed by (43).

•
Constraints of component operation states LS ij,t <= ALS ij,t , ∀i → j, t The operation state variable BS i,t of bus i at time t equals to 1 if the bus is neither damaged nor under repair, as shown in (40). Constraint (41) indicates that the available state variable ALS ij,t of line i→j at time t equals to 1 if the line and its linked buses are normal. Constraint (42) presents that the available state ALS ij,t of the undamaged line i→j is determining by BS i,t and BS j,t . As for the operation state of line i→j, LS ij,t is 0 when ALS ij,t equals to 0 and 1 or 0, according to the decision of system dispatchers when ALS ij,t equals to 1, as expressed by (43).
Besides, the operation state variable GS g,t of unit g at time t equals to 1 if the unit is neither damaged nor under repair, as shown in (44). If the unit is on a damaged bus, the state of the unit follows that of the bus, which is expressed by (45). In addition, constraint (46) presents the start-up state of unit g. If GS g,t is 1 and GS g,t−1 is 0, it means that the unit g is started up at time t, which is designated by SS g,t .

Coupling Constraints
The fault repair route problem and the post-disaster power system operation problem can be extracted from the above constraints. When components are damaged or under repair, they cannot be Energies 2018, 11, 2611 9 of 21 put into operation and will influence the topology of power systems. Therefore, the two problems are coupled by the following coupling constraints.
If BS x,t is 0, all lines connected to bus x are not available and the generators on bus x are out of service.

Transformation of Complex Nonlinear Constraint
The co-optimization model is established from (1)-(48). However, the model is a non-convex and non-linear optimization problem, which is very difficult to solve directly. For the 0-1 variable ALS ij,t , which is the multiplication of two binary variables, it can be linearized by (49). Thus, the complex nonlinear constraint (42) is transformed to linear constraints. As a result, the original model is converted to a MISOCP, which relatively reduces the computational complexity.

Solution
Though the original model is transformed to a MISOCP, it is still a large-scale and complex problem, which is hard for ordinary solvers to solve. The LR method is a decomposition and coordination algorithm and its computational effort increases linearly with the scale of problem, so it is suitable for the solution of the modified model. The MISOCP model can be decomposed into the upper sub problem of routing repair teams and the lower sub problem of power system operation optimization by LR and then calculated by iterative solution. Further, to tackle the concussion and slow convergence of LR, the acceleration strategy is designed to speed up the solving process.

Lagrange Relaxation of the MISOCP Model
In this part, the coupling constraints of the MISOCP model are relaxed and added into the objective function. Only simple constraints (uncoupling constraints) are left and the MISOCP model is decoupled into two sub problems by LR. The objective function is converted to the following form: which equals to λ l,t and λ b,t are Lagrange multipliers of constraints (47) and (48) respectively. It is easy to see that the former part relates to the route and working time of repair teams, while the latter part only relates to states of system components and outputs of units. The upper sub problem of damage repair routing and the lower sub problem of power system operation optimization are further formulated.

•
The power system operation optimization problem Constraints (18) Subsequently, the two problems are solved alternatively and Lagrange multipliers are updated to achieve co-optimization. In this paper, the multipliers are updated by the surrogate gradient method [34] and their updating process follows these rules.
where w up and w down are positive numbers, which represent the increase or decrease steps of multipliers separately. The value of duality gap is used as the LR convergence criterion, that is, where γ is the objective function value of the original model, which can be calculated by constraint (1); β is the objective function value of the relaxed model, which can be calculated by constraint (50).

The Acceleration Strategy
LR tends to oscillate during the convergence procedure, increasing the number of iterations and prolonging computing time. Thus, the acceleration strategy is designed to speed up the solving process, which helps to get an approximate optimal solution quickly.
The upper limit of iteration number (Nul) and the maximum number of violated coupling constraints (Nmv) should be set by decision makers. When the iteration number reaches Nul, or the number of violated coupling constraints reaches Nmv, the acceleration strategy is activated. Under this circumstance, the sub problem of damage repair routing is solved first according to Lagrange multipliers obtained from the last iteration. Then, the sub problem of power system operation optimization is solved based on the results of damage repair routing.

The Algorithm Flow
The flow chart of the acceleration algorithm proposed in this paper is shown in Figure 3. Its solution process is listed as below.

1.
The original co-optimization model is transformed to the modified MISOCP model by linearizing some complex constraints. The MISOCP model is decomposed into the upper sub problem of damage repair routing and the lower sub problem of power system operation optimization.

4.
Solve the two sub problems alternately in each iteration. If the LR convergence criterion is satisfied, output the result and end the calculation or else go to step 5.

5.
If the acceleration convergence condition is satisfied, implement the acceleration strategy and output the result. Otherwise, update the values of Lagrange multipliers and go to step 4.

The Algorithm Flow
The flow chart of the acceleration algorithm proposed in this paper is shown in Figure 3. Its solution process is listed as below. The advantages of the LR-based algorithm are as follows.
(1) The NP-hard problems are difficult to solve due to the optimization of integer variables. LR helps to distribute integer variables to different problems and enhance the solving efficiency. Besides, at the expense of increasing the number of iterations, the overall scale of the problem is reduced, making large-scale problems solvable.
(2) The algorithm makes the physical meanings of the co-optimization more concise. After LR, the model is separated into the damage repair routing problem and the power system operation problem, which are more conformable to the actual work of the maintenance department and the power system dispatching department. For utilities, once the power system is hit by disasters, the maintenance department acquires the information of components states, fault points locations and repair resources and then formulate the best damage repair routes. Power system dispatchers only have to care about the operation cost of generators, the switching of lines and the load shedding situation. Based on independent decision making, the two departments are restricting each other by The advantages of the LR-based algorithm are as follows.
(1) The NP-hard problems are difficult to solve due to the optimization of integer variables. LR helps to distribute integer variables to different problems and enhance the solving efficiency. Besides, at the expense of increasing the number of iterations, the overall scale of the problem is reduced, making large-scale problems solvable.
(2) The algorithm makes the physical meanings of the co-optimization more concise. After LR, the model is separated into the damage repair routing problem and the power system operation problem, which are more conformable to the actual work of the maintenance department and the power system dispatching department. For utilities, once the power system is hit by disasters, the maintenance department acquires the information of components states, fault points locations and repair resources and then formulate the best damage repair routes. Power system dispatchers only have to care about the operation cost of generators, the switching of lines and the load shedding situation. Based on independent decision making, the two departments are restricting each other by Lagrange multipliers to achieve the best economic benefit of the entity. The work efficiency is greatly improved as well.

Case Study
The co-optimization model is tested on the IEEE 57-bus test system [35], to demonstrate the effectiveness of the model. The computer used is with Intel Core i7-6700 3.40 GHz CPU and 8 GB RAM. The problems are modeled in MATLAB R2014a and solved by Gurobi 7.0.2.
In the case study, it is assumed that the power system is hit and disturbed by a severe typhoon. The damage locations of the post-disaster power system are known, as shown in Figure 4. The distance between different fault points, the distance between fault points and repair centers, the needed damage repair resources and time of fault points, the resource amount of repair centers and the economic loss of lost load are displayed from Tables A1-A5 in Appendix A.
The time length of restoration horizon T is 40 h. An hourly simulation step is adopted. From [36,37], the economic loss of lost load is divided into four categories: $10/kWh for the most important load, $6.979/kWh for important load, $3.706/kWh for ordinary load and $0.110/kWh for least important load. Each repair team consists of 5 members and the average wage of each member is assumed to be $70/h. The gas cost for each team is $33/100 km. The average driving speed of each team is 50 km/h. Rating of each transmission line is 100 MVA. Bus voltage limit is [0.94, 1.06]. The maximum number of violated coupling constraints (Nmv) is 6. The weights in the objective function are set to be α 1 = 1, α 2 = 1 and α 3 = 10, to highlight the importance of reducing the power outage loss.
To improve the computational efficiency, the clustering method proposed in [30] is used to acquire repair assignments of each repair center. The objective of the clustering model is to minimize the sum of the distances between repair centers and fault points which they are responsible for. It should be guaranteed that each fault component is repaired and each repair center has enough resources to finish its tasks. The clustering result is listed in Table 1, while 1 indicates the damaged component will be repaired by the repair center and 0 else. Lagrange multipliers to achieve the best economic benefit of the entity. The work efficiency is greatly improved as well.

Case Study
The co-optimization model is tested on the IEEE 57-bus test system [35], to demonstrate the effectiveness of the model. The computer used is with Intel Core i7-6700 3.40 GHz CPU and 8 GB RAM. The problems are modeled in MATLAB R2014a and solved by Gurobi 7.0.2.
In the case study, it is assumed that the power system is hit and disturbed by a severe typhoon. The damage locations of the post-disaster power system are known, as shown in Figure 4. The distance between different fault points, the distance between fault points and repair centers, the needed damage repair resources and time of fault points, the resource amount of repair centers and the economic loss of lost load are displayed from Tables A1-A5 in Appendix A.
The time length of restoration horizon T is 40 h. An hourly simulation step is adopted. From [36,37], the economic loss of lost load is divided into four categories: $10/kWh for the most important load, $6.979/kWh for important load, $3.706/kWh for ordinary load and $0.110/kWh for least important load. Each repair team consists of 5 members and the average wage of each member is assumed to be $70/h. The gas cost for each team is $33/100 km. The average driving speed of each team is 50 km/h. Rating of each transmission line is 100 MVA. Bus voltage limit is [0.94, 1.06]. The maximum number of violated coupling constraints (Nmv) is 6. The weights in the objective function are set to be α1 = 1, α2 = 1 and α3 = 10, to highlight the importance of reducing the power outage loss.
To improve the computational efficiency, the clustering method proposed in [30] is used to acquire repair assignments of each repair center. The objective of the clustering model is to minimize the sum of the distances between repair centers and fault points which they are responsible for. It should be guaranteed that each fault component is repaired and each repair center has enough resources to finish its tasks. The clustering result is listed in Table 1, while 1 indicates the damaged component will be repaired by the repair center and 0 else.

The Advantage Analysis of the Proposed Co-Optimization Model
Two cases are designed in this part to prove the advantages of the proposed co-optimization model. Case 1: the damage repair scheme and the power system operation optimization are co-optimized using the proposed method.
Case 2: the damage repair scheme is formulated independently to realize the minimization of the repair expenses. Then, according to the damage repair scheme and the available states of system components, the power system operation optimization is conducted to minimize the power generation cost and power outage loss.
The repair routes with and without co-optimization are shown in Figures 5 and 6. The arrival moments at each damaged component with and without co-optimization are listed in Tables 2 and 3.

The Advantage Analysis of the Proposed Co-Optimization Model
Two cases are designed in this part to prove the advantages of the proposed co-optimization model.
Case 1: the damage repair scheme and the power system operation optimization are cooptimized using the proposed method.
Case 2: the damage repair scheme is formulated independently to realize the minimization of the repair expenses. Then, according to the damage repair scheme and the available states of system components, the power system operation optimization is conducted to minimize the power generation cost and power outage loss.
The repair routes with and without co-optimization are shown in Figures 5 and 6. The arrival moments at each damaged component with and without co-optimization are listed in Tables 2 and 3.
The economic indices with and without co-optimization are depicted in Table 4. The power outage economic loss and the load loss with and without co-optimization are shown in Figures 7 and 8.     When the power system is restored without co-optimization, the expense of damage repair will be lower whereas the outage loss is higher. Moreover, though the restoration without co-optimization reduces the total load loss, it cannot ensure the quick restoration of the important load. By contrast, if the damage repair is co-optimized with power system operation problem, the power outage loss will be greatly reduced, especially for the most important load, which will bring more obvious social benefits and accord better with requirements of the post-disaster power system restoration. Hence, the results demonstrate that the effectiveness and necessity of the proposed co-optimization model for post-disaster power system restoration.

The Effect Analysis of the Acceleration Algorithm
The computational performance of different solving methods is shown in Table 5. It is observed that since the original model is a large-scale non-convex and non-linear optimization problem and the direct approach cannot converge to a feasible solution. Though the When the power system is restored without co-optimization, the expense of damage repair will be lower whereas the outage loss is higher. Moreover, though the restoration without co-optimization reduces the total load loss, it cannot ensure the quick restoration of the important load. By contrast, if the damage repair is co-optimized with power system operation problem, the power outage loss will be greatly reduced, especially for the most important load, which will bring more obvious social benefits and accord better with requirements of the post-disaster power system restoration. Hence, the results demonstrate that the effectiveness and necessity of the proposed co-optimization model for post-disaster power system restoration.

The Effect Analysis of the Acceleration Algorithm
The computational performance of different solving methods is shown in Table 5. It is observed that since the original model is a large-scale non-convex and non-linear optimization problem and the direct approach cannot converge to a feasible solution. Though the modified MISOCP model can be solved directly, the computation time is long, which does not meet the actual demand for rapid restoration of post-disaster power systems. However, the convergence speed of the MISOCP model can be greatly accelerated when the model is solved by the LR-based algorithm. Besides, the restoration economic loss and cost of the MISOCP model with and without the acceleration algorithm are approximately equal. That is, the acceleration algorithm can ensure the accuracy of the solution and improve the computational efficiency.

The Impacts of Damage Repair Resources Allocation and Adequacy on Restoration
The economic impact of repair resources allocation and adequacy on restoration is analyzed. The resource amount of repair centers is changed. When the resource amount available of the three repair centers is 50, 100 and 60 respectively, the power outage loss rises from $803,380 to $865,624 and the repair cost rises $53,320 from to $57,607, as shown in Table 6. What's more, some damaged components cannot be repaired within 40 h after the beginning of the restoration process and it takes longer time to restore all outage load than the original resource allocation case. Therefore, the adequacy and allocation of repair resources will have significant effect on power system restoration. The inadequate or irrational allocation of resources will increase the loss of power outage and extend the duration of power restoration. Consequently, it is crucial to identify the possible areas disturbed by disasters in advance and distribute enough repair resources to damage repair centers, guaranteeing the post-disaster fast repair of fault points.

The Impacts of Weight Factors in the Objective Function on Restoration
The impacts of weight factors in the objective function on the power system restoration is studied. Four combinations of weight factors are analyzed and the results are listed in Table 7. In the previous analysis, α 1 = 1, α 2 = 1 and α 3 = 10, as mentioned above. With this setting, the restoration work of lost load is prioritized, reducing the power shortage cost. When α 3 is reduced along with the increase of α 1 and α 2 , though the repair cost and operation cost of the restoration process have been reduced, the power outage loss has increased. Especially, if we select α 1 = 10, α 2 = 1 and α 3 = 1 to highlight the expense of power system operation, the system operation cost is decreased by about $380,000; however, the power outage loss is sharply increased by about $1,100,000. Therefore, the values of weight factors in the objective function have great influence on the restoration effect. Generally, the most critical task after disasters is to restore the outage load and decrease power loss of users, so the repair cost and system operation cost should not be important concerns for decision makers. Hence, it is essential to set larger weight for α 3 .

Conclusions
In this paper, a co-optimization model was proposed for the post-disaster power system restoration work and the enhancement of resilience. The model can optimally coordinate the damage repair work with power system operation. Additionally, the original model is converted to a MISOCP model and the algorithm based on Lagrange relaxation is proposed to accelerate the calculation speed. The case study demonstrates the effectiveness of the co-optimization model and the algorithm. Besides, the specific conclusions are drawn as follows.

•
The proposed model can support the formulation of reasonable damage repair scheme, the plan of unit output and the decisions of optimal transmission switching, to minimize the power outage loss with lower cost of damage repair and power system operation.

•
Lagrange relaxation decomposes the original complex model into two small-scale sub problems and the acceleration strategy is implemented to realize the fast solution. For utilities, the work of maintenance department and system dispatching department could be separated and Lagrange multipliers help to coordinate their work. Consequently, the work efficiency will be improved.

•
The adequacy and allocation of damage repair resource can greatly influence the restoration effect and the load loss level. The sufficient resource reserve of repair centers will significantly decrease loss in economy due to power outage.

•
To reduce the power outage loss and realize the enhancement of post-disaster power system resilience, it is significant to highlight the fast restoration of outage loads according to their importance.
Future work will further improve the computational efficiency of the developed model and test the model on real-world and large-scale systems.
Author Contributions: H.Z. and G.L. conceived and designed the study. H.Z. and H.Y. wrote the paper.