Optimal controlled islanding considering frequency-arresting and frequency-stabilising constraints: A graph theory-assisted approach

The optimal controlled islanding of power systems is a practical solution to prevent system blackouts if the boundary of islands and corrective control actions in each island are carefully speciﬁed. This paper establishes a model for the controlled islanding problem by a proposed mixed integer linear program (MILP). The frequency-arresting (FA) and frequency-stabilising (FS) constraints are linearised and incorporated in our FA-FS-constrained MILP model to prevent triggering of load shedding (LS) relays. This achieve-ment makes our model capable of handling low-inertia networks. Intentional LS and stepwise generation curtailment are corrective actions accommodated for frequency control and power mismatch considerations. These corrective actions increase the degree of free-dom and broaden the feasible space of the MILP model under the envisaged tough operating conditions. Our model’s computational efﬁciency is improved using a proposed graph theory-based network reduction technique. The basic groups of coherent generators, determined in the ofﬂine mode, are aggregated as equivalent buses using an extended Steiner tree method. A graph-path determination technique is also proposed to generate disconnection constraints (between equivalent buses of incoherent areas) and bus-allocation constraints. Simulation results on the IEEE 39-bus test system and a 76-bus case study verify the proposed network


Motivation
Power systems, as large-scale dynamic systems, continuously deal with different types of disturbances in their operation. Some of these disturbances are severe and may lead to cascading failures and, eventually, power-system blackouts with significant socio-economic consequences [1]. These blackouts are mostly the result of uncontrolled islands' forced creation with unexpected power imbalances and the loss of synchronism [2]. When a power system is severely disrupted, a series of control measures are needed to prevent the disruption from spreading and the network from collapsing [3]. The August 14, 2003  while there were not enough situational awareness and control actions [4].
Optimal controlled islanding (optimal partitioning of transmission network) is capable of preventing power system blackouts. Optimal controlled islanding can preserve a significant part of the network load and accelerate the system restoration process. The researchers have made extensive efforts to achieve proper islanding strategies based on real data from various power systems worldwide [5][6][7][8].

Review and classification of relevant literature
The controlled islanding procedure has three stages: (i) Proper timing for its initiation (when), (ii) determination of the boundary of the islands (where), and (iii) appropriate corrective control actions in each island to keep it stable (how). For the first stage, some expedient researches resolve the timing problem of controlled islanding [9][10][11][12]. For the second and third stages, optimisation models are proposed to find optimal islands' boundaries and the proper control actions. Controlled islanding optimisation problems are often nondeterministic polynomial time (NP)-hard optimisation problems that include binary variables and non-convex constraints. Consequently, optimal controlled islanding models are computationally challenging, especially for online applications. An optimal controlled islanding model confronts with the following three challenges: (i) How to reduce the underlying transmission network in the model to limit the number of constraints and decision variables while maintaining the solution precision as much as possible? (ii) How to consider the transmission network stability in the view of arresting and stabilising system-frequency after network partitioning? (iii) How to improve the computational performance making the model suitable for online applications?
As per the above challenges, we have prepared a relevant literature review: A graph-based ordered binary decision diagram for network reduction was proposed in [13,14]. Based on the hierarchical spectral clustering method, reference [15] performed a controlled islanding algorithm considering the power network's weighted equivalent graph. The graph spectra methods for controlled islanding of conventional and lowinertia power systems were developed in [16,17]. Reference [16] employed an eigenvector to establish the final islands with coherent groups of generators. The proposed method in [16] was developed to find a reliable index to identify coherent generators, especially for low-inertia resources. In the low-inertia systems, frequency stability issues are much more critical for performing an islanding process. These issues were not incorporated in [16]. A slow coherency technique was proposed in [18], where an eigenvector-based calculation was used for assigning non-generator buses to coherent areas. A multi-attribute approach was developed in [19], which uses the electrical distance concept along with an evolutionary method to extract feasible islanding solutions.
In [20], a mixed-integer linear program (MILP) model with linearised AC power flow and reactive power constraints was suggested for the optimal controlled islanding problem. The alternative MILP model for optimal controlled islanding in [21] has been solved using a search space reduction approach and a recursive process to minimise the power flow disruption. The presented MILP approach in [22] results in multiple valid solutions adapting to various operation conditions, including internal connectivity constraints. As per restoration issues, an MILP-based controlled islanding approach was presented in [23]. Also, the islanding problem is directly linked to the restoration problem by a real-time state estimator. Neither of these approaches has considered the frequency stability of the islands.
As reported in [24], frequency stability constraints can be considered in the controlled-islanding problem's MILP optimisation framework. The proposed method in [25] added transient stability constraints of the obtained islands using equal area criterion. The stepwise generation curtailment (SGC) as an effective and practical means for power mismatch relief is not incorporated in [24,25].
More research papers in the area of controlled islanding can be included in this section. However, the summary of findings obtained from the mentioned literature review in terms of conceptual, technical, and practical is as follows: (i) Given that the controlled islanding problem is a timeconsuming procedure, the graph theory algorithms based on the equivalent graph of power network lead to framework improvement and computational enhancement. (ii) As per post-islanding stability challenges, frequency control constraints, the coherency-based grouping of generators, and necessary control actions are the most reasonable means to achieve stable islands. (iii) Proper distribution of the computational burden of the controlled-islanding approach in offline and online stages makes it more viable and adapt to wide-area monitoring, protection, and control (WAMPAC) capabilities. (iv) In modern power networks with high penetration of renewable resources, including low-inertia issues in the islanding problem produces more pragmatic islands.

Contributions
This paper proposes an MILP framework for optimal controlled islanding where frequency-arresting (FA) and frequencystabilising (FS) constraints are explicitly modelled. Accordingly, our FA-FS-constrained MILP model can also be used in lowinertia networks. We assume that basic coherent groups of conventional generators are successfully identified employing an effective offline method. Also, final coherent groups, including low-inertia resources, are specified by an online coherency detection index. For computational enhancement, each basic coherent group is aggregated as an equivalent bus using an extended version of the Steiner tree algorithm based on the minimum spanning tree concept. A graph-path determination based on Yen's k-shortest paths algorithm is proposed to produce basic paths forming the disconnection constraints and bus-allocation constraints. Thus, the main contributions of the current paper are: (i) We propose a comprehensive FA-FS-constrained MILP model for optimal controlled islanding of transmission networks. The explicit modelling of FA constraints makes our model applicable to low-inertia networks. Accommodating the LS and SGC in our proposed model makes it flexible for severe operating conditions. (ii) A graph theory-based network reduction technique is proposed to improve the proposed FA-FS-constrained

Paper structure
The following sections are arranged as follows. Section 2 establishes the proposed method's outline in seven steps, including four offline steps and three online steps. Section 3 presents the graph simplification technique to network reduction based on an aggregated bus for each basic coherent group's generator buses. Section 4 proposes the path-determination algorithm to identify relevant paths between aggregated buses of incoherent areas as well as non-generator buses and aggregated buses. In Section 5, the FA-FS-constrained MILP formulation is prepared, which minimises the total amount of LS subject to frequency control, disconnection, power flow, bus-allocation, and other related constraints. Section 6 conducts simulation results of the proposed approach, including modelling and time-domain results. Finally, Section 7 concludes the paper.

OUTLINE OF THE PROPOSED METHOD
As shown in Figure 1, the proposed approach is organised in seven sequential steps. The first four steps are offline and can only be carried out once for a network with a given topology. The next three steps are online, making the proposed method applicable in WAMPAC systems.

2.1
Step 1: Coherency detection Following the methodology proposed in [26], the basic coherent groups of generators are determined offline. As a technically sound assumption, certain coherent groups of generators independent of location and severity of disturbances are used [27]. In other words, the most influential factors in determining the basic coherent groups are system topology and loading conditions.
In this step, depending on how the network is configured and loaded, basic coherent groups of generators are achieved. Moreover, this step's framework is founded on an offline contingency analysis using a wide-area measurement system. This step's input database can also be updated to adapt the results with new data. It means that in a given network, the input database is enhanced over time and can result in adaptive outputs based on the most recent topological and operational data. Accordingly, an updated input database is applied to extract the basic coherent groups utilising the proposed approach [26].
In summary, this step is implemented in two stages. First, preparing the updated input database, and second, applying the method of [26] to achieve the basic coherent groups of generators. In online conditions, the last updated basic coherent groups are used, which are already prepared.

2.2
Step 2: Equivalent graph representation The graph theory-based techniques have received growing attention for reduction, simplification, and aggregation of power systems. The weighted adjacency matrix of the equivalent graph is constructed as follows: (i) If there is a line between two buses, the corresponding adjacency matrix element is the amplitude of the respective line impedance as a measure for the electrical distance between buses [28]. In this paper, two graph algorithms, including an extended Steiner tree algorithm and path-determination algorithm, are employed. In these algorithms, all possible connectivity of related buses is evaluated. In the case of parallel transmission lines, the equivalent impedance is used. (ii) Diagonal matrix elements and those non-diagonal elements with no direct connection are set to be infinite.
Accordingly, the adjacency matrix for a non-directional graph representation of a power network is symmetrical. This matrix is employed as an input for graph simplification and graph-path determination algorithms. In a regular test system, this offline step is almost an effortless step that does not considerably affect offline activities' overall computational time.

2.3
Step 3: Graph simplification The original network graph is simplified in this step using graph reduction algorithms. The output of this step is a simplified graph and its adjacency matrix. In this step, using the original adjacency matrix of the network, an extended Steiner tree algorithm is proposed. A searching framework and a minimum spanning tree technique from [29] conclude an aggregated bus for each generators' basic coherent group. Moreover, there are no substantial deviations of power flow results for the reduced network than the main network, applying the least electrical distance concept. In our case studies, the total computational of this offline step is less than 3 min. This step will be elaborated further in Section 3.

2.4
Step 4: Graph-path determination In this step, the way to extract the basic paths as a database in offline mode is explained. This step's output is the basic paths database, employing the simplified graph's adjacency matrix achieved from the previous step. This step's methodology is a graph-path determination procedure that applies Yen's kshortest path algorithm [30]. A minor modification on Yen's algorithm has generated all basic paths between aggregated buses of incoherent areas as well as each aggregated bus and non-generator buses. These paths are essential to precisely formulate disconnection and bus-allocation constraints. It may be noted that this procedure finds the least required paths to limit the number of relevant constraints and variables. Detailed discussion on step 4 will be then given in Section 4.

Step 5: Updating coherent areas and paths
This step is the first online step, and it requires an online coherency detection method, such as the method proposed in [31]. This method employs a real-time index defined over the correlation coefficients of rotor angle swings of network generators. This index is applied to assess the degree of coherency for any combination of network generators by a predetermined setting in a comparative manner. In this paper, this index is utilised to identify the possible combinations of basic coherent groups of generators using online data from wide-area measurements. Naturally, changing the degree of coherency setting can change the final coherent groups of generators. In this paper, we have used the practical setting proposed in [31]. Generally, decreasing the threshold setting leads to more sensitivity and, of course, fewer islands. This step's output is a subset of basic paths that contain final paths corresponding to disconnection and bus-allocation constraints. According to the power network's loading and topological conditions, some of the basic coherent groups of generators could be merged to form larger coherent groups just before the islanding. Hence, paths between these merged groups are removed.
Moreover, in this step, low-inertia resources (such as wind power (WP) stations) are allocated to related basic coherent groups using [31]. Accordingly, each low-inertia resource's output power is added to the coherent group's pertinent aggregated bus' generation power. Besides, bus-allocation binary variables of low-inertia buses are fixed so that these buses indeed lie in their coherent islands. Besides, those paths that are included in switched-off transmission lines would also be removed from the final paths. Briefly, this online step does three main tasks: (i) Finding merged basic coherent groups of conventional generators. (ii) A filtering process on basic paths to acquire final paths. (iii) Allocating low-inertia buses to the relevant coherent groups.
This step is very fast in terms of computational time, which does not affect the proposed approach's online computations requirements.

2.6
Step 6: Proposed FA-FS-constrained MILP Our proposed FA-FS-constrained MILP model minimises the LS objective function subject to a suite of constraints such as power-flow, disconnection, and FA-FS constraints. The FA-FS constraints are conducted in such a way that they incorporate low-inertia issues and include frequency deviation (FD) and rate of change of frequency (ROCOF) constraints. Solving the proposed FA-FS-constrained MILP model yields an islanding solution to implement (step 7). The problem formulation and supplementary discussions are given in Section 5.
In our case studies, the overall computational time for online actions, including 100 ms for transmission switching operation, is less than 1 s.

2.7
Step 7: Implementation The optimal solution of the proposed FA-FS-constrained MILP contains essential information to execute the islanding process. Naturally, the islanding process is executed in two phases. In the first phase, those lines that need to be switched off to isolate the incoherent areas will be tripped. The second phase comprises corrective control actions, including LS and SGC, to stabilise the frequency in the permissible range. Accordingly, the islanding process does not wrongly trigger the conventional underfrequency LS, and ROCOF islanding detection (RID) relays.

GRAPH SIMPLIFICATION
The proposed simplification algorithm arises from the principle that at least one path contains generator buses of each basic coherent group. Representing that path's elements with one single node leads to a reasonable simplification in the graph dimension. A set of elementary graph simplification operations applicable to our problem is demonstrated in [32]. These operations are implemented in the network shown in Figure 2 as a part of the IEEE 39-bus power system. Three equivalent generator buses 6, 10, and 19 for the reduced The final graph simplification step is to achieve an aggregated bus for each basic coherent group of generators. The proposed algorithm for this aggregation is based on the Steiner tree algorithm for a weighted undirected graph. The output of a Steiner tree algorithm for a graph is a sub-tree that contains a number of predetermined nodes [33]. The simplified graph's structure in this section is such that the proposed Steiner tree algorithm's output is the only path spanning generator buses of each basic coherent group. This path is identified using a simple search for those paths associated with the spanning tree and their corresponding values from the adjacency matrix. The total weight for branches of this path is minimal. The algorithm procedure is shown in Figure 3. This algorithm employs a minimum spanning tree algorithm as reported in [29]. In this procedure, all computations are according to the simplified weighted nondirected equivalent graph's adjacency matrix. Figure 4 shows the final simplified network for the system shown in Figure 2. In this example, all buses and branches involved in the path (P) are aggregated as an equivalent bus (bus 6 in Figure 4). As can be seen, after bus aggregation using the proposed algorithm, a different graph with a new adjacency matrix is achieved. Note that the final graph may contain parallel branches between some of the nodes, for example, those between buses 6 and 12 in Figure 4. In such a case, the branch with the least weight is adopted to build the new adjacency matrix. As mentioned before, this kind of simplification aggregates electrically closest buses to preserve the electrical characteristics of the original network as much as possible. In [13], a network reduction method on a node-weighted graph (active power of node as its weight) was proposed that uses the following four rules to retrieve a more manageable network to solve the islanding problem with reasonable complexity: (i) Removal of zero-weighted buses. (ii) Merging two buses by the accumulation of their weights.  This method was aimed to reduce the search space of the islanding problem by reducing the number of related variables. Consequently, it could sacrifice some proper splitting solutions. Moreover, by virtue of rule number 3, this method is not an offline task. It needs to be executed periodically.
Unlike the method in [13], our proposed simplification is focused on reducing the overall paths between incoherent areas. It has some useful features as follows: (i) It is an offline task. (ii) It ensures each island's generator buses' internal connectivity without adding connectivity constraints in the FA-FSconstrained MILP model. (iii) It reduces the number of power flow constraints and their related variables and disconnection constraints (reducing their related paths). (iv) Using coherency and electrical distance concepts, a sensible simplification to have a number of rotor-angle stable islands is accomplished.
As a comparative study, the method in [13] and the proposed method in this paper are used to simplify the IEEE 30-bus power system. The reduced network resulted from the method in [13] without zonal aggregation includes 23 buses and 26 lines. By a simple zonal aggregation, three more buses and five more lines can be removed. Considering tripping line 6 (between buses 3 and 4), the total number of buses will be 20, and the total number of lines will be 20 in the reduced network. Applying larger zonal aggregations can further narrow down the search space, but this certainly degrades the solution's reliability and quality. Also, geographical aggregation of buses cannot effectively reduce the number of paths between incoherent areas. The reduced network that resulted from our proposed method includes 26 buses and 32 lines. Our proposed method is aimed to reduce the number of paths related to disconnection constraints. Thus, this type of aggregation prepares a reduced network consistent with minimal paths for disconnection constraints. The IEEE 30-bus system has 1580 paths related to disconnection constraints. After simplification, the total number of paths will reduce to 144. By considering the same online condition (tripping line 6), the total number of paths will decrease to 18. Accordingly, an impressive reduction in the number of disconnection constraints is achieved (from the original 1580 constraints to the final 18 constraints after simplification). Also, our proposed method, as compared to the method in [13], results in fewer disconnection constraints (i.e. fewer paths) in the final FA-FS-constrained MILP model. The comparative results are stated in Table 1, and our achieved numbers of paths are underlined.
In summary, being offline, ensuring the internal connectivity of coherent areas, and fewer disconnection constraints are three significant advantages of our proposed method in this paper as compared to the method in [13]. Moreover, using the average error (AE) on power flow for preserved lines in the reduced network, the reduced network's power flow deviations, compared to the original network's referenced values, are evaluated. Three loading levels, including light, medium, and high conditions, are incorporated on the 79-bus case study, IEEE 39-bus, and IEEE 30-bus test system. The average error is defined as NL introduces the number of preserved lines. The results for three loading patterns are revealed in Table 2. As shown in Table 2, the AE is less than 6% in all cases. Accordingly, the reduced network based on our reduction technique is a reasonable alternative to the full network from a power flow deviation perspective.

GRAPH-PATH DETERMINATION
The following set of constraints needs to be considered in the optimal controlled islanding problem: (i) Disconnection constraints, which are essential to ensure the separation of incoherent areas. (ii) Bus-allocation constraints, which are essential to tackle the power imbalance for frequency control constraints.
These constraints are formulated based on their respective paths in the final simplified network. The disconnection constraints are formulated based on basic paths between aggregated buses, and bus-allocation constraints are formulated based on basic paths between every aggregated bus and non-generator buses. Figure 5 proposes a procedure for the graph-paths determination that leads to an offline database for basic paths. The proposed approach is based on Yen's k-shortest paths algorithm in [30]. Two main stages of the proposed procedure of step 4 are: (i) Determining all paths between aggregated buses (named as inter-area paths) as well as every aggregated bus and nongenerator buses (called bus-allocation paths) for the final simplified network. The first stage is carried out using a modified Yen's k-shortest paths algorithm. In this modification, we remove the counter and add a conditional break to the main loop. The second stage includes a sifting process on the output of the first stage. The sifting process comes up from the rule that if there is a path that is passed through the other aggregated buses, this path is not a basic one.
As an example, Figure 6 shows a 6-bus simplified power network with three aggregated buses (1, 2, and 4) and three nongenerator buses (3, 5, and 6). According to the proposed procedure, path data for an inter-area case (bus 1 as a source bus and bus 2 as a destination bus) and a bus-allocation case (bus 2 as a source bus and bus 3 as a destination bus) are found in  Table 3. As can be seen, the total number of paths for these cases is 7. By applying the proposed path determination procedure, only two paths (numbers 2 and 5) are determined as basic paths. Other paths are removed from the basic paths database because of passing through other aggregated buses denoted by the bold red font in Table 3.
In summary, the path-determination procedure acts on the simplified network to prepare a set of essential paths to establish disconnection and bus-allocation constraints. These paths are a series of adjacent lines. As mathematically explained in the next section, disconnection constraints are obtained directly organised as a series of linear relationships, while bus-allocation constraints are expressed as a combination of cross-product terms which are linearised using some useful rules.
In step 5, as mentioned before, the last filtering process is carried out according to the network's online status. In online conditions, two factors can reduce the number of basic paths. The first one is the number of paths between merged basic coherent groups, and the second one is the number of switched off transmission lines. Accordingly, this filtering process leads to determining a set of paths among basic paths characterised as final paths. These paths are utilised to establish the respective constraints.

THE FA-FS-CONSTRAINED MILP MODEL
A binary variable is assigned to every line and every generator for their switching status after the islanding. In the literature, minimal LS, minimal power flow disruption, and minimal power imbalance have been used as three different objective functions in the controlled islanding optimisation problem. Minimising each island's power imbalance leads to less FD, while the achieved islands' power flow experiences extreme variations. In such a case, severe power swings may happen, and a challenging post-islanding situation may arise. Minimising the power flow disruption concludes with much more power imbalances. In this case, corrective control actions should ensure the frequency stability of the islands. Our FA-FS-constrained MILP model directly incorporates frequency control constraints. Therefore, minimising the LS subject to FA and FS constraints plus coherency-based grouping improves the post-recovery process. Without a lack of generality, minimal load shedding is adopted here as the objective: The problem constraints are as follows.

Frequency-arresting and frequency-stabilising constraints
In an interconnected power network, islanded coherent areas' frequency responses vary depending on parameters such as inertia constant, time constant, and droop factor of each coherent area [34,35]. The equivalent single-machine swing equation of a coherent area states the coupling between the coherent area's frequency response and its power imbalance just after the islanding. The power systems with inverter-connected renewable generators have decreasing and time-varying inertia. Therefore, the frequency dynamic of these power systems differs from conventional power systems. As reported in [36], a WP generator cannot effectively intervene in frequency control activities. High penetration of low-inertia generators (such as WP generators) can affect the equivalent inertia and governor participation of an isolated power system. For example, the relationship between Nordic frequency response and WP penetration is studied in [37].
Following the study in [37], we assume a linear relationship between the penetration level of low-inertia generation P low m , ΔH m and ΔDR m in our final MILP model. This linear relationship is presented as where C I m and C D m are assumed -1.271 pu/s and -1.294 pu/Hz based on the study in [37], respectively.
Also, piecewise linear differential equations of FA and FS constraints based on Euler's discretisation method are implemented. The SGC variable is included as well. Moreover, the power imbalance of each island is defined through busallocation variables. The link between bus-allocation binary variables and line-switching binary variables is defined as a set of linear equality and inequality constraints. Frequency stability constraints are stated in Equations (5) to (15). By defining nΔt instead of t and using Euler's method for the equivalent single machine swing differential equation, three consecutive discrete equations are achieved as indicated in Equations (5) to (7). These equations construct a piecewise linear framework for FA and FS constraints. These equations' initial conditions are zero, given that there is no FD and governor action before distur-bance. As reflected in Equation (6), the last term of the right side of this constraint models SGC. This kind of formulation is more flexible and practicable than conventional formulation missing SGC. Adding SGC to the modelling results in more flexibility to obtain feasible solutions and less load shedding values in the optimal solution. It is worth noting that the lack of this variable, in some cases, even may lead to the problem infeasibility (which will be shown later).
Additionally, Equations (13) to (15) are bus-allocation constraints that link power imbalance variable after line switching to other related binary variables. In such a way, the important variable of power imbalance for each island is accurately involved in FA and FS constraints. The inertia and governor participation changes are addressed in Equations (6) and (7), respectively. After the coherency detection of WP stations in step 5 of the proposed method, each area's penetration level will be determined. Therefore, inertia and governor participation changes are known parameters in the proposed MILP model. The frequency bound and ROCOF constraints are reflected in Equations (9) and (10). In the following constraints, Δ f m,n represents the frequency deviation, K m,n introduces the rate of change of FD, ΔP mech m,n is the mechanical output power change, P Gpost m,n reflects the generators' output power of the area, and P Gpre e,m is the generators' output power of the area before the islanding: The r,h is bus-allocation binary variable from non-generator bus r to aggregated generator bus h. Also, Γ hr,o reflects multiplication of all switching status binary variables associated with lines forming path o between buses h and r. These cross-product terms can be expressed as a set of linear inequality constraints [38]. N Path h,r is the number of paths between buses h and r. As stated in Equation (15), Y is a linear function that represents a set of linear inequalities and includes Γ hr,o , Z hr,o , and U hr,o . Also, Z hr,o is the set of auxiliary binary variables for linearisation, and U hr,o is the set of switching status binary variables of lines constructing the respective path.
Constraints of Equation (14) ensure that if at least one path between an aggregated bus and a non-generator bus is fully connected (all lines are switched on), then the variable of bus-allocation between these two buses is 1 (otherwise busallocation variable is 0). In this case, the sum of all Γ hr,o related to all paths divided into N Path h,r is forced to be a positive real number (less than 1). Consequently, the relevant bus-allocation binary variable is forced to be 1. The conventional way to linearise a cross-product term of binary variables is the one that is indicated in (a) and (b) in the Appendix. This way is replaced by some useful rules as discussed in [38]. Two of these rules, compared to the conventional way, are stated in the Appendix. Therefore, the sum of several cross-product terms is linearised by fewer constraints than conventional formulation.

Power flow constraints
Power flow constraints are formulated according to the DC load flow model, Equations (16) to (21) (it is assumed that voltage regulation of network is carried out by local compensation of reactive power after the islanding). P F j is power flow of line j, P shed i is load shedding of bus i, u j is the switching status of line j, and P D i is the value of load in bus i. Power balance constraints at each non-generator bus are included as well. Additionally, the maximum rating of in-service lines is incorporated. The aggregated bus of each island is defined as a slack bus with a phase angle of zero:

Parallel lines constraints
As previously stated, in the graph simplification process, there is a possibility to have parallel lines. At least one terminal of each parallel line is connected to an aggregated bus in the network's final graph. Technically, the switching status of all parallel lines between every pair of buses should be identical as stated in Equation (22).

Disconnection constraints
Disconnection constraints are based on the final paths derived in step 5 of the proposed approach. Disconnection constraints force the solution to have the number of islands equal to the number of final coherent groups. Of course, rearranging the disconnection constraints can change the overall number of islands. But, the most reasonable islanding strategy in terms of rotor angle stability is the one that is established based on the coherency concept. These constraints are derived from the fact that at least one transmission line of each path is switched off to isolate all incoherent areas completely. These constraints are simply formulated as inequality constraints implying that the sum of binary variables corresponding to lines constructing a given path is less than or equal to path length minus 1 as formulated in Equation (23). The path length is the total number of transmission lines forming the path. This simple formulation is a graph-based splitting strategy that ensures coherency-based isolation for the network. These constraints are formulated in the simplest possible linear form:

Other constraints
The total number of switching operations for controlled islanding can be limited due to operational and computational reasons as is indicated in Equation (24): We should note that power system restoration constraints like black-start capability, sufficient generation capacity, and unit commitment preservations can also be added to the proposed MILP model [23].

SIMULATION RESULTS
The IEEE 39-bus system, including 10 generators, 34 transmission lines, and 12 transformers, is studied. A 79-bus case study based on real data is also analysed to prove the proposed approach's effectiveness. The FA-FS-constrained MILP model is solved using convex (CVX) in the MATLAB environment. Time-domain simulations are carried out using DIgSILENT PowerFactory software to verify our results. All simulations are executed on a PC with Intel Core™ i5 CPU @2.60 GHz and 8 GB RAM. The commercial MILP solvers have remarkable capabilities in tackling properly modelled optimisation problems. Moreover, these solvers are very stable and can easily find a feasible solution to the optimisation problem and even a near-optimal one in a limited period. Here, the developed MILP model was tackled through the MOSEK solver. This solver uses the branch and cut algorithm to converge to the final-optimal/near-optimal solution [39]. Also, it utilises some speeding-up methods such as applying good initial conditions and relaxing the termination criterion. Considering an MIP gap of 0.2%, the computation time to solve the proposed MILP model for different case studies of simulation networks was less than 1 s, which is compatible with the WAMPAC requirements. Our proposed methodology incorporates four offline steps to reinforce online actions in terms of technical and computational efficiency. The computational time of the offline steps for the IEEE 39-bus case study and 76-bus case study are 9.2 and 13.7 min, respectively. For low-inertia studies, coefficients of linear relationships indicated in Equations (3) and (4) are approximated based on time-domain simulation data from the IEEE 39-bus system. Accordingly, coefficients of Equations (3) and (4) are achieved as C I m = -1.3339 pu/s and C D m = -1.3843 pu/Hz.

IEEE 39-bus system: Scenario I
In this scenario, a three-phase short circuit occurs at the line between buses 26 and 29. Due to the delay in distance relay operation, the short-circuit condition remains for 160 ms. Consequently, the network is pushed towards instability. Using the offline method described in [26], the basic coherent groups of generators, including entire generations for both base and lowinertia cases are stated in Table 4. In the low-inertia case, four WP stations are located in buses 3,11, 15, and 21. Using the online method referred in [31] and according to loading, topological, and fault condition of the network before the islanding, there is no merged group for basic coherent groups of conventional generators in this scenario. Besides, employing this online coherency detection method, WP stations in buses 3, 11, 15, and 21 have belonged to coherent groups of 1, 2, 3, and 4. As a result, the final coherent groups of generators are similar to The islands in scenario I of IEEE 39-bus system for base case those expressed in Table 4 except that each of the four first groups has an extra WP station. In this scenario, the frequency control modelling of islands is performed in a time interval of 10 s with 0.1 s as the time step for piecewise linear modelling. The maximum allowable positive and negative FDs in all islands are considered as ±0.75 Hz to prevent interfering with under-frequency LS relays operation. The minimum value for ROCOFs is assumed -0.5 Hz/s to keep the RID relays inactive. The SGC is included in the frequency response modelling. Running the optimisation model for the base and low-inertia cases, six different islands are specified for the base case shown in Figure 7. In the low-inertia case, boundaries are such that buses 3, 17, and 18 are on island 1, and buses 4, 5, 8, and 14 are on island 6. Other than these changes, both cases are the same regarding territories.
The frequency response for each island is illustrated in Figures 8(a) and (b), respectively. The base case network's total LS is zero, while this value for the low-inertia case is 105.8 MW. Total SGC for the base and low-inertia cases is 1100 MW (G07 and G10) and 1037.1 MW (G02 and G07). All islands, except The frequency response for islands of IEEE 39-bus: (a) scenario I-base case, (b) scenario I-low-inertia case, (c) scenario II-base case, (d) scenario II-low-inertia case island 1, experience a decrease in frequency within the allowable range. Note that each island can capture a sudden power mismatch without intolerable FD and ROCOF. However, over a more extended period, the frequency control system will alleviate these deviations. The comparative results for the base and low-inertia cases of this scenario are shown in Table 5. It can be observed that there are more LS for the low-inertia case, compared to the base case, as a result of the reduction in inertia and governor participation of the areas with a non-zero penetration level of WP stations. As shown in Table 5, the low-inertia case's optimal solution contains two islands with minimum permitted ROCOF (islands 3 and 4). Boundaries of the islands are such that the optimal solution in low-inertia case, in general, includes larger values of FD and ROCOF (in terms of absolute values).
The optimiser termination time for the base case is 0.43 s since this time for the low-inertia case is 0.45 s.
It seems that adding low-inertia issues to the model concludes the islands with more LS and SGC typically. Thus, modern power systems containing several low-inertia resources can be separated into the feasible islands properly only if low-inertia challenges are incorporated in the model.

IEEE 39-bus system: Scenario II
In this scenario, two disturbances occur after each other. First, a three-phase short circuit occurs on the transmission line between buses 29 and 26 near bus 29. Distance relays isolate it after 120 ms. Followed by this event, the line between buses   1 and 39 is switched off because of the overload protection system's function. The controlled islanding is subsequently called to save the system from a blackout. In this scenario, the final coherent groups of generators are reported in Table 6. Three basic coherent groups, including groups 2, 3, and 4, merge and create a larger equivalent coherent group. Also, the WP stations are the same as those in the previous scenario for the low-inertia case and their coherency situation. An important observation is that there is no feasible solution for this scenario without the SGC variable. According to Equations (5) to (7), a large power imbalance (overgeneration) without SGC results in extremely high overfrequency.
The total number of paths between aggregated buses of basic coherent groups in this scenario is 366, including 132 basic paths. Basic paths between merged basic coherent groups and those including a switched-off transmission line would be removed. Consequently, the total number of final paths corresponding to disconnection constraints is 22. Actually, through graph simplification, bus aggregation, and filtering process in paths determination procedure, the number of final paths gets much less than the original paths, assuring the method's outstanding computational efficiency.
In this scenario, the frequency control factors are such that FD cannot remain in the acceptable range. However, SGC can neutralise the power imbalance. With our proposed method, four controlled islands for the base case are specified as shown in Figure 9. The only difference of boundaries between the base and low-inertia cases is the allocation of bus 7. This bus is placed   Table 7 shows comparative results for the optimal solution of base and low-inertia cases of scenario II. Like the previous scenario, the total LS in the lowinertia case is more than that in the base case. In the base case solution, without ROCOF constraints, islands 1, 3, and 4 violate the minimum permitted ROCOF. However, in the low-inertia solution, the ROCOF is restored in permissible range by more LS and this, in turn, changes the boundary for island 4. It can be taken from these results that minimising the total value of LS needs three control actions of transmission switching, LS, and SGC in a coordinated manner. The optimiser termination time for the base case is 0.54 s since this time for the low-inertia case is 0.59 s.
Technically, there is no significant direct relationship between the number of islands and their frequency response. Instead, the factors that directly affect the islands' frequency response are the value of power imbalance, inertia response, governor action, corrective control actions (LS and SGC), and the frequency control system's time constant. All these factors are directly incorporated into the model. Of course, there is more possibility to have larger power imbalances in the more extensive networks with fewer coherent areas, leading to more FD. However, frequency constraints ensure the frequency stability of the islanding solutions.

A 76-bus case study based on real data
A 76-bus case study is evaluated here to assess the performance of our proposed approach. The system serves 2364.2 MW load and has 400, 132, and 63 kV voltage levels. The generation system includes 38 generators in seven power plants (PPs). The transmission system consists of 107 lines. In the low-inertia case, three WP stations are placed on buses 19, 46, and 54. It should be pointed out that the base case's relevant data is originated from a real power system, but the authors assume lowinertia details. It is assumed that a three-phase short circuit on line 23 near bus 8 has occurred. The faulted transmission line is isolated in 0.13 s by the protection system. Afterwards, line 42 experiences a false trip in 0.99 s due to the third zone of distance relay's malfunction. This disturbance provokes the network generators to approach the out-of-step condition as well as transfers the power flow to nearby lines. As shown in Figure 10, some critical generators' out-of-step condition leads to isolated areas with large power imbalance values. Consequently, rotor angle instability is anticipated to launch the controlled islanding process. Five basic coherent groups of generators with total generation for both base and low-inertia cases are given in Table 8. Just before the islanding process, the system has four final coherent groups by merging basic coherent groups of 2 and 3. Using online coherency detection of step 5, WP stations connected to buses 19, 46, and 54 have belonged to coherent groups of 4, 5, and 2. The specified islands for the base case and associated frequency responses for both cases are shown in Figures 11,12(a), and (b), respectively. In the base case's optimal solution, the total LS value of the network is zero, and the total SGC value is 230 MW from PP03 on island 2. Island 2 is a multi-

FIGURE 11
The islands of real 76-bus power network

FIGURE 12
The frequency response of islands in 76-bus power system: (a) base case, (b) low-inertia case ple PP island, while island 3 is a single PP island. The power mismatches of both islands are positive, that is, over-generation status. However, the boundaries of island 3 ensure that its overfrequency is technically tolerable, while island 2 exploits SGC to maintain its frequency in the permissible range. Accordingly, the frequency of island 2 declines but does not trigger the automatic load shedding relays. As shown in Table 9, the total LS of the low-inertia case is 19.87 MW, which is more than the base case, much like the two previous scenarios. In the optimal solution of the low-inertia case, island 3 includes buses 50, 51, 68, 69, 70, 71, and 72, which belong to island 2 in the optimal base case solution. Also, island 1 includes bus 14, which is placed on island 2 in the base case. Other than these changes, all islanding boundaries are the same in both cases. Islands 2 and 4 have the highest level of penetration for low-inertia generations. In the low-inertia case's optimal solution, the ROCOF is preserved in the permissible range by more LS (island 4) and boundary change (island 2). Accordingly, low-inertia issues like inertia and governor participation changes, as well as ROCOF constraints, make our model more complicated and decrease the quality of the optimal solution, compared to the model without low-inertia issues. Thus, in today's power systems with high penetration of low-inertia resources, including these issues seems necessary.
The optimiser termination time for the base case is 0.54 s since this time for the low-inertia case is 0.83 s. Some additional operations can reduce the termination time of optimiser to be more applicable to WAMPAC systems as follows: (i) The predetermination of switching status of some transmission lines as well as some vital generators (ii) Pre-allocation of some non-generator buses to the coherent areas according to network operators' operational experiences and other relevant technical limits. (iii) More network reduction based on additional bus aggregation operations. (iv) Increase the time step to such an extent that the quality of the solution remains acceptable.

Verification with time-domain simulations
In this section, time-domain simulation results for the 76-bus case study's controlled islanding are obtained using PowerFactory software. Note that the simulations are with no simplification in component models. Figure 13 compares time-domain simulations (solid lines) with those captured in the optimisation model (dotted lines). It is observed that the approximations made in the proposed FA-FS-constrained MILP do not introduce significant errors.

CONCLUSION
This paper proposes a comprehensive MILP model for optimal controlled islanding. The FA and FS constraints are linearised and incorporated in our MILP model. Accordingly, our proposed FA-FS-constrained MILP model can effectively handle the optimal controlled islanding in the low-inertia power systems. To improve the computational burden of our MIP