New Cycle-based Formulation, Cost Function, and Heuristics for DC OPF Based Controlled Islanding

This paper presents a new formulation for intentional controlled islanding (ICI) of power transmission grids based on mixed-integer linear programming (MILP) DC optimal power flow (OPF) model. We highlight several deficiencies of the most well-known formulation for this problem and propose new enhancements for their improvement. In particular, we propose a new alternative optimization objective that may be more suitable for ICI than the minimization of load shedding, a new set of island connectivity constraints, and a new set of constraints for DC OPF with switching, and a new MILP heuristic to find initial feasible solutions for ICI. It is shown that the proposed improvements help to reduce the final optimality gaps as compared to the benchmark model on several test instances.


I. INTRODUCTION
Intentional controlled islanding (ICI) is an important system integrity protection scheme (SIPS) aiming to prevent system collapse due to wide-area instability by separating the system into a set of non-interacting islands. A typical situation requiring such a control action is loss of synchronism between some parts of the system, but it can also be used to limit the spread of cascading outages, or to isolate unstable parts of the system. The high relevance of this control action is confirmed by the increasing number of blackouts around the globe [1].
From the mathematical perspective, computing the splitting boundary to separate the system into a number of selfsustained islands equates to an OPF problem with switching. Such problems are notoriously hard to solve, as OPF is nonlinear, and switching decisions are discrete, which amounts to a mixed-integer nonlinear program (MINLP). Moreover, ICI must be solved rapidly in real-time as a corrective control action for instability mitigation. For this reason, solving ICI as a MINLP is not feasible, and MILP reformulations of the original MINLP are often seen as an acceptable trade-off between accuracy and computation time.
Although multiple MILP approximations for ICI and optimal transmission switching (OTS) exist in literature (e.g., [2], [3]), this paper is focused on the basic DC OPF ICI model.
Our goal is to resolve the fundamental issues associated with the DC OPF ICI model to devise the efficient computational enhancements that could also be applied to more advanced MILP-based ICI models [2]- [4].
We first observe that the conventional ICI objective of minimizing the total post-splitting load shedding [5], [6] may often result in a poor MIP optimality gap progression if the exact optimal load shedding is close to zero. To alleviate this issue, we propose to minimize the total load-generation imbalance in the computed islands while possibly keeping the load shedding as the secondary optimization objective. Minimizing the total load-generation imbalance should limit load shedding and additionally limit the value of the rate of change of frequency (ROCOF) following the system splitting.
Our second contribution aims to improve the modeling of the islands' connectedness requirement. The previously published papers on optimization-based ICI either use artificial commodity flows to impose connectivity constraints on islands [6] or neglect this requirement [2], [7]. Our proposal is to model connectivity through directed spanning trees, which aims to avoid the computationally inefficient big-M coefficients that accompany the commodity flow based formulation. This technique has previously been applied to enforce radiality in distibution networks [8], but our approach differs from [8].
Our third contribution addresses the big-M coefficients that are normally present in the DC OPF constraints related to switching decisions. We reformulate these constraints as Kirchhoff's voltage law (KVL) equations with additional slack variables to eliminate the uncertain voltage angle differences across an open line otherwise modeled as large big-M constants. Cycle inequalities involving KVL have been previously proposed in [9], but our work differs from [9] in several aspects. The goal of [9] was to strengthen the LP relaxation of the conventional DC optimal transmission switching (OTS) model, while our goal is the elimination of the big-M coefficients. This motivates a different modeling approach.
Our fourth contribution is a new MILP heuristic to find initial feasible solutions for DC OPF ICI. For some MIP problems, a feasible solution can be obtained trivially or in polynomial time, but the OPF ICI problem combining power flow physics and discrete decisions does not allow for a simple initialization. Obtaining an early feasible solution is crucial for the progress of a modern MIP solver, as there are efficient heuristics that are able to progressively improve the best feasible solution (e.g., the RINS heuristic). Unfortunately, we could not identify any published MILP heuristics for OPF ICI, which prompted the design of a simple yet effective method based on the LP relaxation of DC OPF ICI.
The following sections describe the above contributions in more detail. Section II-B introduces the used notation and the conventional benchmark DC OPF ICI model. Section III details our proposed model and cost function. Section IV introduces the MILP heuristic for OPF ICI and explains how MILP is solved. Section V describes the test cases and presents the computational results. Section VI concludes the paper.

A. Notations
A power network consisting of n nodes and m branches is represented by a graph with set of nodes V, set of edges E (|E| = m), and set of arcs A (|A| = 2m). Power demand and power generation at node i prior to splitting are denoted as P s L,i and P s G,i . Load demand and generation shedding to satisfy the post-splitting conditions are denoted as P L,i and P GS,i . Power flow through branch (i, j) ∈ E is defined as p i,j , its value prior to splitting is p s i,j , and its limiting magnitude is p max i,j . Admittance of branch (i, j) is denoted as b i,j . As the used model is DC OPF, reactive power is not modeled.
The splitting configuration is modeled by binary variables x i,k (x i,k = 1 implies that node i is assigned to island k) and y i,j (y i,j = 1 implies that branch (i, j) ∈ E is open). The total number islands is K. An important ICI requirement that promotes transient stability after splitting is to assign generators from the same coherent group to one island. Coherent groups are denoted as G 1 , . . . , G K (i.e., each coherent group is supposed to form a separate island).
Indices i and j run over network nodes and pairs of nodes (edges or arcs). Index k runs over islands (k = 1, . . . , K). Superscripts s, * , LR denote pre-splitting values, integral solutions of MILP programs, and linear relaxation (LR) solutions of MILP programs respectively.

B. Basic MILP Model for DC OPF ICI
Probably the first paper proposing a MILP-based ICI strategy was [5]. Since then several related models appeared in the literature, of which the model in [7] is the most compact one. Unlike [5], [10], [11], the model in [7] does not use auxiliary binary variables that model the membership of edge (i, j) in island k. Instead, the edge status variables y i,j are used directly, which eliminates O(mK) extra binary variables and constraints. However unlike most other references, the models in [2], [7] neglect the islands' connectedness requirement. To avoid the deficiencies of any single existing model, a combined DC OPF ICI model is introduced below: The objective of (1) combines the relevant objectives from [5], [10], [11] with different weight factors. The factor β is associated with minimal load shedding, which is of great importance for ICI. While minimizing load shedding, it is important to avoid arbitrary large values of shed or disconnected generation, which is achieved by introducing the minimization of generation shedding into (1a) with the weight γ. Additionally, [10], [11] use the minimization of total power flow disruption as their cost function. In our experience, this choice of cost function allows for an easier convergence to the optimum, but it cannot well enough represent the power balance within islands, which is the main challenge of ICI. However, a large power flow disruption is undesirable [12], which explains its presence in (1a) with the weight µ.
Switching constraints (1b)-(1f) ensure the separation of buses belonging to different islands from each other through open branches (y i,j = 1). Here (1f) ensures that lines connecting buses in the same island are closed, thus precluding line switching inside islands. Because of (1f), the formulation in (1) becomes easier to solve, as its feasible region becomes more constrained. From an operational perspective, requiring both the optimal splitting cutset and optimal line switching inside of each island may be too complex and impractical during an emergency condition.
DC OPF constraints (1g)-(1m) describe the physics of DC power flow. Kirchhoff current law (KCL) is modeled by (1g), which represents the power balance at each node. Constraints (1j)-(1m) model Ohm's law under the presence of switching decisions y i,j , which requires the introduction of the node potential variables ϕ i , ∀i ∈ V. When branch (i, j) is closed, p i,j is governed by Ohm's law. If branch (i, j) is open, p i,j = 0 and the angle difference (ϕ i − ϕ j ) can be arbitrary. Therefore, the big-M constants M ϕ i,j are needed to allow (ϕ i −ϕ j ) to take some rather large values. Additionally, the lower and upper bounds ϕ min and ϕ max in (1n) can also be seen as big-M constants, as their values often cannot be estimated exactly.
Single commodity flow constraints are used to ensure the connectedness of each island in the same way as the switching constraints (1b)-(1e) ensure the separation of islands from one another. For each coherent group G k , a root node r k ∈ G k , r k ∈ R is selected, whereby R is the union of root nodes of all groups. Constraint (1r) requires that the number of units of an artificial commodity produced at r k equals the total number of nodes in island k corresponding to G k , while constraints (1s) demand one unit of commodity to be consumed at each non-root node. The satisfaction of flow balance (1r)-(1s) requires each island to be connected. However, the upper bounds of commodity flows for each branch are not known, which prompts the usage of another big-M constant n − 1 in (1p)-(1q). This big-M constant could be somewhat lowered through exact calculations, but n − 1 is shown here to simplify the representation. Unlike [5], [10], (1) does not require a separate set of commodity flow variables for each island, as it appears that a single set of m commodity flows f ij can be used for all islands.

A. DC OPF ICI using Cycle-based KVL Constraints
A major issue with (1) is the presence of big-M constants M ϕ i,j in (1l)-(1m). While there are methods to bound these coefficients when the solution is required to contain a single connected component (e.g., [13]), no bound strengthening method is known for the disconnected networks resulting from ICI. Thus, a possibly unrealistically large value needs to be assumed for M ϕ i,j in order to preserve all feasible solutions. An alternative to Ohm's law form of power flow equations in (1l)-(1m) is to apply KVL to the entire network, while considering the switching decisions. The KVL equations are written with respect to a set of linearly independent cycles in the network, which is known as cycle basis. For the reasons that become clear below, it is desirable to compute the cycle basis as quickly as possible, which is achieved by using as cycle basis fundamental cycles induced by the network's minimal spanning tree (MST) [14]. Given cycle basis CB, the KVL alternatives of (1l)-(1m) can be formulated as follows: where C ∈ CB is a cycle formed by undirected edges (i, j) ∈ E, and function sgn(i, j) takes the plus sign if edge (i, j) is aligned with the clockwise direction of C and the minus sign otherwise.
As (1f) precludes line switching inside islands, each cycle is either connected or has at least two open lines. Therefore, i,j , M C can be computed exactly, but it may still be quite large for long cycles. Thus, to additionally strengthen the model, we are adding some additional KVL constraints (2) associated with cycles shorter than a certain length. All such cycles can be found in polynomial time by applying to each node techniques based on recursive node neighborhood traversal. Short cycles (e.g., up to length 7) can be found quickly before the start of the MIP solver, but enumerating longer cycles becomes increasingly computationally inefficient.
Obviously, any cycle basis of the disconnected network contains less cycles than the cycle basis of the original network. As only a small fraction of edges is opened for ICI, the majority of fundamental cycles of the disconnected network coincide with the initial fundamental cycles. However, in some cases the original fundamental cycles do not fully describe the cycle basis of the disconnected network. Thus, for each integer solution obtained during the solution process, cycles violating (2) need to be identified, and the corresponding inequalities (2) added to the model by using MIP solver callbacks. A possible computationally efficient way to achieve this is as follows: 1) Compute the MST of the original network and the associated fundamental cycle basis CB 0 . 2) Assign close-to-zero weights to the network edges belonging to the MST. 3) Once an integer solution is found, compute its minimal spanning forest and the associated fundamental cycle basis CB i (i is the solution number). The previously assigned small edge weights should promote the alignment of the newly computed cycle basis with the original one. 4) Check the solution for conformity with (2) (i.e., bi,j sgn(i, j) = 0, ∀C ∈ CB i ). 5) If any cycles in CB i violate (2), use them to add new KVL constraints (2) to the model. Besides strengthening the model, the abovementioned KVL constraints based on extra short cycles also improve the satisfaction of (2) for the intermediate integer feasible solutions.

B. Island Connectivity Constraints
The single commodity flow constraints that are most commonly used to enforce islands' contiguity [5], [6], [10] have the drawback of introducing at least m auxiliary variables and large big-M constants that may lead to loose LP relaxations. In [11], two alternative methods based on multicommodity flows and network cutsets are proposed. However, these alternatives may often underperform the simple single commodity flow approach. At least the multicommodity flow method has largely theoretical importance due to the very large number of auxiliary variables that it requires [15]. The cutset-based method requires solving many max-flow min-cut problems to identify cutsets that violate connectivity, which has the time complexity of O(mn 2 ) for a single max-flow min-cut run.
In this section, a different approach to islands' connectivity is proposed that does not involve big-M coefficients and allows to quickly find strong inequalities that ensure network connectivity. This approach is related to the studies on radiality constraints in distribution networks (e.g., [8], [16]), but has some application-specific features and computational improvements. One of conditions for network connectedness is that there exists a spanning tree with n − 1 edges that spans all of its nodes. In general, for a network with k connected components there exists k spanning trees (one for each component) forming a spanning forest that includes n−k edges. This observation motivates the following connectivity formulation based on rooted directed spanning trees: where z i,j and z j,i represent the status of arcs (i, j) and (j, i) ∈ A (z i,j = 1 iff (i, j) is enabled).
The equalities in (3b)-(3c) define the root nodes of K coherent generator groups as root nodes of K trees each spanning one island. According to (3b)-(3c), each root node should have no incoming arcs and at least one outgoing arc. Conversely, each non-root node should have one incoming arc according to (3d). Constraints (3e)-(3f) link z and y variables. For the edges that participate in network cycles (i.e., in the cycle basis described in Section III-A), (3e) is valid. For the rest of the edges, the stronger version (3f) can be used instead of (3e). The number of enabled arcs is constrained by (3a).
Combined constraints (3a)-(3f) can define k connected components spanned by k directed trees, or more than k connected components. In the latter case, the additional components will form directed cycles to satisfy (3a)-(3f). Therefore, cycle breaking constraints (3g) are needed to only allow spanning forests as solutions, where D is the set of all directed cycles in the network and C is a sequence of arcs in A forming a directed cycle. As there are exponentially many cycles in the network, all constraints in (3g) cannot be simultaneously added to the model. Instead, only the violated inequalities (3g) are added in MIP solver callbacks that are available in most of modern MIP solvers (see Section III-E).

D. Minimization of Load-Generation Imbalance
Minimal load shedding is often chosen as the cost function of ICI, as it is highly important for utilities to minimize the loss of load. However, shedding an exact amount of load at exact bus at exact moment in time to realize the outcome of (1) may be sometimes quite challenging with the available equipment. Moreover, the optimal load shedding can be zero, in which case it will be impossible to provide a non-zero lower bound, and the integrality gap of the MIP solver will be close to 100 % during the whole optimization.
Among practitioners, the notion of balanced islands with good load-generation balance is no less popular than the idea of minimal load shedding. Islands with good power balance would naturally limit load shedding without the need of picking specific buses, while also limiting the ROCOF value after splitting. At the same time, the total power imbalance of all islands is always non-zero and can be much easier bounded from below. Based on (1), the minimization of load generation imbalance can be formulated as follows: where P ∆,k represents the power imbalance of island k. The inequalities in (5b)-(5c) are constraining the absolute value | i∈V (P s G,i − P s L,i )x i,k | from below for each k, while the minimization requirement (5a) constrains it from above.

E. Overall Solution Process
The ideas from Sections III-B-III-A can be combined together into the following novel DC OPF ICI formulation: (3a) − (3h), (4), (5b) − (5c) The above formulation requires MILP callbacks for cycle breaking constraints (3g) to be solved efficiently. In our implementation, the violated cycle breaking constraints are easily separated using Algorithm 1 for each new integer solution.
Given an integer solution z * i,j , it can be considered as directed graph G. For this graph, the spanning forest F can be C ← P ∪ a i,j 10: C ← C ∪ C 11: end for 12: return C computed by using the standard spanning tree algorithms with loglinear time complexity. If F has k connected components, the solution is accepted, and the best integer solution gets updated. Otherwise, Algorithm 1 finds arcs that are in G, but not in F . Each such arc corresponds to a directed cycle, which allows to retrieve the full set of cycles C that violate (3g) in loglinear time.
In addition to cycle breaking constraints (3g), we are adding directed cutset constraints between the root node of each group and the remaining generators of the same group using the same z i,j variables. These constraints are convenient for enforcing connectivity of a defined set of terminal nodes to the root node. A detailed description of these constraints including improvements like back cuts and creep flow can be found in [15] and several other references.
The cycle breaking constraints in (3g) also combine well with the cycle-based DC OPF ICI formulation in (2). To apply (2), a cycle basis and optionally some additional network cycles need to be found. From these precomputed cycles, the shorter ones (e.g., with the length less than 6-8 edges) can be used to strengthen the initial model formulation with some constraints of type (3g).

IV. MILP HEURISTICS
When solving (6), no feasible solution could be found within the optimization time limit for several larger instances. In practice, it is highly beneficial to obtain an initial feasible solution as early as possible because this would increase the available time to improve the upper bound of (6) using various efficient techniques available in the modern MIP solvers. As the practice has shown that multiple MILP heuristics available in Gurobi [18] often cannot retrieve a feasible solution, a new problem-specific heuristic had to be proposed to improve the solvability of (1) and (6). This heuristic is outlined in Algorithm 2.
The LR of x i,k is taken as input and analyzed in lines 2-7 to identify the edges that are considered to be closed by the current LR. Then these edges are used to build an undirected graph G in line 8, and the connected components of that graph are identified. Next, each connected component is analyzed for the presence of a root node r ∈ R in it. If a root node is present end for 7: end for 8: CC ← connected components(G(V, E on )) 9: P ← ∅ // Set of partial islands 10: for CC ∈ CC do 11: r ← R ∩ CC 12: if |r| = 0 then continue 13: else if |r| ≥ 2 then return ∅ Fix x i,k , y i,j , z i,j according to P and solve the reduced MILP (e.g., (1)) with more than 80% fixed buses. 19 is the component, it is saved to the set of partial islands P. Finally, the x i,k variables of the nodes belonging to contiguous partial islands are fixed, and the values of the remaining node variables are obtained by solving the initial problem as MILP with the majority of discrete variables being fixed. Because of a high degree of reduction of discrete variables, the reduced MILP is usually solved in a fraction of a second. However, it is important to find connected partial islands, as just rounding the LR will very often result in infeasibility. The heuristic in Algorithm 2 is only run until a feasible solution becomes available, as modern MILP solvers are able to apply more efficient heuristics to consistently improve the existing solution.

A. Test Setup
This section 1 presents a comparison of the proposed formulation in (6) to the known benchmark model in (1). The results are produced for a number of test networks from the MATPOWER toolbox [19]. As publicly available power system test cases do not include coherent generator groupings that are required in (1c), those were obtained with the generator coherency algorithm from [20] 10. The reported solution times include the presolve and MILP solver time, but exclude the model preparation time, which is polynomial. The time limit for the test cases with less than 500 buses has been set to 480 seconds and to 720 seconds for the test cases with more than 500 buses. In all the experiments, the default Gurobi solver settings were used. The maximal time to be spent in the MILP heuristic from Section IV is set to 3 % of the time limit. For simplicity, the auxiliary MILP model in Algorithm 2 is based on (1) for all test cases. The solution was considered optimal if its MIP optimality gap was less than 1%.

B. Minimization of Load Shedding
The first set of test results is related to minimizing the known objective in (1a). The objective weights are set as follows: α = 0, β = 1, γ = 0.01, µ = 0.1. The relatively high value of µ is selected mainly to ensure non-zero lower bounds in all cases and thus a more informative progression of the MIP optimality gap during the optimization. The big-M coefficients M ϕ i,j in (6) are set to 2π, ϕ min is set to −π, ϕ max is set to π , and the power flow limitation p max ij is assumed to be equal to b i,j π/4, ∀(i.j) ∈ E.
First, the results for the baseline formulation in (1) are presented in Table I. The table header can be read as follows: U B is the best MIP objective value (i.e., the upper bound), g is the final optimality gap, T is the solution time, P LS , P ∆ , P GS , p Σ ij are the final values of load shedding, total power imbalance, generator shedding, and power flow cut respectively. In Table I, g and T are included interchangeably:  if g is less than the optimality tolerance then T is included, otherwise g is shown, as T equals to the time limit. Next, in the formulation in (1) the flow-based connectivity constraints (1p)-(1s) are replaced with (3) and the MILP heuristic is enabled. The results obtained after this modification are summarized in Table II. Finally, the results of our complete formulation in (6) with the MILP heuristic enabled for the given cost function are given in Table III. As it can be seen, the results in Tables II-III in terms of  the MIP optimality gap and solution time usually exceed those  in Table I, which demonstrates the efficacy of the proposed improvements.

C. Minimization of Total Power Imbalance
To complement Section V-B, the same test cases are recomputed with the main objective of power imbalance minimization. The objective weights are set as follows: α = 1, β = 0.01, γ = 0.01, µ = 0.01. Now the highest priority is given to generation-load imbalance minimization, with load shedding being included with a small weight. As the initial power imbalance is nearly always substantially higher than zero after the network is split, non-zero lower bounds can be produced without assigning an increased value to µ.
The results comparison follows the same route as in Section V-B. The results of minimizing the new objective with the formulation in 1 are listed in Table IV. The outcome of the exchange of (1p)-(1s) with (3) is illustrated in Table V. Finally, the results of our formulation in (6) are shown in Table VI.
As it can be seen from Table IV, with the change of the objective function more cases fail to compute at least a feasible solution with the basic formulation (1). Thus, the role of feasibility MILP heuristics becomes larger. In fact, the proposed MILP heuristic is enabled for the results in Tables V-VI, and feasibility is achieved for all the test cases. The performance is also improved with the proposed enhancements, which can be seen by comparing the 4 th columns of Tables IV-VI. An additional important observation can be made by comparing the upper bounds of case89pegase across Tables IV-VI. Although the optimal solution could be reached in all cases, the optimal value is better with our complete formulation (6). This is explained by the absence of big-M constraints in (6), the values for which must be selected. Thus, (6) can be certain to achieve the true optimal solution. In the present case study, it is possible to achieve the same optimum using (1) if the angular big-M constraints are substantially increased (e.g., by a factor of 10). 2

VI. CONCLUSIONS
This paper has proposed several computational enhancements to the DC OPF ICI problem. The main idea was to eliminate the big-M constants that are present in the most existing formulations in order to tighten the linear relaxation of DC OPF ICI. The big-M constants present in island connectivity constraints have been removed by replacing the popular single commodity flow based connectivity model with the new model based on directed spanning forests, which does not require any large coefficients. The big-M constants associated with ICI switching decisions and DC OPF in general have been excluded by replacing the DC OPF switching constraints based on Ohm's law with the new ones based on KVL.
Besides of handling of all the three types of big-M constants in the standard DC OPF ICI formulation, it has also been observed that the existing ICI objective functions may not fully correspond to the practical requirements associated with system splitting. To ameliorate this, the new objective function minimizing the load-generation imbalance after splitting has been introduced. With this objective, situations when the MIP solver cannot find an initial feasible solution have become more common. To fix this, a new MILP heuristic for DC OPF ICI has been proposed as well. In general, the solution process was noticeably influenced by the type of objective function. Thus, it could be possible to choose "convenient" objective functions to quickly find feasible initial solutions.
Finally, the proposed enhancements have been tested against the compact and efficient baseline DC ICI OPF model from the literature, and it has been found that the novel ICI model produces better results in most of cases. In the future, it is planned to apply the findings of this paper to more complex and accurate ICI models and possibly to other transmission switching problems.   (6)