Propagation-Based Network Partitioning Strategies for Parallel Power System Restoration With Variable Renewable Generation Resources

Grid resiliency is defined as the ability to prepare for and adapt to changing conditions and withstand and recover rapidly from disruptions. Following large-scale disruptions of the power grid such as a complete blackout, parallel power system restoration accelerates the recovery process by allowing for simultaneous restoration in each island through network partitioning, thus enhancing grid resiliency. However, existing network partitioning strategies have not sufficiently considered variable renewable resources yet, and may not be flexible or computationally efficient to accommodate additional requirements from the evolving bulk power grid. To bridge these gaps, this study proposes a propagation-based optimization approach for power grid network partitioning. The requirements of blackstart resources, power balancing, synchrocheck relays for ties lines, and increased ramp rate requirements due to variable renewable resources are modeled in the constraints. Through propagation on the physical connectivity of grid assets, the network partitioning issue is formulated as a mixed integer linear programming (MILP) problem. The linearity ensures highly computational performance and optimality of the solutions from the proposed approach to identify the cutset edge possibilities for bulk power grid applications. Graph reduction strategies are proposed to further improve the computational performance for online or nearly online applications. Case studies based on two IEEE benchmark systems show that the proposed MILP model is able to determine the optimal cutset for network partitioning in approximately 0.28 seconds, and the proposed network reduction strategy is able to further improve the computational performance by approximately 9%.


I. INTRODUCTION
Efficient power system restoration following a blackout improves the electricity supply quality for over 330 million people in North America and over 7 billion people worldwide. The recent blackout in Texas has revealed certain vulnerabilities of the energy infrastructure, and some people lost power for up to 4 days in Feb. 2020 [1]. It has been reported that the entire power grid in Texas was 4.37 minutes away from a complete blackout [1]. In that event, it is estimated that it would have taken up to one full month of effort to fully restore the Texas energy infrastructure back to normal operation [1]. The devastating effects from extreme events, such as the one The associate editor coordinating the review of this manuscript and approving it for publication was S. Srivastava.
in Texas, highlight the need for the enhancement of grid resiliency [2].
According to the U.S. Department of Energy, grid resiliency is defined as the ability to prepare for and adapt to changing conditions and withstand and recover rapidly from disruptions [3]. The methods, challenges, and opportunities for measuring smart grid resiliency have been reviewed in [4]. The grid asset hardening strategies for improved grid resiliency have been studied in [5]. The improvements and challenges for a more resilient bulk power system are discussed in [6] from perspectives including grid operations, infrastructure planning, markets, and cyber and physical security. By enabling faster power recovery for interrupted customers from a given disruption, advanced power system restoration methods provide an effective way for enhanced grid resiliency. Ref. [7] provides a comprehensive review of VOLUME 9, 2021 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ the status, challenges, and opportunities for restoration of smart grids. ''Generic restoration milestones'' are proposed to outline the restoration process in [8]. Two restoration strategies, i.e., top-down and bottom-up, can be adopted by utilities for accelerated restoration of bulk power systems [9]. Bottom-up restoration is generally deemed to be more efficient for large-scale bulk power systems, and parallel grid restoration within islands accelerate the recovery processes following a large scale blackout [9]. To facilitate parallel grid restoration and safeguarding the power grid from potential blackouts, North American Electric Reliability Corporation (NERC) developed standards, including Emergency Operations Planning (EOP)-005-2 [10], which mandates an updated restoration plan for utilities and Independent System Operator/Regional Transmission Operator (ISO/RTO). Per these requirements from the NERC, each island needs to maintain sufficient resources for power network partitioning [10]. These constraints include: 1) each island shall have at least one black-start generator to ensure that the island can be restored independently and the connectivity within the island needs to be maintained; 2) each partition shall have sufficient voltage control resources to maintain system voltage stability; 3) the tie lines between disjointed islands shall have the necessary equipment for resynchronization; 4) each partition shall have sufficient power generators to maintain the frequency within an acceptable range by matching power generation and consumption; 5) each island shall be monitored by the control center to ensure its security and correct operation.
The NERC standards provide guidance to assist in parallel power system restoration, and a computational tool to automate the process for network partitioning for large-scale bulk power system applications is critically needed. The complexity resulting from thousands of grid assets in the service territory for a regional transmission operator poses a challenge for efficient network portioning. This challenge is further compounded by the changing generation mix with conventional generators replaced by variable renewable generation resources. The intermittency and uncertainties of these variable renewable resources necessitate an increased ramping capability from conventional generators to maintain the grid frequency.
To solve the power network partition issue for accelerated power system restoration, numerous technologies have been previously developed. A methodology based on Ordered Binary Decision Diagram (OBDD) is used to find possible partitioning strategies that satisfy black-start constraints and generation/load balance constraints [11]. Complex network theory is used in [12] to divide restoration subsystems in which the number of tie-lines between islands is minimized. The first eigenvector of the Laplacian Matrix of the graph representing the power grid has been used in [13] for network partitioning. Reference [14] proposed a spectral-clusteringbased methodology to form islands that have strong internal but weak external connections. PMUs are used in [15] to ensure the observability of each island. An agglomerative clustering-based network partitioning method is proposed to determine island boundaries with heuristics [16]. Other methodologies, such as mixed-integer programming [17], [18], multilevel k-means [19] and linear AC power flow [20], are also proposed for network partitioning. Utilities and ISOs, such as Pennsylvania-New Jersey-Maryland (PJM) Interconnection [21], have also developed parallel restoration strategies. These strategies are designed based on factors such as asset ownership, historical operating experience or geographic regions.
Previous state-of-the-art methodologies however have not effectively incorporated renewable energies in the grid restoration process, which may lead to the suboptimal network partitioning and subsequently impede parallel restoration efforts following a complete blackout. Moreover, existing methodologies such as the clustering-based approach or OBDD may not be sufficiently flexible to easily incorporate the new constraints from the evolving grids or grid code requirements. Even though the author's previous study in [22] has resulted in a tool to assist system planners in developing detailed power system restoration plans, numerous constraints, such as generator restarting time requirements and the transmission path search, make development of feasible restoration plans a complex process [23]. The network partitioning provides a means to facilitate the development of system restoration plans by reducing the system size, and this benefit can be better achieved when the island is optimally determined.
To bridge these gaps, this study proposes a computational methodology for network partitioning to sectionalize the entire system into different islands for parallel power system restoration. Through propagation on the physical connectivity of grid assets, the network partitioning issue is formulated as an MILP problem. The increased ramping requirements from variable renewable generation resources have been considered in the proposed algorithm. The contributions of this paper are summarized as • A propagation-based optimization model is proposed to partition the network for parallel power system restoration with an objective to minimize the power imbalance within each island while satisfying the requirements of subsystem connectivity, ramping need, synchrocheck relays for tied lines, and blackstart resources. The propagation model fully incorporates the physical connectivity of grid assets such as transmission lines, transformers and buses in network partitioning for improved computational performance.
• Variable renewable generation resources are incorporated into network partitioning. The ramping requirement and renewable generation are modeled in the constraint of the proposed optimization to form disjointed islands for parallel power system restoration.
• Network reduction principles are proposed to further improve the computational performance for network partitioning. These reduction principles include removing double lines without impacting the graph connectivity and aggregating buses when feasible. The paper is organized as follows. Section II provides the problem formulation of network partitioning for parallel power system restoration. In Section III, the propagation based optimization model is presented to search for the cutset transmission lines for network partitioning. Section IV presents the two case studies based on IEEE Benchmark systems. The conclusions are given in Section V.

II. PROBLEM FORMULATION
The power grid network is first represented by a graph. And the network partitioning problem is subsequently formulated based on the graph in this section.  Fig. 2, v 2 and v 11 are marked as red to denote that the buses are connected to a blackstart resource; v 7 and v 10 are marked as green to denote that the buses are through nodes; and v 8 is marked as grey to denote that the node is with a degree of one. For the graph, the graph diameter and degree of a vertex are defined as • Diameter L of a graph is defined to be the number of edges in the shortest path between the most distant vertices, which can be expressed as where L is the diameter of the graph quantified in the number of edges. Path {i, j} denotes the shortest path from vertex i to vertex j.
• Degree of a vertex D v i is defined as the number of edges connecting with a vertex v i , which can be mathematically expressed as where D v i denotes the degree of vertex v i ; e indicates the set of edges in the graph. e k−v i indicates the edges between vertex k to vertex v i . In Fig. 2, the diameter of the graph is 5 (for example, the path between vertices 1 and 8 is v 12 8 . Therefore, the diameter of the graph is 5). Take v 5 as an example. There are four edges {e 2 , e 5 , e 7 , e 10 } connected to the vertex. As such, D v 5 is 4. For power system restoration, only blackstart resources have self-starting capabilities, and other generators rely on cranking power from the network to restart. In the modified IEEE-14 Bus system, there are two blackstart resources, connected to Buses 2 and 11, respectively.

B. PROBLEM FORMULATION FOR NETWORK PARTITIONING
Network partitioning is used to facilitate the simultaneous restoration of multiple islands, thus minimizing the outage duration. Network partitions are subject to a set of requirements as outlined in NERC Standard EOP-005-2 [10]. Each of those main requirements are modeled as constraints in order to determine the optimal cutset edge of the graph representing the power grid network.

1) CUTSET CONSTRAINTS
Assume that there will be M islands for parallel power system restoration. In this context, the network partitioning issue is converted to identify the cutset edges so that the entire graph is separated into M disjoint islands. There will be 1) no shared vertices among two islands, 2) all vertices within an island are connected, and 3) the vertices from all islands together VOLUME 9, 2021 comprise the vertices in the graph. That is ) denotes an arbitrary path between vertices m and k within Island i.

2) BLACKSTART RESOURCE CONSTRAINTS
Assume each island has one and only one blackstart unit. Indeed, if an island has multiple blackstart units, the island may be restored in each subgroup and synchronization is carried out to integrate different subgroups, which has a similar procedure for island synchronization. Therefore, the subgroup can be regarded as a small island and it can be assumed that each island has only one blackstart unit. This can be expressed as where {v BSU } indicates the set of vertices connected with blackstart resources.

3) CONSTRAINTS OF CRITICAL LINES AND TIE LINES WITH SYNCHROCHECK RELAYS
For power system operation, some critical transmission lines denoted as {e CL } may help maintain the voltage stability of the grid. Excluding these critical transmission lines may lead to instability issues. Therefore, critical transmission lines should not be included in the cutset {e cutset } that separates the power grid network into islands. Moreover, when synchronizing different islands through tie lines, synchrocheck relays are needed and the cutset edges have to be equipped with these devices. Use {e M } to denote the set of transmission lines equipped with synchrocheck relays. These constraints can be represented as

4) CONSTRAINTS OF POWER IMBALANCE WITHIN EACH ISLAND
In order to maintain the stability of the power system under load variations, a network partitioning strategy must reserve sufficient generation power in each of the partitions and it may be desirable to find cutset edges that minimize the generation-load imbalance of the disjointed partitions as v island where d island k is the power imbalance of Island k.

5) CONSTRAINTS OF RAMPING REQUIREMENT FROM RENEWABLES
The variation of renewable generation resources necessitates the provision of the reserve and ramping capability in each island to maintain the grid frequency.  (7) In summary, network partitioning for parallel power system restoration can be formulated as determining the cutset edges of the graph representing the power network in such a way that the entire graph is separated into disjointed islands with each island having sufficient blackstart and generation resources. The issue by nature is highly nonlinear and combinatorial, which is computationally expensive.

III. PROPAGATION-BASED OPTIMIZATION FOR NETWORK PARTITIONING
After a blackout, it is critical to determine the status of electrical assets and identify the location of blackstart resources. This study assumes that such information has been available a priori. System planners/operators use this information as an input to sectionalize the entire power grid into disjointed islands. This section introduces a computational propagationbased optimization approach as a solution to the network partitioning issue. Note that the propagation in this study is more to capture the connectivity of grid assets instead of actually energizing the buses and transmission lines in the restoration process.

A. OBJECTIVE FUNCTION
The objective function is to determine the cutset edges to separate the entire graph into disjointed islands v island where d island k is the decision variable representing the power imbalance of Island k.
The min-max objective function in (8) increases the complexity for optimization, which can be converted into a single level optimization without the loss of accuracy as Subject to where d is the decision variable to characterize the potential maximum power imbalance among islands.

B. CONSTRAINTS FOR NETWORK PARTITIONING 1) PROPAGATION CONSTRAINTS FOR TRANSMISSION BRANCHES AND BUSES
Cutset edges for network partitioning have the following properties: 1) if the cutset edges are removed from the graph, the vertices within an island are still able to be propagated from the vertex representing the blackstart resources in the particular island, and 2) a blackstart resource has to go through the cutset edge in order to reach a vertex in another island. The reachability of vertices from the blackstart resources can be mathematically modeled as a propagation problem. The propagation starts with each blackstart resource, and at each step, the neighboring transmission edges and vertices can be propagated if needed. Any arbitrary two buses in the network can be reached within L propagation steps. These constraints are modeled as • If a bus can be propagated at step N, it is required that at least one of its connected transmission lines has been propagated at step N.  (14) where N BSU is the number of BSUs in the system.
• For each transmission line, it can only be propagated by at most one blackstart resource.  (16) where z BSUi linemn is a decision variable to denote if one of the buses of Line mn is propagated by BSU i. In this study, M is equal to N BSU . The cutset edges are the ones with connected buses propagated by different blackstart units, which can be mathematically expressed as Note that this constraint is not directly imposed in the optimization for network partitioning. After solving the optimization, the results are processed by applying (17) to find the cutset transmission branches.
• For a power grid, the critical lines to maintain grid stability cannot be in the cutset for network partitioning. Therefore, the two buses connected to a particular transmission line are expected to be propagated by the same blackstart resource. The constraint is mathematically expressed as

2) ELECTRICAL CONSTRAINTS IN ISLANDING
Electrical constraints such as power imbalance and system ramping requirement due to variable renewable resources are considered for network partitioning. Network partitioning results serve as an input for system operators/planners VOLUME 9, 2021 to develop power system restoration plan, i.e., provision of blackstart resources to crank non-blackstart units, sequentially energize buses & transmission lines, and pick up loads.
• The power imbalance of each island is the mismatch of the generation capacity and load demand. During power system restoration, system operators will need to pick up sufficient load to stabilize the island as the generators have a given stable operating region and to manage the overvoltage issues due to the reactive power when energizing lightly-loaded transmission lines. Therefore, a significant amount of load is expected to be restored before synchronization, and the operating levels of generation & load powers within an island need to be balanced.
where P D (v m ) and P G (v m ) are the load demand and generation at vertex v m , respectively.
• The ramping capability of conventional generators needs to meet the ramping requirement from the variable renewable generation resources within each island. m∈{1,...c,|v|} where Rr G (v m ) is the ramping capability of the generation connected to vertex v m if any and Rr REW (v m ) is the ramping need due to the variable renewable generation resource at vertex v m if any. Note that in (19), the generation from variable renewable resources can be the average value in the projected outage horizon. Current multi-day weather forecasting provides a means to estimate the available renewable generation and ramping need due to its variability even though it may have some deviation from the actual output. It has been reported that the average day-ahead forecasting error can achieve 8%-10% for wind farm and 4% for system ramp errors [24]. Therefore, averaging forecasted generation from each renewable plant can be incorporated in (19) to define power imbalance.

C. GRAPH SIMPLIFICATION TO REDUCE THE SEARCHING SPACE FOR NETWORK PARTITIONING
The optimization based network partitioning needs to model the status of each vertex and edge using binary decision variables. For a large power grid, the number of decision variables increases linearly with the number of buses, the number of transmission lines, and the diameter of the graph representing the grid network. Graph reduction cuts down the number of decision variables so as to reduce the searching space of the optimization and improve the computational performance. Some graph reduction principles are proposed and implemented as follows.

1) REDUCE DOUBLE EDGES BETWEEN VERTICES
The graph is meant to capture the spatial connectivity among nodes without explicitly considering the capacity of each edge for network partitioning. Removing the redundant edges between any two neighboring vertices reduces the number of edges while maintaining the graph connectivity. In Fig. 2, edges e 21 and e 1 between v 1 and v 2 are the double edges, which can be reduced to a single edge e 1 as shown in Fig. 3.

2) AGGREGATE VERTICES WITH DEGREE OF 1
The vertices with Degree of 1 are connecting to the neighboring vertex through a single edge. Therefore, the vertices can only be propagated by blackstart resources through the neighboring vertex for network partitioning, which lays the foundation for aggregation without impacting the partitioning results. For example, v 8 in Fig. 2 has a Degree of 1 and it can be aggregated with its neighboring vertex v 7 by removing v 8 and the connected edge e 11 as shown in Fig. 4. For the aggregated vertex after graph reduction, the property of the aggregated vertex is updated as (21) where binary status BSU (v i ) indicates if vertex v i is connected to a blackstart resource.

3) REMOVE THROUGH VERTICES
Through vertices are defined as the ones without connecting to any generators, loads, or critical lines. Therefore, regardless of the blackstart resource to propagate the through vertices, the power imbalance of each island will not be impacted. However, removing the through vertices improves the computational performance by reducing the number of decision variables in optimization. After reduction of the through vertices and the associated edges, the neighboring vertices need to be fully connected so as to maintain the connectivity of the original graph. Assume that the through vertex to be reduced has Degree of n. When performing the through vertices reduction, n neighboring edges will be removed and n * (n-1)/2 edges will be added. In this study, it is preferred to avoid an increased number of edges after removing the through vertices and the connected edges. Therefore, n * (n-1)/2≤n is expected to hold, and n is inferred to be less or equal to 3. As such, this study is focusing on removing through vertices with Degree of 2 or 3. Take v 10 in Fig. 2 as an example. v 10 is a through vertex with Degree of 2. The reduction of v 10 is shown in Fig. 5.

D. REDUCTION OF PROPAGATION STEPS FOR NETWORK PARTITIONING
In the proposed optimization-based network partitioning, the total number of propagation steps are equal to the diameter L, which is 5 of the power network graph in Fig. 2. In network partitioning for power system restoration, the propagation model is focused on the reachability of the vertices from blackstart units. The diameter for network partitioning can be redefined as the maximum length of the path from a blackstart unit to other vertices in the graph for power system restoration. The path from vertices with BSUs, such as v 2 and v 11 in Fig. 2 Fig. 2, v 2 and its connected edges {e 2 , e 3 , e 4 , e 5 } are removed as shown in Fig. 6. The longest path from v 11 to v 3 is v 11 → v 9 → v 4 → v 3 with a length of 3. The same procedure is implemented to search for the path from other vertices to v 2 . The longest path is identified as v 2 → v 4 → v 9 → v 14 with a length of 3. Therefore, the propagation step is reduced to 3 from the original diameter of the graph, which is 5 in Fig. 2. The reduced propagation step further reduces the computational complexity of the proposed optimization.

IV. CASE STUDY
The propagation-based optimization approach for network partitioning is implemented in IBM CPLEX 20.1.0.
Simulations have been conducted on a computer with 1.9-GHz Processor, 8 GB Memory on Windows 10 Operating System.

A. MODIFIED IEEE 14-BUS SYSTEM
The generation, wind, and load data for the modified IEEE 14-Bus System are given in Table 1. By performing graph reduction with the proposed principles aforementioned for the modified IEEE-14 Bus system in Fig. 2, v 8 as a through vertex is removed; v 10 as a vertex with Degree of 1 is reduced; and e 21 as a redundant edge between v 1 and v 2 is removed. The reduced graph is shown in Fig. 7 and the parameters of the graph before and after reduction are summarized in Table 1. n v and n e denote the number of vertices and edges, respectively.  v BSU = {v 2 , v 11 } and the necessary propagation steps for the graph after reduction are 3. The proposed optimization model is carried out to determine the cutset transmission lines VOLUME 9, 2021 for network partitioning. Note that the wind power is the average of the forecasted value during the power system restoration time horizon. The ramping rate of wind captures the worst ramping demand during the time horizon. In this case study, it is assumed that all transmission lines are equipped with synchrocheck relays. The proposed partitioning method identifies cutset transmission lines e 4 ,e 6 , e 7 and e 10 , which is the Cutset 1 in Fig. 7. This partitions the entire grid into two disjointed islands with the maximum power unbalance of 323.5 MW. Buses 1, 2, 3, and 5 are in Island 1, and Buses 4 & 6-14 are in Island 2. Note that the system has a total generation of 1,591.5MW while load demand is 953MW, making the system generation rich. The network partitioning results are summarized in Table 3.  A second scenario with only transmission lines e 8 , e 9 and e 10 equipped with synchrocheck relays is further studied. The proposed methodology determines Cutset 2 in Fig. 7 to partition the entire network into two disjointed islands for parallel power system restoration. The scenario with only edges e 8 , e 9 and e 10 equipped with synchrocheck relays also serves as a test case to validate the efficacy of the proposed MILP model as they comprise the only feasible cutset for network partitioning. The proposed optimization takes ∼1ms to solve in CPLEX and the feasible solution is successfully identified. The maximum power imbalance in this scenario is 567.2 MW.
For the network partitioning in this test scenario, results show that Buses 1, 2, 3, 4, and 5 are in Island 1, and Buses 6, 7, 8, 9, 10, 11, 12, 13, and 14 are in Island 2. Island 1 has the blackstart resource 1 connected to vertex 2, and Island 2 has the blackstart resource 2 connected to vertex 11. Vertices in both islands are connected within each island. These results are in line with the first requirement from the NERC standard, which requires that at least one blackstart unit is within an island. The power imbalance is minimized to ensure that each island has sufficient generation resources to maintain the grid frequency during restoration, which is the fourth requirement from the NERC standard. The two islands are interconnected by the edges of e 8 , e 9 and e 10 , which are equipped with synchrocheck relays, to ensure that the monitoring equipment for resynchronization, which is the third requirement from the NERC standard. This study assumes that the control room has full controllability of the grid and there are sufficient reactive power dispatch resources to meet the second requirement of the NERC standard. Indeed, reactive power compensation is generally a local issue. The impact of reactive power resources on the network partitioning will be further investigated in future work. By comparing the network partitioning results with the scenario of all lines equipped with synchrocheck relays, results show that the placement of synchrocheck relays has an impact on the network partitioning for parallel grid restoration.
The improved resiliency of the grid from parallel restoration can be measured as the reduced restoration time. Note that this study is aimed for network partitioning for parallel power system restoration instead of developing detailed restoration procedures. Assume that it takes a time of t 1 = 12 hours to use the blackstart unit at v 2 for grid restoration and t 2 = 16 hours for the blackstart unit at v 11 to restore the power grid. The parallel restoration times for each island are 7 hours for island 1 and 6 hours for island 2 respectively. The improved grid resiliency can be quantified by the reduced restoration time t = (t 1 + t 2 ) 2 − max (t Island1 , t Island2 ) = 14−7 = 7 hours. The detailed quantification of the improved grid resiliency will be performed in the future work with network partitioning integrated into development of power system restoration plans.

B. IEEE-118-BUS SYSTEM
The proposed approach for network partitioning was also tested on the IEEE 118-bus system. The system has 118 buses, 186 branches, 54 generators, and 91 loads. The detailed system information can be found in [25]. In this study, the generators connected to buses 31 and 87 are assumed to be blackstart units. The diameter of the graph is 14. With the system parameters, the propagation-based optimization is solved to determine the cutset transmission lines for the grid network partitioning. The scenario of an additional BSU installed at BUS 59 is also investigated. The results for the two scenarios are summarized in Table 4 and in Fig. 8. The solution time for the network partitioning is 0.279 seconds for the scenario with two BSUs and is 0.069 seconds for the scenario with three BSUs. The solution process of the optimization in CPLEX is shown in Fig. 8 for the scenario with two blackstart units. From the CPLEX solution process, it shows that the optimization converges rapidly to the optimal results.
When the system has only two BSUs at Buses 31 and 87, the optimization model identifies the transmission lines 77-82, 80-96, 96-97, 98-100, and 99-100 as the cutset edges to separate the entire networks into disjointed islands with the maximum power mismatch to be 2440.9 MW. Note that the maximum generation output from a generator is used in the optimization based network partitioning model. The edges of Cutset 1 are also highlighted in red in Fig. 9. When an additional BSU is installed at Bus 59, the transmission lines that separate the entire power grid into three islands are Lines 77-82, 80-96, 96-97, 98-100, 99-100, 15-33, 19-34, 30-38, 69-70, 69-75, 75-118, and 75-77. These form Cutset 1 and Cutset 2 as shown in Fig. 9. The results show the impact that additional blackstart resources have on network partitioning. Note that the optimization based islanding model is solved in under 1 second for both scenarios. This computational performance shows the potential of the proposed method for applications to bulk power systems that have hundreds of buses and branches.   Table 5. The proposed principles reduce the number of buses from 118 to 105, and reduces the number of branches from 186 to 169. The solution time is compared for the system with 2 BSUs with and without graph reduction. Simulation results show that the computation time reduces from 0.279 seconds to 0.254 seconds after the graph reduction, which is a 9% improvement.

C. IEEE 118-BUS SYSTEM WITH RENEWABLE GENERATION RESOURCES
Conventional generators connected at Buses 82, 85, 89, 92, 99, 100, 111, 112, and 113 are replaced by wind farms with the same generation and ramping requirement. The renewable generation and ramping requirement come from the forecasting results for the estimated restoration horizon, e.g., 24 hours. Note that the case is modified for the purpose of testing the proposed approach for network partitioning with variable renewables. With the renewable energy farms, the propagation-based optimization model is carried out with two blackstart resources installed at Bus 31 and Bus 87 respectively to determine the cutset edges for network partitioning. The solution results show that transmission lines 15-33, 19-34, 30-38, 69-70, 69-75, 75-118, and 75-77 are the cutset edges to separate the entire grid into two disjointed islands. This scenario is compared to the one without renewable energies in the system, and the results are summarized in Table 6. The results show that with conventional generators replaced by renewable energy sources, the network partitioning is significantly changed in this scenario. The primary cause of this change is that the island propagated by the blackstart resource at Bus 87 needs additional ramping generation resources to deal with the ramping need from renewable wind farms.

D. FLEXIBILITY OF THE METHODOLOGY
The proposed propagation-based optimization for network partitioning offers flexibilities for end-users. If a particular constraint does not apply or some additional constraints need to be factored in/activated, the proposed methodology is able to accommodate these needs. Take the increased ramping constraint for an example. If the end-users deem that the ramping constraint is not needed, the right side term of (20) can be replaced by a large negative number to ensure that the constraint will not bind when solving the optimization. In the meantime, if the system operators/planners prefer to control the size of the islands after network partitioning, additional constraints can be added to impose it by limiting the nearness of the transmission lines and buses to the blackstart resources within each island. For example, if the island wants to limit the grid assets to be within N Line lines from the blackstart resources, the added constraint can be expressed as where Path k, v BSU i is the distance from Bus k to blackstart unit i quantified by the number of transmission lines in between.
For a given network, the graph can be searched offline to quantify the distance from each bus to the blackstart resources, and the resulting information can serve as an input for (22) for network partitioning. The constraint is able to help limit the size of islands after network partitioning.

E. COMPARISON WITH THE STATE OF THE ART METHODOLOGIES
The proposed approach for network partitioning has been compared to the state of the art methodologies in terms of computational complexity, optimality, and flexibility. The comparison results are summarized in Table 7. The proposed methodology has a high computational performance resulting from the linearity of the optimization model, the graph reduction strategies, and the propagation principles on physical connectivity of grid assets. Through formulation of the network partitioning issue as an MILP problem, the optimality of the solution is achieved. In addition, the proposed approach is easily extended to incorporate other constraints such as controlling of the island sizes or deactivating constraints for improved flexibility for the end-users.

V. CONCLUSION
This research proposes a new propagation-based optimization approach for network partitioning in parallel power system restoration. The optimization model is designed to determine the cutset transmission lines to separate the entire power grid into disjointed islands while meeting the requirement of maintaining sufficient blackstart resources, network connectivity, power imbalance, ramping requirement, and synchrocheck relays for tie lines. The proposed network partitioning approach fully takes advantage of the physical network connectivity of grid assets for propagation from the blackstart resources to determine the cutset for improved computation performance. Principles for graph reduction are proposed to further reduce the calculation complexity. The simulation results with the IEEE 14-Bus and 118-Bus systems with different scenarios validate the efficacy of the proposed approach for network partitioning for parallel power system restoration. The following conclusions are arrived based on the case studies: • The proposed approach is able to identify the cutset transmission lines for network partitioning within a second (0.279 seconds for IEEE 118-Bus System).
• The proposed graph reduction strategy improves the computational performance through reducing the solution time by 9% in the second case study; • The proposed MILP model is able to achieve the optimal solution with its linear formulation of the problem. The cutset transmission lines leading to the minimal power imbalance are identified. The study provides a solution for network partitioning for power grids with and without variable renewable generation resources. The optimal network partitioning accelerates power system restoration for fast power recovery to enhance grid reliability and resiliency.