Since developing partitioning strategies of a power system is to classify buses or select tie lines, the power system can be simplified by applying the graph theory. According to the topological and physical characteristics, a complex power system is abstracted as an undirected sparse connected graph by adopting the following simplified standards Sun et al., 2016:

1 The transmission network in a power system is considered, whereas distribution networks are treated as equivalent loads.

The calculation of edge weight includes P _ij and M _ij .

P _ij is the transmission active power of edge i − j and can be calculated from the power flow or obtained from the pre-blackout monitor. In power networks, edges with vast transmission power are the critical edges. Hence, the physical characteristic of power systems can be evaluated by the active power value of edges.

The edge betweenness M _ij is usually used to evaluate the hub of edge i − j in a network Lin et al., 2011 The larger M _ij of edge i − j, the greater role that edge plays in the network. It reflects the importance of edge i − j in the network and can be calculated by M i j = ∑ i ≠ j σ i j e σ i j (1) where σ _ij denotes the number of shortest paths from node i to node j and σ _ij(e) represents the number of shortest paths from node i to node j when edge i − j is removed.

Thus, considering the topological and physical characteristics of edges for partitioning, this paper improves the weight B _ij of edge i − j as: B i j = ω 1 ⋅ | P i j | + ω 2 ⋅ M i j (2) where ω ₁ and ω ₂ are the weight coefficient, which can be obtained by the entropy weight method.

The node importance is assessed by S _i , S ₀, and N _di .

S ₀ is the base capacity. S _i is the sum of the rated capacities of all connecting edges between nodes i and j, and S _ij is the rated capacity of all connecting edges between nodes i and j.In power systems, the node capacity ratio s _i is a significant physical factor for evaluating node i’s importance. The large value of s _i means that node i has a high transmission power ratio and plays an important role in the power system. It can be obtained by the power flow calculation or the pre-outage monitor, i.e., s i = S i / S 0 (3) S i = ∑ i , j S i j (4)

N _di represents node contraction. In power systems, N _di refers to the degree of contraction between nodes connected in the power system. The larger the N _di , the greater the competitiveness of nodes in the network, reflecting the importance of nodes in ensuring network connectivity in the power system. N _di involves the connectivity and topological structure between nodes. It can be obtained by: N d i = 2 T i k i k i − 1 (5) where T _i denotes the actual number of edges between the neighboring nodes of node i, and k _i represents the number of neighboring nodes of node i.

Thus, considering the topological and physical characteristics of nodes, the improved node importance of node i is described as: C = ω 3 ⋅ s i + ω 4 ⋅ N d i (6) where ω ₃ and ω ₄ are the weight coefficient, which can be obtained by entropy weight method.

In general, B _ij reflects the tightness and the energy flow capacity between the corresponding nodes. It is used to calculate the weighted adjacency matrix W . C _i considers the node degree and the power capacity. It is used to compute the node subsystem influence indicator K _is . Then, W and K _is are utilized for iterative computation of the label propagation matrix P . Hence, the topological and the physical characteristics of power systems are considered comprehensively.

By using the label propagation algorithm, some buses occur oscillation, which results in the inability to determine the boundaries of subsystems. Which subsystem the oscillation bus should be divided into can be modeled as a cooperative game problem.

Cooperative game theory can determine the cooperation between players and coalitions (groups of players). The cooperation between players and coalitions can in-crease the coalitions’ total profit. A hierarchical clustering algorithm for the cooperative game was proposed to reveal the network community Zhou et al., 2015. A community formation game-theoretic framework with the gain function and the loss function was proposed to evaluate the possibility of nodes joining a community Chen et al., 2011.

In a coalition, the profit is allocated to different players according to the distribution agreement. Shaply value is the only definite solution for transferable profit with coalitions. It provides a unique and fair solution to the cooperation among the players in a cooperative game, considering the relative importance of each player to the game. Thus, the Shapley value of a player accurately reflects the marginal value that the player contributes to the game. The larger the Shapley value of the player, the greater contribution of that to the coalition.

A cooperative game strategy C _g =(n, v) consists of n players and a characteristic function v. The Shapley value of player i in coalition s, i.e., X (s,i), is represented as: X s , i = ∑ W | s | v s − v s \ i (12)

W | s | = n − | s | ! | s | − 1 ! n ! (13) where |s| is the number of players in the coalition s; v(s) is the cooperation profit of coalition s including player i, whereas v(s\i) is that without player i; W (|s|) is the average contribution factor. “!” represents the calculation of factorial. The more contributions of player i, the more profits can be obtained by player i. Thus, the Shapley value of play i can be used as the average of the marginal contributions of player i joining in each possible coalition.

For the partitioning problem of a power system in this paper, buses are regarded as the players, and subsystems are treated as the coalitions. Each bus tries to join in each subsystem to cooperate with the other buses in the corresponding subsystem. The Shapley value represents the closeness of a bus and the other buses in the same subsystem. Thus, it indicates the contribution of buses for system partitioning. If the Shapley value of a bus in one of all subsystems is larger than that of the bus in the other subsystems, the bus should be divided into the corresponding subsystem.

According to the partitioning result in each iteration of labels propagating, this paper uses the dynamic load rate R _i,s and the node contraction N _di for evaluating the profit of each subsystem: v s = 0 , if  | s | ≤ 1 ∑ i ∈ s ∑ j ∈ s , j ≠ i f i ⋅ a i j , if  | s | ≥ 2 (14) where a _ij denotes the connection between buses i and j. If there is an edge between i and f, a _ij = 1; otherwise, a _ij = 0. The importance degree f(i) of bus i in subsystem s is defined as: f i = R i , s + N d i (15) R i , s = P D i L s (16) where L _s is the total loads in subsystem s; P _Di is the load at bus i. If bus i is not in subsystem s, R _i,s equals 0.N _di can be obtained by Formula 5.

Here is a small test system for explaining how boundaries of subsystem are determined by the Shapley values. Figure 1 shows that the IEEE 9-bus test system is divided into two subsystems, including 3 generators, 9 buses, 3 transformers, and 6 lines. According to the label propagation algorithm, buses 1, 3, 4, 5, 6, and 9 are located in subsystem S1, while buses 2, 7, and 8 are located in subsystem S2. The X values of buses are X 1 = 1,1.185 , 3.105 , 1.105 , 2.185 , 2.37 and X 2 = 1.395 , 2.79 , 1.395 , respectively.

In the next iteration of label propagation, the X values of buses 1, 2, 3, 4, and 6 do not increase further. However, if bus 9 is divided into subsystem 2, its X value is 1.175, which is smaller than that of 2.37 in subsystem 1. Thus, bus 9 should return subsystem 1. Moreover, if bus 5 is divided into subsystem 2, its X value is 1.29, which is greater than that of 1.105 in subsystem 1. Thus, bus 5 should be divided into subsystem 2. Ultimately, the X values of all buses do not increase with the propagation of labels. The partitioning strategy is S1 = 1,3,4,6,9 and S2 = 2,5,7,8 , as shown in Figure 2.

Based on the label propagation algorithm and the cooperative game theory, this paper proposed a novel partitioning method for parallel restoration. The label propagation algorithm with the characteristic of low computational complexity can reduce the solution space of the cooperative game, and avoid huge subsystems or meaningless subsystems. Especially, the Shapley value of buses can deal with the label oscillation, which is a common disadvantage of the label propagation algorithm. Hence, the proposed partitioning method for parallel restoration has low complexity, excellent classification ability, and fast convergence speed.

The proposed partitioning approach is shown in Figure 3, including 6 steps.

Step 1

Initialize the edge weight B _ij , node importance C _i , and weighted adjacency matrix W . Data about the topology and the operating condition before the blackout are input to construct the label propagation matrix P ₀ and the subsystem discriminant matrix F ₀. Initial subsystems are determined based on the BS units.

Step 2

Initialize the X value of all buses. In the initial stage of labels propagating, only the buses of BS units have the subsystem labels. It means that each initial subsystem has a bus. The BS unit-bus in each subsystem cannot cooperate with other buses out of the corresponding subsystem. Thus, their X values are 0, i.e., X(s0,i) = 0. For other buses, they do not belong to anyone subsystem. Hence, the initial X value of every bus except the BS unit-bus is also equal to 0.

Step 3

Construct the label propagation strategy base on Formula 11.

Step 4

A bus marked by the label can cooperate with each subsystem in each label propagation process. The X value of each bus is calculated based on Formula 12. According to the maximum X value of a bus cooperating with different subsystems, the affiliation of the corresponding bus is determined.

Step 5

The feasibility of the generating partitioning strategy is evaluated by partitioning constraints, including BS constraint, power balance constraint, subsystem size constraint, observability constraint, and network topology connectivity constraint. If all constraints are met, go into step 6; otherwise, return to step 3 for the next iteration without updating P .

Observability constraint: Since the system condition should be monitored during the process of power system restoration, all buses should be observed by the phasor measurement units (PMUs). Although the outage system is divided into multiple subsystems for parallel restoration, all buses in each subsystem should be also observed by PMUs. The observability constraint is ∑ z i j r j + r i ≥ 1 (17) r i = 1 , if there exists a PMU at  i 0 , otherwise (18) where r _i indicates whether there is a PMU at bus i. If edge i-j is within the subsystem, z _ij = 1; otherwise, z _ij = 0.

Network topology connectivity constraint: Each generated subsystem should be a connected network; otherwise, it cannot be called a subsystem. J = 0 (19) J = F ⊕ F ∧ F A F T (20) where J is an S× S matrix. A = [ a i j ] n × n represents the adjacency matrix of the network. If bus i and bus j are connected, a _ijn×n = 1; otherwise, a _ijn×n = 0. ⊕ represents the logical operation of exclusive disjunction and ∧ is the logical operation of conjunction.

Subsystem size constraints: Parallel restoration is to partition the outage system into multiple small subsystems for service restoration. The effect of parallel restoration relies on the subsystem with the longest restoration time. Huge subsystems need plenty of time to restore loads. Thus, to avoid the occurrence of huge subsystems, the subsystem size should be limited: max m l i ≤ U max (21) where m _li represents the shortest path from the black start power source to bus i, U _max is the maximum allowed size of the subsystem, and N _D is the set of nodes updated in each iteration process.

Power balance constraint: For a generating unit, the minimum output active power is one of the key characteristics and is determined by technical conditions. Generating units cannot be operated stably at an output lower than their minimum output power. Thus, the sum of the generating units’ minimum output in each subsystem should be less than the total load. ∑ i = 1 n α P G i − ∑ i = 1 n P D i ≤ 0 , i = 1,2 , … , n (22) where α represents the minimum technical output coefficient of unit at i, P _Gi denotes the rated output power of unit at i, and P _Di represents the load at i.

Step 6

If F _t is converged, output the node sets for each subsystem; whereas, update P _t and return to continue iteration in step 3.

The IEEE 39-bus system includes 10 generators, 39 buses, 12 transformers, and 34 lines Zimmerman et al., 2011. BS units are located at buses 30, 31, and 34, respectively. The initial X values of buses 30, 31 and 34 are 0. Thus, this system is divided into three subsystems. Buses 30, 31, and 34 are in subsystems 1, 2, and 3, respectively. The constraint on the subsystem size is set as U _max = 8 Li et al., 2019. The optimal PMU locations are shown in Figure 4 Demetriou et al., 2019. The minimum technical output coefficient α of the hydropower unit is 0, while that of the thermal power unit is 0.35. The values of edge-weight in (2) and node-weight in (6) are calculated by the entropy weight method. That is, ω 1 = 0.4096, ω 2 = 0.5904, ω3 = 0.6750, and ω4 = 0.3250. The partitioning strategy by using the proposed method is shown in Figure 4.

From the 1st step to the fifth step iteration, the three types of labels, representing three subsystems, propagate from buses 30, 31, and 34 to other buses. The X values of all buses are increasing. In the 6 ^th label propagation, the X values of all buses except bus 9 do not increase. In this step, the X value of bus 9 in subsystem 1 is 1.009. If bus 9 is divided into subsystem 2, X (2, 9) is 1.032, which is greater than X (1, 9) = 1.009. Therefore, bus 9 should be divided into subsystem 2. In the end, all buses’ X values do not increase as the labels propagate.

Edges 3–4, 9–39, 14–15, and 16–17 are the tie lines of three subsystems. The buses in each subsystem are in Table 1. The X values of some buses are shown in Table 2. The sizes of subsystem 1, subsystem 2 and subsystem 3 are 5, 4, and 6, respectively. The size of each subsystem is roughly balanced.

Since the power system has obvious small-world network characteristics, the index of modularity, namely Q, can be used to evaluate the rationality of the generating subsystems Newman and Girvan, 2004. In order to reflect the electrical connection of the power system, this paper takes the active power of buses as the weight to construct the weighted graph. For an actual network, the value of Q is generally between 0.3–0.7 Lin et al., 2011.The Q value of this partitioning strategy is 0.638. It indicates that the partitioning strategy is reasonable. Figure 5 shows that each subsystem has sufficient generation capacity to match the load consumption. All buses in each subsystem are entirely observable. Thus, the three formed subsystems in this strategy satisfy the partitioning constraints. The partitioning strategy obtained by this method is beneficial to the parallel restoration after a power outage.

Figure 6 and Table 3 shows the comparison of the proposed method and the other three methods in Lin et al., 2011, Quirós-Tortós et al., 2014; Sun et al., 2016 References Lin et al., 2011and Quirós-Tortós et al., 2014 both employed algorithms based on complex network theory, focusing solely on the partitioning from the network’s topological characteristics. Reference Lin et al., 2011 utilized the GN algorithm to obtain the partitioning results with relatively high modularity, while reference Quirós-Tortós et al., 2014 used the cut-set theory to achieve the partitioning strategies with fewer interval connection lines. Reference Sun et al., 2016, on the other hand, employed the mathematical optimization method, considering power flow calculations, to obtain the partitioning results with the minimal interval transmission exchange power. In Figure 6, references Lin et al., 2011 and Quirós-Tortós et al., 2014represent comparison method 1, while reference Sun et al., 2016 represents comparison method 2. The number of tie lines obtained in this paper is the same as Lin et al., 2011; Quirós-Tortós et al., 2014, but smaller than that of Sun et al., 2016. A small number of tie lines indicates that the complexity of reconnecting subsystems is low. It helps improve the success of the third stage of parallel restoration. The Q value of the partitioning strategy in this paper is 0.6380, which is the largest Q of the strategies obtained by using methods in Figure 6, indicating that the partitioning characteristics of the strategy developed by the proposed method are obvious. The small power exchange among different subsystems helps the power balance of each subsystem Qing, 2015, Demetriou et al., 2019. The power exchange between subsystems is less, the connection between subsystems is weaker. The absolute value of power exchange among subsystems obtained by the proposed method is 339.04 MW, the smallest value of the three strategies. Overall, the proposed method has a good performance of system partitioning.

The IEEE 118-bus system has 54 generators, 9 transformers, 118 buses, and 186 branches. The BS units are located at buses 12, 25, 55, 66, and 100, respectively. The initial X values of buses 12, 25, 55, 66, and 100 are equal to 0. The IEEE 118-bus system is sectionalized into 5 subsystems. Buses 12, 25, 55, 66, and 100 are in subsystems 1, 2, 3, 4, and 5, respectively. The optimal PMU locations for this test system are shown in Demetriou et al., 2019.

The label oscillation among buses in the label propagation algorithm is shown in Table 4. For instance, buses 19, 30, and 31 are marked by subsystem 1-based label or subsystem 2-based label repeatedly. The label propagation algorithm is not converged, and the partitioning strategy cannot be generated. Thus, the buses with label oscillation should be merged into suitable subsystems to develop the final partitioning strategy.

The X values of some buses are shown in Table 5. For bus 74, it is marked by subsystem 2- or 5-based labels repeatedly. Thus, there are two types of label oscillation in bus 74. By applying the proposed cooperative game strategy, X (2, 74) and X (5, 74) are 0.998 and 0.912, respectively. Since X (2, 74) is larger than X (5, 74), bus 74 is merged into subsystem 2. The other buses with label oscillation are cooperated with each subsystem and divided into the subsystem with the maximum Shapley value.

Figure 7 shows the final partitioning strategy. Edges 8–30, 17–30, 17–31, 17–113, 19–20, 30–38, 34–43, 37–38, 37–39, 37–40, 49–54, 49–51, 50–57, 59–60, 59–61, 63–64, 69–70, 70–75, 74–75, 68–81,69–75, and 69–77 divide the system into 5 subsystems. The absolute value of power exchange among subsystems is 1474.03 MW. The number of buses in each subsystem is not much different. The sizes of subsystems 1, 2, and 5 are all equal to 5, and that of subsystems 3 and 4 are 3 and 4, respectively, as shown in Table 6. The size of each subsystem is roughly same. The Q value of the partitioning strategy is 0.605, indicating that the partitioning characteristics of the strategy developed by the proposed method are obvious. Figure 8 shows that each subsystem has sufficient generation capacity to match the load consumption. All the generated subsystems are completely observable. The results indicate that the method proposed in this paper is effectively applied to the partitioning problem of large power systems.

An urban power grid in South China is used to verify the effectiveness of the proposed method on practical power systems. Only 500 kV and 220 kV networks are considered, as shown in Figure 9. This grid has 7 power plants with a total installed capacity of 5100 MW. The total load is 3507.7 MW. The BS units are located at plant-buses 10, 22, and 29, respectively. The initial X value of each bus is 0.

By using the proposed method, the grid is partitioned into 3 subsystems. The tie lines are edges 2–3, 2–6, 7–35, 9–17, 20–27, 20–31, and 32–27, as shown in Figure 9. The Q value of the partitioning strategy is 0.413, indicating that the partitioning characteristics of the strategy developed by the proposed method are obvious. The size of subsystem 1 is 5, subsystem 2 is 3, and subsystem 3 is 5respectively, as shown in Table 7. The size of each subsystem is almost the same as that of others. Each subsystem contains a BS unit and has sufficient generation capacity to match the load consumption. The results indicate that the method proposed in this paper can be effectively applied to the partitioning problem of practical power systems.