# A Hybrid Submodular Optimization Approach to Controlled Islanding with Post-
Disturbance Stability Guarantees

Luyao Niu1∗, Dinuka Sahanbandu1∗, Andrew Clark2, and Radha Poovendran1 This
work was supported by the Air Force Office of Scientific Research (AFOSR)
through grants FA9550-23-1-0208, FA9550-22-1-0054, and
FA9550-20-1-0074.*Authors contributed equally to this work.1Luyao Niu, Dinuka
Sahabandu, and Radha Poovendran are with the Network Security Lab, Department
of Electrical and Computer Engineering, University of Washington, Seattle, WA
98195-2500<EMAIL_ADDRESS>Clark is with the Department
of Electrical and Systems Engineering at Washington University in St. Louis,
St Louis, MO 63130-4899<EMAIL_ADDRESS>
###### Abstract

Disturbances may create cascading failures in power systems and lead to
widespread blackouts. Controlled islanding is an effective approach to
mitigate cascading failures by partitioning the power system into a set of
disjoint islands. To retain the stability of the power system following
disturbances, the islanding strategy should not only be minimally disruptive,
but also guarantee post-disturbance stability. In this paper, we study the
problem of synthesizing post-disturbance stability-aware controlled islanding
strategies. To ensure post-disturbance stability, our computation of islanding
strategies takes load-generation balance and transmission line capacity
constraints into consideration, leading to a hybrid optimization problem with
both discrete and continuous variables. To mitigate the computational
challenge incurred when solving the hybrid optimization program, we propose
the concepts of hybrid submodularity and hybrid matroid. We show that the
islanding problem is equivalent to a hybrid matroid optimization program,
whose objective function is hybrid supermodular. Leveraging the
supermodularity property, we develop an efficient local search algorithm and
show that the proposed algorithm achieves $\frac{1}{2}$-optimality guarantee.
We compare our approach with a baseline using mixed-integer linear program on
IEEE 118-bus, IEEE 300-bus, ActivSg 500-bus, and Polish 2383-bus systems. Our
results show that our approach outperforms the baseline in terms of the total
cost incurred during islanding across all test cases. Furthermore, our
proposed approach can find an islanding strategy for large-scale test cases
such as Polish 2383-bus system, whereas the baseline approach becomes
intractable.

###### Index Terms:

Controlled islanding, hybrid submodular, cascading failure, power system
restoration, blackstart allocation

##  I Introduction

Power systems are often operated close to their stability and capacity limits.
When power systems incur disturbances such as cyber attacks [1], natural
disasters [2], and spiking growth in demand [3], some transmission lines and
generators may experience outages, leading to overflow at neighboring
transmission lines and hence cascading failures [4]. Cascading failures can
lead to disastrous consequences. For example, the 2023 Pakistan blackout [5]
left approximately $220$ million people without power.

One approach to improve the resilience of power system and mitigate cascading
failures is through _controlled islanding_ [6]. Controlled islanding
deliberately trips a subset of transmission lines to partition the power
system into a collection of disjoint, internally stable and connected islands.
Currently, metrics such as power flow disruption, generator coherency, and
load-generation imbalance are widely adopted when computing controlled
islanding strategies [7, 8, 9, 10, 11, 12, 13, 14, 15, 16].

In addition to the aforementioned metrics, the post-disturbance stability
needs to be taken into consideration when computing the controlled islanding
strategies. For example, if a controlled islanding strategy cannot guarantee
the post-islanding power flow to be within the capacity of each transmission
line, controlled islanding may have to be executed again within one newly
formulated island to avoid cascading failure inside the island. In the worst-
case, such islanding strategies may have to be executed multiple times,
leading to widely-spread blackout. Furthermore, the controlled islanding
strategy should partition the power system in a way such that the tripped
connections can be re-established and thus the whole power system can be
restored efficiently [17, 11, 18]. To allow efficient restoration, each island
is required to have sufficient availability to blackstart generators.

At present, computing a controlled islanding strategy that jointly (i)
optimizes the widely-adopted metrics such as power flow disruption and
generator coherency, (ii) takes post-disturbance stability into consideration,
and (iii) guarantees restoration after islanding, has been less studied.
Synthesizing such an islanding strategy involves both continuous variables
(post-islanding power flow and load shedding at each load bus for load-
generation balance) and discrete variables (choices of tripping transmission
lines and allocation of blackstart generators). Such hybrid optimization
programs are generally NP-hard. Finding exact solutions to these programs is
computationally intensive, and hence does not scale to large-scale power
systems. Although there exist heuristic solution algorithms to compute
controlled islanding strategies [19, 20, 21], they omit some constraints from
(i)-(iii) and cannot provide provable optimality guarantees.

In this paper, we investigate how to compute a controlled islanding strategy
that jointly satisfies three goals. The first goal is to optimize metrics
including power flow disruption and generator coherency. Second, the islanding
strategy considers post-disturbance stability, which is captured by satisfying
load-generation balance to each island and ensuring the post-islanding power
flow on each transmission line to remain with capacity limit. The third goal
of the islanding strategy is to guarantee efficient restoration of power
system after islanding by ensuring sufficient availability of energized
blackstart generators within each island. We formulate a hybrid optimization
to compute such a controlled islanding strategy. To solve the hybrid
optimization problem, we propose a concept named _hybrid submodularity_ as a
generalization of discrete submodularity. We prove that the objective
functions and constraints are hybrid submodular, and develop an efficient
algorithm to compute the islanding strategy with provable optimality
guarantees. To summarize, this paper makes the following contributions.

  * •

We formulate the problem of synthesizing a controlled islanding strategy that
jointly optimizes power flow disruption, dynamical stability, post-disturbance
stability, and blackstart generator allocation for restoration.

  * •

We translate the islanding problem to a matroid optimization problem. We prove
that transmission line capacity constraint for post-islanding power flow and
blackstart generator allocation can be encoded by hybrid monotone and
supermodular functions.

  * •

We present a local search algorithm to compute the islanding strategy. We
develop a $\frac{1}{2}$-optimality bound for the algorithm based on hybrid
submodularity property.

  * •

We evaluate our proposed approach on IEEE 118-bus, IEEE 300-bus, ActivSg
500-bus, and Polish 2383-bus systems. Compared with a baseline approach using
mixed-integer linear program (MILP)-based formulation, our results show that
the proposed solution approach outperforms the baseline in terms of both total
cost and amount of load shedding. Moreover, our approach scales well to large-
scale test cases whereas the baseline approach becomes infeasible.

The present paper generalizes the approach proposed in [22] in the following
aspects. This paper incorporates both the post-disturbance stability and
availability to blackstart generators in each island when computing controlled
islanding strategies, which are not considered in [22]. We prove that both
constraints on transmission line capacity and blackstart generator allocation
are hybrid supermodular in the islanding strategies, and develop an efficient
solution algorithm with $\frac{1}{2}$-optimality guarantee. We implement the
proposed solution approach to large-scale power systems, namely IEEE 118-bus,
IEEE 300-bus, ActivSg 500-bus, and Polish 2383-bus systems.

The remainder of this paper is organized as follows. Section II reviews
related literature. In Section III, we present preliminary background on
submodularity and matroids. Section IV formulates the controlled islanding
problem. Section V translates the islanding problem to a matroid optimization
problem, and develops a local search algorithm based on hybrid submodularity
property with $\frac{1}{2}$-optimality guarantee. In Section VI, we evaluate
our proposed approach using IEEE 118-bus, IEEE 300-bus, ActivSg 500-bus, and
Polish 2383-bus systems. Section VII concludes the paper.

##  II Related Work

Power systems are often operated at a stringent operating point. Therefore,
disturbances such as cyber attacks [1] and natural disasters [2] could
destablize power systems and cause cascading failures. Controlled islanding
has been shown to be an effective approach to mitigate cascading failures
following disturbances [8, 20, 23].

A widely adopted solution to computing controlled islanding strategy is based
on slow coherency [7, 8, 9]. These category of approaches groups generators
using slow coherency analysis. Then an islanding strategy can be found by
partitioning the power system so that non-coherent generators are
disconnected. Slow coherency analysis-based computation of islanding
strategies requires slow eigenbasis analysis, and may not scale to power
systems of large sizes.

Spectral clustering controlled islanding [23] has been proposed to improve the
computation efficiency to obtain the islanding strategies. This approach
follows two steps, where the first step groups the generators based on their
dynamics, and the second step trips the transmission lines using metrics such
as power flow disruption or load-generation imbalance. An alternative class of
approaches to improve the scalability of controlled islanding algorithms is to
applying ordered binary decision diagram (OBDD) methods on simplified graph
representations of the power systems [19, 20]. Such approaches may overly
simplify the power systems, and do not provide provable optimality guarantees
on the obtained islanding strategies.

To obtain the exact islanding strategies, mixed-integler linear programs have
been formulated to incorporate different constraints such as power flow
disruption, load-generation imbalance, generator coherency, and power system
restoration [10, 11, 12, 13, 14, 15, 16]. However, solving NP-hard MILPs
formulated on large-scale power systems is computationally expensive. A linear
program-based microgrid formation is investigated in [24]. However, this
approach is network dependent, and may not be readily generalized to arbitrary
large-scale power systems.

Submodularity-based approaches [25, 22, 26] have recently been proposed to not
only efficiently compute controlled islanding strategies for large-scale power
systems, but also provide provable optimality guarantees for the obtained
islanding strategies. In [26], metrics including power flow disruption,
generator coherency, and load-generation imbalance are shown to be monotone
and submodular, leading to efficient local search solution algorithms with
$\frac{1}{2}$ optimality bounds. Submodular optimization has also been adopted
in other application scenarios such as the placement of energy storage units
[27], voltage control [28], and distribution network configuration [29]. The
concept of hybrid submodularity in this paper can benefit these solutions to
incorporate continuous variables such as capacities of storage units into
consideration to improve the optimalities of solutions [27, 28, 29].

##  III Preliminary Background

This section introduces preliminary background on submodularity and matroids.
Consider a finite set $\mathcal{V}$. A function
$f:2^{\mathcal{V}}\rightarrow\mathbb{R}_{\geq 0}$ is monotone nondecreasing if
$f(\mathcal{S})\leq f(\mathcal{T})$ for any $\mathcal{S}\subseteq\mathcal{T}$.
Function $f$ is said to be submodular [30] if

$f(\mathcal{S}\cup\\{v\\})-f(\mathcal{S})\geq
f(\mathcal{T}\cup\\{v\\})-f(\mathcal{T})$

holds for any $\mathcal{S}\subseteq\mathcal{T}\subseteq\mathcal{V}$ and
$v\in\mathcal{V}\setminus\mathcal{T}$. A function $f$ is supermodular if $-f$
is submodular.

We present the definition for matroids, which give rise to a class of
submodular functions, as follows.

######  Definition 1 (Matroid).

A matroid $\mathcal{M}$ is a pair $(\mathcal{V},\mathcal{I})$, where
$\mathcal{V}$ is a finite set, and $\mathcal{I}\subseteq 2^{\mathcal{V}}$ is a
collection of subsets of $\mathcal{V}$ satisfying (i)
$\emptyset\in\mathcal{I}$, (ii) if
$\mathcal{S}\subseteq\mathcal{T}\subseteq\mathcal{V}$ and
$\mathcal{T}\in\mathcal{I}$, then $\mathcal{S}\in\mathcal{I}$, and (iii) if
$\mathcal{S},\mathcal{T}\in\mathcal{I}$ and $|\mathcal{T}|>|\mathcal{S}|$,
then there exists $v\in\mathcal{T}\setminus\mathcal{S}$ such that
$\mathcal{S}\cup\\{v\\}\in\mathcal{I}$.

Each set in $\mathcal{I}$ is called an independent set of matroid
$\mathcal{M}$. A maximal independent set of $\mathcal{M}$ is a basis of
$\mathcal{M}$. The set of bases of matroid $\mathcal{M}$ is denoted as
$\mathcal{B}(\mathcal{M})$. The rank function
$\rho_{\mathcal{M}}:2^{\mathcal{V}}\rightarrow\mathbb{Z}_{\geq 0}$ is defined
as
$\rho_{\mathcal{M}}(\mathcal{T})=\max\\{|\mathcal{S}|:\mathcal{S}\subseteq\mathcal{T},\mathcal{S}\in\mathcal{I}\\}$.
The rank function of matroid $\mathcal{M}$ is nondecreasing and submodular
[31].

In this paper, we consider graphic matroid $\mathcal{M}_{G}$ induced by a
graph $G=(\mathcal{V},\mathcal{E})$, where $\mathcal{V}$ is a finite set of
vertices, and $\mathcal{E}\subseteq\mathcal{V}\times\mathcal{V}$ is the set of
edges. Graphic matroid $\mathcal{M}_{G}$ is defined as
$\mathcal{M}_{G}=(\mathcal{E},\mathcal{I})$. Each independent set in
$\mathcal{I}$ of graphic matroid $\mathcal{M}_{G}$ is acyclic. For a connected
graph $G$, the bases $\mathcal{B}(\mathcal{M}_{G})$ are the spanning trees of
graph $G$ [31].

##  IV System Model and Islanding Formulation

We consider a power system and model it as a graph, denoted as
$G=(\mathcal{V},\mathcal{E})$, where $\mathcal{V}$ is the set of vertices and
$\mathcal{E}\subseteq\mathcal{V}\times\mathcal{V}$ is the set of edges.
Without loss of generality, we assume that graph $G$ is connected. Here we use
the set of vertices $\mathcal{V}=\\{1,\ldots,n\\}$ to represent the set of
buses, and use the set of edges $\mathcal{E}$ to represent the collection of
transmission lines. Among the buses, we denote the set of generator and load
buses as $\mathcal{G}\subseteq\mathcal{V}$ and
$\mathcal{L}\subseteq\mathcal{V}$, respectively, where
$\mathcal{G}\cap\mathcal{L}=\emptyset$. For each generator bus
$j\in\mathcal{G}$, we denote its generation capacity as $\bar{g}_{j}$. The set
of generators $\mathcal{G}$ can be classified based on whether a unit has
blackstart capability or relies on cranking power from the system to restart.
If a generator bus $j$ is attached with a blackstart generator, we represent
it as $\mathbb{I}_{j}=1$. Otherwise $\mathbb{I}_{j}=0$. We further group the
set of generators $\mathcal{G}$ into $m$ coherence groups based on their
coherence [8]. We select a reference generator $r_{k}$ within each coherence
group. For each load bus $j\in\mathcal{L}$, we denote its maximum load as
$\bar{d}_{j}$. For each transmission line $(j,j^{\prime})\in\mathcal{E}$, we
denote the power flow from bus $j$ to $j^{\prime}$ through transmission line
$(j,j^{\prime})$ as $P_{jj^{\prime}}$. The power flow $P_{jj^{\prime}}$ should
satisfy $P_{jj^{\prime}}\leq\bar{P}_{jj^{\prime}}$ to avoid overflow and
transmission line outage, where $\bar{P}_{jj^{\prime}}$ is the transmission
line capacity of $(j,j^{\prime})$.

Our goal is to partition the power system $G$ into $m$ disjoint islands,
denoted as $I_{1},\ldots,I_{m}$, to mitigate cascading failures. Each island
$I_{k}$ can be represented as

$I_{k}=(\mathcal{V}_{k},\mathcal{E}_{k}),$

where $\mathcal{V}_{k}\subset\mathcal{V}$ is a subset of buses and
$\mathcal{E}_{k}=(\mathcal{V}_{k}\times\mathcal{V}_{k})\cap\mathcal{E}$ is a
subset of transmission lines. To ensure the islands to be disjoint, we have
that $\mathcal{V}_{k}\cap\mathcal{V}_{k^{\prime}}=\emptyset$ for all $k\neq
k^{\prime}$. Note that for each reference generator $r_{k}$, it must be
contained within island $I_{k}$ so that $r_{k}$ and $r_{k^{\prime}}$ are
disconnected for any $k\neq k^{\prime}$. In addition, each island $I_{k}$ is
connected for all $k=1,\ldots,m$. In the remainder of this paper, we will
denote $H=\cup_{k=1}^{m}\mathcal{E}_{k}$ as the set of transmission lines that
remain connected in the power system after formulating the islands.

In this paper, we partition the power system into $m$ number of islands by
optimizing generator coherency and power flow disruption. We formulate each
metric in detail below.

We denote the coherence matrix as
$\bar{A}\in\mathbb{R}^{m\times|\mathcal{G}|}$ (see Appendix for a detailed
derivation of $\bar{A}$.) The generator coherency is given as

$F_{1}(H)=\|\bar{A}-A(H)\|_{F},$ (1)

where $\|\cdot\|_{F}$ represents the Frobenius norm, and
$A(H)\in\mathbb{R}^{|\mathcal{G}|\times m}$ whose $(j,k)$-th entry is defined
as $\left[A(H)\right]_{jk}=1$ if $j\in\mathcal{V}_{k}$ and zero otherwise.

Power flow disruption measures the total amount of power flow being disrupted
due to tripping transmission lines to form islands. The power flow disruption
can be computed as

$F_{2}(H)=\sum_{(j,j^{\prime})\in\mathcal{E}\setminus
H}\frac{|P_{jj^{\prime}}|}{2},$ (2)

where $P_{jj^{\prime}}$ is the pre-islanding power flow.

In addition to generator coherency and power flow disruption, the formed
islands $I_{1},\ldots,I_{m}$ need to guarantee post-disturbance stability. In
particular, we will ensure load-generation balance and transmission line
capacity constraints on the post-islanding power flow.

Load-generation balance requires that the amount of power generation within
each island should be sufficient to meet the total power demand. Due to the
islanding operation, there may exist deficit in power generation, requiring
load shedding from some buses. We define
$c_{j}:\mathbb{R}\rightarrow\mathbb{R}$ as the cost function of load shedding
for load bus $j$. We assume that $c_{j}$ is monotone increasing and convex.
Then given set $H$, the loads $d=[d_{1},\ldots,d_{|\mathcal{L}|}]^{\top}$ of
buses in $\mathcal{L}$ can be found by

$\displaystyle F_{3}(H,d)=$
$\displaystyle\min_{d}\quad\sum_{j\in\mathcal{L}}c_{j}(\bar{d}_{j}-d_{j})$
(3a)
$\displaystyle\text{s.t.}\quad\sum_{j\in\mathcal{L}\cap\mathcal{V}_{k}}d_{j}\leq\sum_{j\in\mathcal{G}\cap\mathcal{V}_{k}}\bar{g}_{j},~{}\forall
k=1,\ldots,m$ (3b) $\displaystyle\quad\quad\quad
d_{j}\in[0,\bar{d}_{j}],~{}\forall j\in\mathcal{L}.$ (3c)

Any feasible solution to optimization program (3) guarantees load-generation
balance. In the remainder of this paper, we define $(H,d)$ as an islanding
strategy, which specifies the set of transmission lines that will remain in
the system as well as the amount of load shedding incurred by the load buses.

We next focus on the transmission line capacity constraint for post-islanding
power flow, denoted as $\tilde{P}_{jj^{\prime}}$. The post-islanding power
flow should satisfy the conservation law

$\sum_{j^{\prime}:(j,j^{\prime})\in\mathcal{E}}\tilde{P}_{jj^{\prime}}+d_{j}-g_{j}=0,~{}\forall
j\in\mathcal{V}.$ (4)

Furthermore, $\tilde{P}_{jj^{\prime}}$ needs to satisfy the transmission line
capacity constraint given as below

$\tilde{P}_{jj^{\prime}}\leq\bar{P}_{jj^{\prime}},~{}(j,j^{\prime})\in H,$ (5)

in order to ensure the post-disturbance stability.

To allow the power system to be restored, each island should contain at least
one blackstart generator. This constraint is formulated as follows

$\sum_{j\in\mathcal{V}_{k}}\mathbb{I}_{j}\geq 1,~{}\forall k=1,\ldots,m.$ (6)

Let $\alpha_{i}\in(0,1)$ for $i=1,2,3$ be constants modeling the trade-off
parameters among metrics $F_{1}(H)$, $F_{2}(H)$, and $F_{3}(H,d)$,
respectively. We then formulate the following optimization program to compute
the islanding strategy $(H,d)$

$\displaystyle\min_{E,d}\quad$
$\displaystyle\alpha_{1}F_{1}(H)+\alpha_{2}F_{2}(H)+\alpha_{3}F_{3}(H,d)$ (7a)
s.t. $\displaystyle r_{k}\in I_{k},~{}\forall k=1,\ldots,m$ (7b)
$I_{1},\ldots,I_{m}$ form a set of islands (7c) Eqn. (3b), (3c), (4), (5), and
(6) (7d)

##  V Hybrid Submodularity and Proposed Solution

Solving optimization problem (7) is computationally intractable for large-
scale power systems since it involves coupled continuous and discrete
variables. Furthermore, the solution space of the discrete variable $H$ grows
exponentially with respect to the size of graph $G$.

In this section, we present a computationally efficient solution approach to
compute an islanding strategy with provable optimality guarantee. We first
propose a concept named hybrid submodularity. Next, we prove that the metrics
($F_{1}(H)$, $F_{2}(H)$, and $F_{3}(H,d)$) and the constraints (Eqn. (3b),
(3c), (5), and (6)) satisfy hybrid submodular properties. We finally present a
local search algorithm to compute the islanding strategy and prove the
optimality guarantee based on the hybrid submodular property.

###  V-A Hybrid Submodularity

This subsection introduces the concepts of hybrid monotonicity and hybrid
submodularity.

######  Definition 2 (Hybrid Monotonocity).

Let $\mathcal{V}$ be a finite set and $\mathcal{D}\subset\mathbb{R}^{n}$. A
function
$f:2^{\mathcal{V}}\times\mathcal{F}(\mathcal{D})\rightarrow\mathbb{R}$ is
hybrid monotone nondecreasing if, for any finite sets
$\mathcal{S}\subseteq\mathcal{T}\subseteq\mathcal{V}$ and
$\Lambda\subseteq\Lambda^{\prime}\subseteq\mathcal{D}$, we have

$f(\mathcal{S},\Lambda)\leq f(\mathcal{T},\Lambda^{\prime}),$ (8)

where $\mathcal{F}(\mathcal{D})$ represents the collection of finite subsets
of $\mathcal{D}$. A function $f$ is hybrid monotone nonincreasing if $-f$ is
hybrid monotone nondecreasing.

We next define hybrid submodularity which generalizes the classic
submodularity property defined over discrete sets by incorporating infinite
ground set.

######  Definition 3 (Hybrid Submodularity).

Let $\mathcal{V}$ be a finite set and $\mathcal{D}\subset\mathbb{R}^{n}$. A
function
$f:2^{\mathcal{V}}\times\mathcal{F}(\mathcal{D})\rightarrow\mathbb{R}$ is
hybrid submodular if, for any
$\mathcal{S}\subseteq\mathcal{T}\subseteq\mathcal{V}$ and
$\Lambda\subseteq\Lambda^{\prime}\subseteq\mathcal{D}$, the following
properties hold:

  1. 1.

For any $j\in\mathcal{V}\setminus\mathcal{T}$, we have

$f(\mathcal{S}\cup\\{j\\},\Lambda)-f(\mathcal{S},\Lambda)\geq
f(\mathcal{T}\cup\\{j\\},\Lambda^{\prime})-f(\mathcal{T},\Lambda^{\prime}).$

  2. 2.

For any $\lambda\in\mathcal{D}\setminus\Lambda^{\prime}$, we have

$f(\mathcal{S},\Lambda\cup\\{\lambda\\})-f(\mathcal{S},\Lambda)\geq
f(\mathcal{T},\Lambda^{\prime}\cup\\{\lambda\\})-f(\mathcal{T},\Lambda^{\prime}).$

We say a function $f$ is hybrid supermodular if $-f$ is hybrid submodular. We
finally present hybrid matroid as a generalization of matroids to discrete
sets.

######  Definition 4 (Hybrid Matroid).

A hybrid matroid is a tuple
$\mathcal{M}=(\mathcal{V},\mathcal{D},\mathcal{I})$, where $\mathcal{V}$ is a
finite set, $\mathcal{D}\subseteq\mathbb{R}^{n}$, and $\mathcal{I}$ is a
collection of subsets of $\mathcal{V}\times\mathcal{D}$ such that

  1. 1.

$(\emptyset,\emptyset)\in\mathcal{I}$.

  2. 2.

if $(\mathcal{T},\Lambda^{\prime})\in\mathcal{I}$,
$\mathcal{S}\subseteq\mathcal{T}$, and $\Lambda\subseteq\Lambda^{\prime}$,
then $(\mathcal{S},\Lambda)\in\mathcal{I}$.

  3. 3.

if $(\mathcal{S},\Lambda),(\mathcal{T},\Lambda^{\prime})\in\mathcal{I}$, where
$|\mathcal{S}|<|\mathcal{T}|$ and $|\Lambda|<|\Lambda^{\prime}|$, then there
exists $v\in\mathcal{T}\setminus\mathcal{S}$ and
$\lambda\in\Lambda^{\prime}\setminus\Lambda$ such that
$(\mathcal{S}\cup\\{v\\},\lambda\cup\\{\lambda\\})\in\mathcal{I}$.

If set $(\mathcal{S},\Lambda)\in\mathcal{I}$, we then say
$(\mathcal{S},\Lambda)$ is an indpendent set of hybrid matroid $\mathcal{M}$.
A maximal independent set of $\mathcal{M}$ is a basis of $\mathcal{M}$. We
denote the set of bases of $\mathcal{M}$ as $\mathcal{B}(\mathcal{M})$.

Let a function $f$ be hybrid monotone nondecreasing and submodular and
$\mathcal{M}$ be a hybrid matroid. Then we can derive a local optimality
guarantee for function $f$, as stated below.

######  Proposition 1 ([22]).

Suppose that function
$f:2^{\mathcal{V}}\times\mathcal{F}(\mathcal{D})\rightarrow\mathbb{R}$ is
hybrid monotone nondecreasing and submodular. Let $\mathcal{M}$ be a hybrid
matroid and $(\mathcal{S},\Lambda)\in\mathcal{B}(\mathcal{M})$. For any
$v\in\mathcal{S}$, $t\notin\mathcal{S}$, $\lambda\in\Lambda$, and
$\lambda^{\prime}\notin\Lambda$ satisfying
$(\mathcal{S}\setminus\\{v\\}\cup\\{t\\},\Lambda\setminus\\{\lambda\\}\cup\\{\lambda^{\prime}\\})$,
if we have

$f(\mathcal{S},\Lambda)\geq
f(\mathcal{S}\setminus\\{v\\}\cup\\{t\\},\Lambda\setminus\\{\lambda\\}\cup\\{\lambda^{\prime}\\}),$

then for any $(\mathcal{T},\Lambda^{\prime})\in\mathcal{B}(\mathcal{M})$,
$f(\mathcal{S},\Lambda)\geq\frac{1}{2}f(\mathcal{T},\Lambda^{\prime})$.

###  V-B Reformulating Eqn. (7) as Hybrid Matroid Optimization

This subsection relaxes the formulation in Eqn. (7), and converts the problem
to a hybrid matroid optimization program. We show that the objective function
of the hybrid matroid optimization program is hybrid supermodular.

We augment the power system $G$ to
$\bar{G}=(\bar{\mathcal{V}},\bar{\mathcal{E}})$, where
$\bar{\mathcal{V}}=\mathcal{V}\cup\\{a\\}$ is obtained by introducing a
supernode $a$, and

$\bar{\mathcal{E}}=\mathcal{E}\cup\\{(a,r_{k}):k=1,\ldots,m\\}.$ (9)

The set of edges $\\{(a,r_{k}):k=1,\ldots,m\\}$ connects the supernode $a$
with each reference generator $r_{i}$.

Since the reference generators belong to different coherence groups, they
should be partitioned into different islands, i.e., there exist no $r_{k}$ and
$r_{k^{\prime}}$ such that $r_{k},r_{k^{\prime}}\in\mathcal{V}_{k}$ for any
$k\neq k^{\prime}$. Based on this insight, we have the following preliminary
result.

######  Lemma 1.

If set $H\cup\\{(a,r_{k}):k=1,\ldots,m\\}$ forms a spanning tree of $\bar{G}$,
then islands $I_{1},\ldots,I_{m}$ induced by $H$ on power system $G$ satisfy
constraints (7b) and (7c).

###### Proof.

We prove the lemma by contradiction. Suppose that set
$H\cup\\{(a,r_{k}):k=1,\ldots,m\\}$ forms a spanning tree of $\bar{G}$ while
there exist reference generators
$r_{k},r_{k^{\prime}}\in\mathcal{V}_{k^{\prime\prime}}$ for some $k\neq
k^{\prime}$. Since $r_{k},r_{k^{\prime}}\in\mathcal{V}_{k^{\prime\prime}}$ and
island $I_{k^{\prime\prime}}$ is connected to be an island, we have that there
exists some path $r_{k},j_{1},\ldots,j_{t},r_{k^{\prime}}$ that connects
reference generators $r_{k}$ and $r_{k^{\prime}}$ without going through
supernode $a$, i.e., $j_{l}\neq a$ for all $l=1,\ldots,t$. By the construction
of graph $\bar{G}$, reference generators $r_{k}$ and $r_{k^{\prime}}$ are also
connected by another path $r_{k},a,r_{k^{\prime}}$ that goes through supernode
$a$. Hence, paths $r_{k},a,r_{k^{\prime}}$ and
$r_{k},j_{1},\ldots,j_{t},r_{k^{\prime}}$ form a cycle, which contradicts the
definition of spanning tree. Therefore, we have that if set
$H\cup\\{(a,r_{k}):k=1,\ldots,m\\}$ forms a spanning tree of $\bar{G}$, then
constraint (7b) is satisfied.

Since set $H\cup\\{(a,r_{k}):k=1,\ldots,m\\}$ forms a spanning tree of graph
$\bar{\mathcal{G}}$, we have that $H$ is a spanning forest for power system
$G$ containing $m$ trees. If $H$ contains more that $m$ trees, then
$H\cup\\{(a,r_{k}):k=1,\ldots,m\\}$ cannot form a spanning tree since there
must exists some nodes that are disconnected from the spanning tree. If $H$
contains less that $m$ trees, then there must exist some reference generators
$r_{k},r_{k^{\prime}}\in I_{k^{\prime\prime}}$, which has been falsified
earlier in the proof. We thus have that
$\mathcal{V}_{k}\cap\mathcal{V}_{k^{\prime}}=\emptyset$ for any $k\neq
k^{\prime}$. Furthermore, we have that each tree $I_{k}$ contained in the
spanning forest induced by $H$ must be connected by using the definition of
spanning tree. Combining these arguments yields the result that the islands
$I_{1},\ldots,I_{m}$ satisfy constraints (7c). ∎

Lemma 1 allows us to encode constraints (7b) and (7c) using a graphic matroid
constraint $E\in\mathcal{B}(\mathcal{M}_{\bar{G}})$, where
$\mathcal{M}_{\bar{G}}$ represents the graphic matroid of $\bar{G}$ and $E$ is
the set of edges in the spanning tree. We define

$\bar{\mathcal{D}}=\prod_{k=1}^{m}\prod_{j\in\mathcal{L}}[0,\bar{d}_{j}],$

and a hybrid matroid
$\bar{\mathcal{M}}=(\bar{\mathcal{E}},\bar{\mathcal{D}},\bar{\mathcal{I}})$.
We can thus restrict ourselves to $(E,\Lambda)\in\bar{\mathcal{I}}$ to search
for islanding strategies, where
$E\in\mathcal{B}(\mathcal{M}_{\bar{\mathcal{G}}})$, $d\in\Lambda$, and
$|\Lambda|\leq 1$. Therefore, we translate optimization program (7) to the
following constrained matroid optimization program

$\displaystyle\min_{E,d}\quad$
$\displaystyle\alpha_{1}F_{1}(E)+\alpha_{2}F_{2}(E)+\alpha_{3}F_{3}(E,d)$
(10a) s.t. $\displaystyle(E,d)\in\mathcal{B}(\bar{\mathcal{M}})$ (10b) Eqn.
(3b), (3c), (4), (5), and (6) (10c)

###  V-C Hybrid Sumodularity-based Algorithm Development

Solving the constrained matroid optimization in Eqn. (10) is still challenging
due to the presence of coupled discrete and continuous variables. In what
follows, we prove that the objective function and constraint (10c) are hybrid
supermodular in the islanding strategy. Using the hybrid supermodularity
result, we then develop a local search algorithm to efficiently compute the
islanding strategy. We conclude this section by presenting the optimality
guarantee of our developed algorithm.

We define an indicator function $\chi_{kj}(E)$ for each $k=1,\ldots,m$ and
$j\in\mathcal{V}$ as follows

$\chi_{kj}(E)=\begin{cases}1,&\mbox{ if bus $j$ is in island $I_{k}$,}\\\
0,&\mbox{ otherwise}.\end{cases}$

We have the following preliminary result.

######  Lemma 2.

Let $\mathcal{M}_{kj}$ be the graphic matroid of graph
$(\mathcal{V},E\cup\\{(r_{k},j)\\})$. There is a nonincreasing supermodular
function
$\bar{\chi}_{kj}(E)=n-m-\rho_{\mathcal{M}_{kj}}(E\cup\\{(r_{k},j)\\})$ such
that $\chi_{kj}(E)=\bar{\chi}_{kj}(E)$ for any
$E\in\mathcal{B}(\mathcal{M}_{\bar{G}})$. Furthermore, there are nonincreasing
supermodular functions $\bar{F}_{1}(E)$ and $\bar{F}_{2}(E)$ defined as

$\displaystyle\bar{F}_{1}(E)$
$\displaystyle=\sum_{k=1}^{m}\sum_{j\in\mathcal{G}}(\bar{A}_{kj}+(1-2\bar{A}_{kj})\Omega_{kj}(E)),$
$\displaystyle\bar{F}_{2}(E)$ $\displaystyle=\sum_{k^{\prime}\neq
k}\sum_{(j,j^{\prime})\in\mathcal{E}}\phi_{k,k^{\prime},j,j^{\prime}}(E)|P_{jj^{\prime}}|,$

such that $\bar{F}_{1}(E)=F_{1}(E)$ and $\bar{F}_{2}(E)=F_{2}(E)$ for all
$E\in\mathcal{B}(\mathcal{M}_{\bar{G}})$, where

$\displaystyle\Omega_{kj}(E)$
$\displaystyle=\begin{cases}\bar{\chi}_{kj}(E),&\mbox{ if
$\bar{A}_{kj}\leq\frac{1}{2}$},\\\ 1-\sum_{k^{\prime}\neq
k}\bar{\chi}_{k^{\prime}j}(E),&\mbox{ if
$\bar{A}_{kj}>\frac{1}{2}$}.\end{cases}$
$\displaystyle\phi_{k,k^{\prime},j,j^{\prime}}(E)$
$\displaystyle=\begin{cases}1,&\mbox{ if $j\in I_{k}$ and $j^{\prime}\in
I_{k^{\prime}}$},\\\ 0,&\mbox{ otherwise}.\end{cases}$

###### Proof.

The monotonicty and supermodularity of function $\bar{\chi}$ follow from the
monotonicity and submodularity of the rank function of a matroid [31]. When
$E\in\mathcal{B}(\mathcal{M}_{\bar{G}})$, then we have that
$\chi_{kj}(E)=1-(\rho_{\mathcal{M}_{kj}}(E\cup\\{(r_{k},j)\\})-\rho_{\mathcal{M}}(E))$,
where $\mathcal{M}$ is the graphic matroid of $G$ [26]. Note that
$\rho_{\mathcal{M}}(E)=n-m-1$ since the number of edges in a spanning tree is
equal to the number of nodes minus one. We thus have that
$\bar{\chi}_{kj}(E)=\chi_{kj}(E)=1-(\rho_{\mathcal{M}_{kj}}(E\cup\\{(r_{k},j)\\})-\rho_{\mathcal{M}}(E))$.

Let $\bar{A}$ be the coherence matrix [32]. We then show that there exists a
nonincreasing and supermodular function $\bar{F}_{1}(E)$ such that
$\bar{F}_{1}(E)=F_{1}(E)$ for all $E\in\mathcal{B}(\mathcal{M}_{\bar{G}})$.
Since $\bar{\chi}_{kj}(E)$ is supermodular, we have that $\Omega_{kj}(E)$ is
also supermodular when $\bar{A}_{kj}\leq\frac{1}{2}$, making
$(1-2\bar{A}_{kj})\Omega_{kj}(E)$ further supermodular. When
$\bar{A}_{kj}>\frac{1}{2}$, then $-(1-2\bar{A}_{kj})>0$, making
$(1-2\bar{A}_{kj})\Omega_{kj}(E)=(1-2\bar{A}_{kj})-(1-2\bar{A}_{kj})\sum_{k^{\prime}\neq
k}\bar{\chi}_{k^{\prime}j}(E)$ supermodular. The statement that
$\bar{F}_{1}(E)=F_{1}(E)$ can be verified by the definition of Frobenius norm
and the fact that $\sum_{k=1}^{m}\chi_{kj}(E)=1$ for any set $E$ that forms
islands.

One can verify that function $\bar{F}_{2}(E)$ is nonincreasing and supmodular
when $E\in\mathcal{B}(\mathcal{M}_{\bar{G}})$ by showing the supermodularity
of $\phi_{k,k^{\prime},j,j^{\prime}}(E)$, which can be found in [26]. ∎

Lemma 2 indicates that metrics $F_{1}(E)$ and $F_{2}(E)$ can be rewritten
using monotone and supermodular functions $\bar{F}_{1}(E)$ and
$\bar{F}_{2}(E)$, respectively. In the following, we focus on metric
$F_{3}(E,d)$ which involves both continuous and discrete variables. We define
a set of auxiliary variables $d_{kj}$ for each $k=1,\ldots,m$ and
$j\in\mathcal{L}$ to model the amount of load attached to bus $j$ when
contained in island $I_{k}$. By the definition of variable $\chi_{kj}(E)$ and
$d_{kj}$, we have that $d_{kj}=0$ must hold when $\chi_{kj}(E)=0$. As a
consequence, we relax metric $F_{3}(E,d)$ by using a function
$\bar{F}_{3}(E,d)$ as

$\bar{F}_{3}(E,d)=\min_{d\in\Lambda}\sum_{j\in\mathcal{L}}c_{j}(\bar{d}_{j}-\sum_{k=1}^{m}d_{kj}).$
(11)

We next show that function $\bar{F}_{3}(E,d)$ is hybrid monotone and
supermodular in the islanding strategy.

######  Lemma 3.

Function $\bar{F}_{3}(E,d)$ is hybrid monotone nonincreasing and hybrid
supermodular in islanding strategy $(E,d)$.

###### Proof.

The lemma holds by the fact that any function of the form
$f(\mathcal{S})=\min_{i\in\mathcal{S}}c_{i}$ is monotone nonincreasing and
supermodular. ∎

In what follows, we focus on the constraints given in Eqn. (3b), (3c), (4),
(5), and (6). We relax the constraints by encoding them into penalty
functions, and show that the penalty functions are hybrid monotone and
supermodular.

We first consider load-generation balance constraint (3b). We define a penalty
function $F_{4}(E,d)$ given as below

$F_{4}(E,d)=\sum_{k=1}^{m}\left\\{\sum_{j\in\mathcal{L}}d_{kj}+\sum_{j\in\mathcal{G}}\sum_{k^{\prime}\neq
k}\chi_{k^{\prime}j}(E)\bar{g}_{j}-\sum_{j\in\mathcal{G}}\bar{g}_{j}\right\\}_{+}$
(12)

to penalize the amount of load that exceeds the total power generation within
each island, where $\\{\cdot\\}_{+}$ denotes $\max\\{\cdot,0\\}$.

We then consider constraint (3c). We relax constraint (3c), by introducing a
penalty function

$F_{5}(E,d)=\sum_{k=1}^{m}\sum_{j\in\mathcal{L}}\left\\{d_{kj}+\sum_{k^{\prime}\neq
k}\chi_{k^{\prime}j}(E)\bar{d}_{j}-\bar{d}_{j}\right\\}_{+},$ (13)

which imposes a positive penalty when the load attached to each load bus $j$
exceeds the limit $\bar{d}_{j}$.

We next consider constraint (5). We modify $\bar{G}$ by further introducing a
sink node $s$ and a set of edges $\mathcal{A}=\\{(j,s):j\in\mathcal{L}\\}$.
Edge set $\mathcal{A}$ connects each load bus $j\in\mathcal{L}$ with the sink
node $s$. We further define the capacities for the edges in $\mathcal{A}$ to
be infinity, and thus no capacity constraint can be violated by any
$(j,s)\in\mathcal{A}$. We then have the following result.

######  Proposition 2.

Consider the augmented graph $\bar{G}$ with sink node $s$ and edge set
$\mathcal{A}$. Given the load $d_{j}\leq\bar{d}_{j}$ for all
$j\in\mathcal{L}$, we have that there exists some post-islanding power flow
satisfying $\tilde{P}_{jj^{\prime}}\leq\bar{P}_{jj^{\prime}}$ for all
$(j,j^{\prime})\in\bar{\mathcal{E}}$ and conservation law given in Eqn. (4).

###### Proof.

Suppose that there exists some bus $j\in\mathcal{V}$ such that constraints (4)
and (5) cannot be satisfied simultaneously. We denote the post-islanding power
flow that creates capacity constraint violation as $\hat{P}_{jj^{\prime}}$.
Without loss of generality, we consider that bus $j$ is inlcuded in island
$I_{k}$. We let $j^{\prime\prime}\in\mathcal{V}_{k}\cap\mathcal{L}$ be some
load bus that is contained in island $I_{k}$. Since each island is connected,
there must exist some path $\eta$ from bus $j$ to $j^{\prime\prime}$. With a
slight abuse of notation, we use $(t,j^{\prime})\in\eta$ to represent that a
transmission line $(t,j^{\prime})$ is on path $\eta$.

By introducing the sink node $s$, we can ensure that constraints (4) and (5)
are satisfied by first choosing power flow $\tilde{P}_{jj^{\prime}}$
satisfying (5) for all $j^{\prime}$ along the path from $j$ to
$j^{\prime\prime}$, and then letting
$\tilde{P}_{j^{\prime\prime}s}=\sum_{(t,j^{\prime})\in\eta}(\hat{P}_{tj^{\prime}}-\tilde{P}_{tj^{\prime}})$
such that Eqn. (4) is met. Such power flow $\tilde{P}_{j^{\prime\prime}s}$
always exists since $\tilde{P}_{j^{\prime\prime}s}\in\mathbb{R}$. ∎

Proposition 2 allows us to verify the satisfaction of constraint (5) by
verifying whether $\sum_{j\in\mathcal{L}}|P_{js}|=0$ holds or not. We observe
that the sink node and edge set $\mathcal{A}$ do not exist in the power system
$G$, and hence we can view all auxiliary transmission line being tripped by
the islanding strategy. Therefore, we have that the transmission line capacity
constraint can be penalized by the total amount of power flow disruption over
the edges in $\mathcal{A}$. We modify function $\bar{F}_{2}(E)$ as

$\bar{F}_{2}(E)=\sum_{k^{\prime}\neq
k}\sum_{(j,j^{\prime})\in\mathcal{E}}\phi_{k,k^{\prime},j,j^{\prime}}(E)|P_{jj^{\prime}}|+\beta\sum_{j\in\mathcal{L}}|\tilde{P}_{js}|$
(14)

to jointly penalize (i) the power flow disruption due to tripping transmission
lines, and (ii) any violation of the capacity constraints in Eqn. (5), where
$\beta>0$ is a constant modeling the trade-off between (i) and (ii).

Consider constraint (6) for the blackstart generator allocation. We define a
penalty function

$F_{6}(E)=\Big{\\{}\sum_{k^{\prime}\neq
k}\sum_{j\in\mathcal{V}_{k^{\prime}}}\chi_{k^{\prime}j}(E)\mathbb{I}_{j}+\sum_{k^{\prime}\neq
k}\sum_{j\in\mathcal{V}_{k^{\prime}}}\chi_{k^{\prime}j}(E)\\\
-\sum_{j\in\mathcal{V}}\mathbb{I}_{j}+1\Big{\\}}_{+},$ (15)

to penalize the scenarios where some island contains no blackstart generator.

Using the penalty functions given in Eqn. (12) to (15), we reformulate the
constrained matroid optimization program in Eqn. (10) as the following
unconstrained program

$\displaystyle\min_{E,d}\quad$
$\displaystyle\alpha_{1}\bar{F}_{1}(E)+\alpha_{2}\bar{F}_{2}(E)+\alpha_{3}\bar{F}_{3}(E,d)+\alpha_{4}F_{4}(E,d)$
$\displaystyle\quad\quad+\alpha_{5}F_{5}(E,d)+\alpha_{6}F_{6}(E)$ (16a) s.t.
$\displaystyle(E,d)\in\mathcal{B}(\bar{\mathcal{M}})$ (16b)

where parameters $\alpha_{1}$ to $\alpha_{6}$ are positive constants modeling
the trade-off among metrics $\bar{F}_{1}(E)$, $\bar{F}_{2}(E)$,
$\bar{F}_{3}(E,d)$ as well the constraints in Eqn. (3b), (3c), (5), and (6).

In what follows, we convert the matroid optimization program in Eqn. (16) to
an equivalent hybrid submodular optimization problem. We define

$\displaystyle\bar{F}_{4}(E,d)$
$\displaystyle=\sum_{k=1}^{m}\Big{\\{}\min_{d\in\Lambda}\sum_{j\in\mathcal{L}}d_{kj}+\sum_{j\in\mathcal{G}}\sum_{k^{\prime}\neq
k}\chi_{k^{\prime}j}(E)\bar{g}_{j}-\sum_{j\in\mathcal{G}}\bar{g}_{j}\Big{\\}}_{+}$
$\displaystyle\bar{F}_{5}(E,d)$
$\displaystyle=\sum_{k=1}^{m}\sum_{j\in\mathcal{L}}\Big{\\{}\min_{d\in\Lambda}d_{kj}+\sum_{k^{\prime}\neq
k}\chi_{k^{\prime}j}(E)\bar{d}_{j}-\bar{d}_{j}\Big{\\}}_{+}$
$\displaystyle\bar{F}_{6}(E)$
$\displaystyle=\sum_{k=1}^{m}\Big{\\{}\sum_{k^{\prime}\neq
k}\sum_{j\in\mathcal{V}_{k^{\prime}}}\chi_{k^{\prime}j}(E)\mathbb{I}_{j}$
$\displaystyle\quad\quad\quad\quad+\sum_{k^{\prime}\neq
k}\sum_{j\in\mathcal{V}_{k^{\prime}}}\chi_{k^{\prime}j}(E)-\sum_{j\in\mathcal{V}}\mathbb{I}_{j}+1\Big{\\}}_{+}$
$\displaystyle\bar{F}(E,\Lambda)$
$\displaystyle=\alpha_{1}\bar{F}_{1}(E)+\alpha_{2}\bar{F}_{2}(E)+\alpha_{3}\bar{F}_{3}(E,d)$
$\displaystyle\quad\quad\quad\quad+\alpha_{4}\bar{F}_{4}(E,d)+\alpha_{5}\bar{F}_{6}(E,d)+\alpha_{6}\bar{F}_{6}(E).$

Given the definition of function $\bar{F}(E,\Lambda)$, we reformulate the
optimization problem in Eqn. (16) as

$\displaystyle\min\quad$ $\displaystyle\bar{F}(E,\Lambda)$ (17a) s.t.
$\displaystyle(E,\Lambda)\in\mathcal{B}(\bar{\mathcal{M}})$ (17b)

The equivalence between optimization programs (17) and (16) is established as
follows.

######  Proposition 3.

An islanding strategy $(E,d)$ is an optimal solution to Eqn. (16) if and only
if $(E,\Lambda)$ is an optimal solution to Eqn. (17), where $\Lambda=\\{d\\}$.
Furthermore, $\bar{F}(E,\Lambda)=F(E,d)$.

###### Proof.

The equivalence between $(E,d)$ and $(E,\Lambda)$ is established by
$\Lambda=\\{d\\}$. One can verify that $\bar{F}(E,\Lambda)=F(E,d)$ when
$\Lambda=\\{d\\}$ by using the one to one correspondence between each term in
$\bar{F}(E,\Lambda)$ and Eqn. (16a). ∎

We further establish the following monotonicity and supermodularity properties
for $\bar{F}(E,\Lambda)$ as given below.

######  Theorem 1.

The function $\bar{F}(E,\Lambda)$ is hybrid monotone nonincreasing and
supermodular in $(E,\Lambda)$.

###### Proof.

We prove the theorem by showing that each term in $\bar{F}(E,\Lambda)$ is
hybrid monotone nonincreasing and supermodular.

The monotonicity and supermodularity of functions $\bar{F}_{1}(E)$ and
$\bar{F}_{2}(E)$ hold by Lemma 2. The monotonicity and supermodularity of
$\alpha_{3}\bar{F}_{3}(E,d)$ follows from Lemma 3.

Using Lemma 2 and Definition 3, we have that
$\min_{d\in\Lambda}\sum_{j\in\mathcal{L}}d_{kj}+\sum_{j\in\mathcal{G}}\sum_{k^{\prime}\neq
k}\chi_{k^{\prime}j}(E)\bar{g}_{j}-\sum_{j\in\mathcal{G}}\bar{g}_{j}$ is
hybrid monotone nonincreasing and supermodular. In addition, we have that
function $\max\\{f(S),c\\}$ is monotone nonincreasing and supermodular for any
constant $c$ if $f(S)$ is monotone nonincreasing and supermodular. We thus
have that

$\alpha_{4}\sum_{k=1}^{m}\Big{\\{}\min_{d\in\Lambda}\sum_{j\in\mathcal{L}}d_{kj}+\sum_{j\in\mathcal{G}}\sum_{k^{\prime}\neq
k}\chi_{k^{\prime}j}(E)\bar{g}_{j}-\sum_{j\in\mathcal{G}}\bar{g}_{j}\Big{\\}}_{+}$

is monotone nonincreasing and supermodular. Similar arguments can be used to
show that the term

$\alpha_{5}\sum_{k=1}^{m}\sum_{j\in\mathcal{L}}\Big{\\{}\min_{d\in\Lambda}d_{kj}+\sum_{k^{\prime}\neq
k}\chi_{k^{\prime}j}(E)\bar{d}_{j}-\bar{d}_{j}\Big{\\}}_{+}$

is hybrid monotone nonincreasing and supermodular.

The hybrid monotonicity and supermodularity of the term

$\alpha_{6}\sum_{k=1}^{m}\Big{\\{}\sum_{k^{\prime}\neq
k}\sum_{j\in\mathcal{V}_{k^{\prime}}}\chi_{k^{\prime}j}(E)\mathbb{I}_{j}+\sum_{k^{\prime}\neq
k}\sum_{j\in\mathcal{V}_{k^{\prime}}}\chi_{k^{\prime}j}(E)\\\
-\sum_{j\in\mathcal{V}}\mathbb{I}_{j}+1\Big{\\}}_{+}$

follows by Lemma 2.

Applying the argument that linear combinations of monotone and supermodular
functions with positive weights are still monotone and supermodular completes
the proof. ∎

Algorithm 1 Local search algorithm for controlled islanding

1:Set $\epsilon\in[0,0.5)$

2:Initialize $(E,d)$ to define a valid islanding strategy

3:$found\leftarrow 1$

4:while $found==1$ do

5: $found\leftarrow 0$

6: for $(j_{1},j_{1}^{\prime})\in\mathcal{E}\setminus E,\
(j_{2},j_{2}^{\prime})\in E$ do

7: if
$E\cup\\{(j_{1},j_{1}^{\prime})\\}\setminus\\{(j_{2},j_{2}^{\prime})\\}\in\mathcal{B}(\mathcal{M}_{0})$
then

8: Compute $d^{\prime}$ and $\tilde{P}_{jj^{\prime}}$ using Eqn. (18)

9: if
$\overline{F}(E\cup\\{(j_{1},j_{1}^{\prime})\\}\setminus\\{(j_{2},j_{2}^{\prime})\\},\\{d^{\prime}\\})<(1-\epsilon)\bar{F}(E,\\{d\\})$
then

10: $E\leftarrow
E\cup\\{(j_{1},j_{1}^{\prime})\\}\setminus\\{(j_{2},j_{2}^{\prime})\\}$,
$d\leftarrow d^{\prime}$

11: $found\leftarrow 1$

12: Break

13: end if

14: end if

15: end for

16:end while

17:return $E,d$

The monotonicity and supermodularity properties established by Theorem 1 allow
us to develop an efficient local search solution algorithm, as shown in
Algorithm 1, with provable optimality guarantee. The algorithm first
initializes an islanding strategy $(E,d)$ in line 2. Note that the initial set
$E$ should not include any transmission line from the auxiliary transmission
line set $\mathcal{A}$. The algorithm then proceeds in an iterative manner. At
each iteration, the algorithm generates a new islanding strategy by including
an unselected transmission line $(j_{1},j_{1}^{\prime})\in\mathcal{E}\setminus
E$ and excluding a selected transmission line $(j_{2},j_{2}^{\prime})$, as
shown in line 7. Note that transmission line $(j_{1},j_{1}^{\prime})$ will
connect two disjoint islands, and $(j_{2},j_{2}^{\prime})$ will be used to
ensure the $m$ islands are generated with all reference generators being not
connected with each other. Then the algorithm computes the continuous variable
$d^{\prime}$ and $P$ by solving the following linear program

$\displaystyle\min_{d^{\prime},\tilde{P}}\quad$ $\displaystyle
F_{3}(E,d^{\prime})+\gamma\sum_{j\in\mathcal{V}}|P_{js}|$ (18a) s.t. Eqn.
(3b), (3c), (4) and (5) (18b)

where $\gamma>0$ captures the weight assigned to load shedding cost and
violations of capacity constraints. Line 9 to line 13 evaluates the
performance of the islanding strategy $(E,d^{\prime})$ using function
$\bar{F}(E,\Lambda)$ with $\Lambda=\\{d^{\prime}\\}$. The algorithm finally
decides whether the islanding strategy should be updated (if the condition in
line 9 holds) or not (line 12).

We conclude this section by presenting the optimality guarantee provided by
Algorithm 1.

######  Theorem 2.

Let $Z$ be a sufficiently large positive number such that
$Z-\bar{F}(E,\Lambda)\geq 0$ for all $(E,\Lambda)$ where $|\Lambda|\leq 1$. As
parameter $\epsilon\rightarrow 0$, Algorithm 1 returns an islanding strategy
$(E,d)$ such that

$Z-\bar{F}(E,\Lambda)\geq\frac{1}{2}(Z-\bar{F}(E^{\prime},\Lambda^{\prime}))$

for all $(E^{\prime},\Lambda^{\prime})\in\mathcal{B}(\bar{\mathcal{M}})$,
where $\Lambda=\\{d\\}$ and $|\Lambda^{\prime}|\leq 1$.

###### Proof.

By Theorem 1, we have that $\bar{F}(E,\Lambda)$ is hybrid supermodular.
Therefore, $Z-\bar{F}(E,\Lambda)$ is hybrid submodular. By applying
Proposition 1 yields the theorem. ∎

Denote $E_{t}$ and $d_{t}$ as the set of edges and load shedding obtained
after $t$ iterations by using Algorithm 1. We have
$\overline{F}(E_{t},d_{t})<(1-\epsilon)^{t}\overline{F}(E_{0},d_{0})$. Hence
we have that the algorithm terminates within
$\lceil\frac{\log{\left\\{\frac{\overline{F}(E_{t})}{\overline{F}(E_{0})}\right\\}}}{\log{(1-\epsilon)}}\rceil$
iterations and each iteration has worst case $O(|\mathcal{E}|^{2}\xi)$
complexity, yielding an upper bound of

$O\left(\lceil\frac{\log{\left\\{\frac{\overline{F}(E_{t})}{\overline{F}(E_{0})}\right\\}}}{\log{(1-\epsilon)}}|\mathcal{E}|^{2}\xi\rceil\right)$

on the complexity, where $\xi$ is the computational complexity of solving the
linear program in Eqn. (18).

##  VI Numerical Study

This section presents three case studies to compare our proposed approach and
the state-of-the-art MILP-based solution. We first present necessary
background on the MILP formulation. We then present a comparison of results
from the MILP baseline and our approach on the IEEE 118-bus system [33], IEEE
300-bus system [34], ActivSg 500-bus system [35], and Polish 2383-bus system
[36]. We consider that at time $t=0$, buses 10, 69, 18, and 10 incur three-
phase faults for IEEE 118-bus, 300-bus, ActivSg 500-bus, and Polish 2383 test
cases, respectively, which necessitates controlled islanding. The fault is
cleared at time 0.4s.

###  VI-A Baseline: Mixed Integer Linear Program

This subsection introduces the baseline approach that utilizes mixed-integer
linear program (MILP) to solve the islanding problem. The baseline approach is
developed based on [10, 11].

The baseline approach takes the coherent generator groups, denoted as
$\mathcal{C}_{1},\ldots,\mathcal{C}_{m}$, as input, and computes the load
shedding at each load bus, the transmission lines to trip to partition the
power system, and the post-islanding power flow. We summarize the notations
and variables used in the MILP in Table I. The MILP given as below minimizes
the costs incurred by load shedding and power flow disruption.

$\displaystyle\min\quad$
$\displaystyle\alpha_{2}F_{2}(I_{1},\ldots,I_{m})+\alpha_{3}F_{3}(I_{1},\ldots,I_{m})$
(19a) s.t.
$\displaystyle\sum_{j\in\mathcal{L}}d_{kj}\leq\sum_{j\in\mathcal{G}}\bar{g}_{j}x_{kj},~{}\forall
k=1,\ldots,m,~{}\forall j\in\mathcal{V}$ (19b) $\displaystyle 0\leq
d_{kj}\leq\bar{d}_{j}x_{kj},~{}\forall k=1,\ldots,m,~{}\forall
j\in\mathcal{V}$ (19c) $\displaystyle x_{kj}=v_{kj^{\prime}},~{}\forall
j,j^{\prime}\in\mathcal{C}_{k}$ (19d) $\displaystyle
w_{k,j,j^{\prime}}\in\\{0,1\\},~{}\forall
k=1,\ldots,m,\forall(j,j^{\prime})\in\mathcal{E}$ (19e) $\displaystyle
x_{kj}\in\\{0,1\\},~{}\forall k=1,\ldots,m,\forall j\in\mathcal{V}$ (19f)
$\displaystyle
z_{jj^{\prime}}\in\\{0,1\\},~{}\forall(j,j^{\prime})\in\mathcal{E}$ (19g)
$\displaystyle\sum_{k=1}^{m}x_{kj}\leq 1,~{}\forall j\in\mathcal{V}$ (19h)
$\displaystyle w_{k,j,j^{\prime}}\leq x_{kj},~{}w_{k,j,j^{\prime}}\leq
x_{kj^{\prime}},~{}\forall k,(j,j^{\prime})\in\mathcal{E}$ (19i)
$\displaystyle
z_{jj^{\prime}}=\sum_{k=1}^{m}w_{k,j,j^{\prime}},z_{jj^{\prime}}=z_{j^{\prime}j},~{}\forall(j,j^{\prime})\in\mathcal{E}$
(19j) $\displaystyle 0\leq l_{k,j,j^{\prime}}\leq
Zz_{jj^{\prime}},~{}\forall(j,j^{\prime})\in\mathcal{E}$ (19k) $\displaystyle
v_{kj}\sum_{j\in\mathcal{V}}x_{kj}-x_{kj}+\sum_{j\in\mathcal{N}(j^{\prime})}f_{k,j,j^{\prime}}=\sum_{j\in\mathcal{N}(j^{\prime})}f_{k,j^{\prime},j},$
$\displaystyle\quad\quad\quad\quad\quad\quad\forall
j\in\mathcal{V},k=1,\ldots,m$ (19l)
$\displaystyle\tilde{P}_{jj^{\prime}}\leq\bar{P}_{jj^{\prime}}z_{jj^{\prime}},~{}\forall(j,j^{\prime})\in\mathcal{E}$
(19m)
$\displaystyle\sum_{j\in\mathcal{N}(j^{\prime})}\tilde{P}_{jj^{\prime}}+g_{j^{\prime}}-d_{j^{\prime}}=0,~{}\forall
j^{\prime}\in\mathcal{V}$ (19n) $\displaystyle\sum_{j}x_{kj}\mathbb{I}_{j}\geq
1,~{}\forall k=1,\ldots,m$ (19o)

In order to retain linearity, we note from Eqn. (19a) that the MILP cannot
optimize over generator coherency in Eqn. (1) as our proposed approach does.
In the MILP, function $F_{3}(I_{1},\ldots,I_{m})$ can be represented as
$F_{3}(I_{1},\ldots,I_{m})=\sum_{(j,j^{\prime})\in\mathcal{E}}(1-z_{jj^{\prime}})\frac{|P_{jj^{\prime}}|+|P_{j^{\prime}j}|}{2}$.
Constraint (19b) and (19c) require the admissible loads at the each load bus
to satisfy the load-generation balance. Constraint (19d) specifies the
reference generator of each island $I_{k}$. Constraints (19e) to (19j) ensure
that variables $x$, $w$, and $z$ define a set of disjoint islands. Constrains
(19k) and (19l) define an auxiliary flow $l_{k,j,j^{\prime}}$ initiating from
the reference generator to ensure the flow conservation law and connectivity
within each island. Constraints (19m) and (19n) define the post-islanding
power flow $\tilde{P}_{jj^{\prime}}$ on each transmission line. Constraint
(19o) assigns the blackstart generators to ensure efficient restoration. The
MILP formulated in Eqn. (19) involves
$3m|\mathcal{V}|+2m|\mathcal{E}|+2|\mathcal{E}|$ number of decision variables
and $(4m+2)|\mathcal{V}|+(3m+7)|\mathcal{E}|+m$ number of constraints.

TABLE I: This table contains the set of variables used in the MILP formulation. The first column gives the notation of each variable. The last column presents the interpretation of each variable. Notation | Type | Interpretation
---|---|---
$x_{kj}$ | Binary | Whether bus $j$ is contained in island $I_{k}$ ($x_{kj}=1$) or not ($x_{kj}=0$)
$w_{k,j,j^{\prime}}$ | Binary | Whether transmission line $(j,j^{\prime})\in\mathcal{E}$ is contained in island $I_{k}$ ($z_{k,j,j^{\prime}}=1$) or not ($z_{k,j,j^{\prime}}=0$)
$z_{jj^{\prime}}$ | Binary | Whether transmission line $(j,j^{\prime})$ will be tripped ($z_{jj^{\prime}}=0$) or not ($z_{jj^{\prime}}=1$)
$v_{kj}$ | Binary | Indicator for the reference generator $j$ in island $I_{k}$
$l_{k,j,j^{\prime}}$ | Nonnegative real | Auxiliary flow on transmission line $(j,j^{\prime})$ when included in island $I_{k}$
$Z$ | Positive real | A sufficiently large positive number
$\mathcal{N}(j)$ | Set | The neighboring buses $\\{j^{\prime}:(j^{\prime},j)\in\mathcal{E}\\}$ of $j$

We solve the MILP formulated in (19) by using two solvers: intlinprog provided
by Matlab [37] and mixed-integer program solver provided by Gurobi [38]. We
refer to the former as MILP-M, while the latter as MILP-G in the rest of this
section. Since MILP-G and MILP-M cannot optimize the generator coherency, we
compute the generator coherency for MILP-M and MILP-G beforehand for the
purpose of comparison.

(a) IEEE 300-bus, $m=3$ using MILP.

(b) IEEE 300-bus, $m=3$ using Algorithm 1.

(c) ActivSg 500-bus, $m=3$ using MILP.

(d) ActivSg 500-bus, $m=3$ using Algorithm 1.

Figure 1: Comparison of islanding strategies given by MILP-based approach and
Algorithm 1 for IEEE 300-bus and ActivSg 500-bus test cases with $m=3$ islands
(in blue, green, and purple colors). Each node in the graph represents a bus
of power system. The edges correspond to the set of transmission lines. The
set of transmission lines tripped in each case is marked in red, whereas the
those remain in the power system are in grey color. Isolated buses are marked
in black color.

###  VI-B Simulation Results

In our case study, we choose the load shedding cost function $c_{j}(x)$ as
$c_{j}(x)=\zeta_{j}x_{j}$, where $\zeta_{j}$ is randomly generated following a
uniform distribution within range $(1,100)$. Trade-off parameters
$\alpha_{1},\ldots,\alpha_{6}$ are chosen to be $1$. We choose parameter
$\epsilon=0$. We compare our proposed approach in Algorithm 1 with the
baseline approach in two settings for the evaluation purpose, where the first
setting partitions the test cases into $m=2$ islands, and the second setting
generates $m=3$ islands.

TABLE II: This table presents the set of reference generators given to
Algorithm 1 and the number of transmission lines tripped by islanding
strategies obtained using MILP-G (MILP using Gurobi), MILP-M (MILP using
intlinprog), and Algorithm 1 (hybrid submodular approach). The islanding
strategy given by Algorithm 1 trips fewer transmission lines in all test cases
including IEEE 118-bus, IEEE 300-bus, ActivSg 500-bus, and Polish 2383-bus
systems with different settings, and thus the proposed approach is more
practical for implementation.

Test Case |  Number of Islands |  Ref. Generators (for Algorithm 1) | Number of Transmission Lines Tripped
---|---|---|---
MILP-G |  MILP-M | Algorithm 1
IEEE 118-bus | 2 | 12, 100 | 9 | 9 | $\mathbf{5}$
3 | 26, 65, 80 | 14 | 14 | $\mathbf{8}$
IEEE 300-bus | 2 | 91, 198 | 14 | 5 | $\mathbf{4}$
3 | 92, 124, 198 | 18 | 20 | $\mathbf{12}$
ActivSg 500-bus | 2 | 9, 16 | 34 | 42 | $\mathbf{11}$
3 | 9, 82, 231 | 45 | 32 | $\mathbf{6}$
Polish 2383-bus | 2 | 41, 1726 | NA | NA | $\mathbf{3}$
3 | 45, 125, 1106 | NA | NA | $\mathbf{11}$

Table II presents the numbers of transmission lines removed by the baseline
and our approach. We observe that the islanding strategies given by Algorithm
1 trip fewer transmission lines compared to the baseline in all test cases,
and hence are more practical to be implemented. Fig. 1 compares the islanding
strategies given by MILP and Algorithm 1 for IEEE 300-bus and ActivSg 500-bus
test cases when the desired number of islands is set to $m=3$.

We summarize the amount of load shedding, the values of objective function
$\bar{F}(E,\Lambda)$, and run time of MILP-M, MILP-G, and Algorithm 1 in Table
III. We note that for Polish 2383-bus system, MILP-M and MILP-G run out of
memory (“NA” entries in Table III), and thus do not give an islanding
solution. We observe that our approach outperforms MILP-M and MILP-G in terms
of the total cost $\bar{F}(E,\Lambda)$ for all test cases. The amount of load
shedding required by Algorithm 1 is no larger than MILP-M or MILP-G for all
test cases. From the last column of Table III, we note that the run time of
the MILP-M grows fast as the scales of test cases increase, making the
approach not applicable for Polish 2383-bus system. Although MILP-G takes the
least run time for small-scale (IEEE 118-bus and IEEE 300-bus) and medium-
scale (ActivSg 500-bus) test cases, it is not applicable to large-scale test
case (Polish 2383-bus). Our proposed approach computes a controlled islanding
solution for large-scale test cases within a reasonable amount of time. We
further remark that neither MILP-M nor MILP-G could optimize over generator
coherency and thus involve fewer decision variables compared with our
approach.

For the ActivSg 500-bus test case with $m=2$ islands, the run time of MILP-G
and MILP-M are about $12$ times less than our approach. However, in this
scenario the islanding strategies given by MILP-G and MILP-M trip about three
to four times as many transmission lines compared to Algorithm 1 (see Table
II). Furthermore, the islanding strategies given by MILP-G and MILP-M
disconnect some buses from the islands, whereas our approach ensures that each
bus is connected within one island. We note that ActivSg 500-bus system
contains a collection of star subgraphs. To improve the run time of Algorithm
1, we can view each star subgraph as one node, and implement Algorithm 1 on
this reduced graph to approximate the islanding strategy.

Finally, we evaluate how the choices of $\alpha_{1},\alpha_{2},\alpha_{3}$
impact the controlled islanding strategy. Our choices of
$\alpha_{1},\alpha_{2},\alpha_{3}$ assign higher weight to the generator
coherency. The controlled islanding strategy did not change as we increase the
values of $\alpha_{2}$ and $\alpha_{3}$. However, by further increasing
$\alpha_{1}$ to be $\alpha_{1}>3$, we observe that the controlled islanding
strategy tripped more transmission lines to improve generator coherency.

TABLE III: This table summarizes the amount of load shedding, value of the
objective function $\bar{F}(E,\Lambda)$, and run time for IEEE 118-bus, IEEE
300-bus, ActivSg 500-bus, and Polish 2383-bus systems by using MILP-G (MILP
using Gurobi), MILP-M (MILP using intlinprog), and Algorithm 1 (hybrid
submodular approach) to generate $m\in\\{2,3\\}$ islands. “NA” in the table
represents that the baseline approaches run out of memory for Polish 2383-bus
system.

Test Case |  Number of Islands | Amount of Load Shedded | Objective Function Value | Run Time
---|---|---|---|---
MILP-G | MILP-M | Algorithm 1 | MILP-G | MILP-M | Algorithm 1 | MILP-G | MILP-M | Algorithm 1
IEEE 118-bus | 2 | 1691.00 MW | 1691.00 MW | 1526.00 MW | $4.9942\times 10^{8}$ | $4.9942\times 10^{8}$ | $\mathbf{4.9941\times 10^{8}}$ | $\mathbf{2.35}$ | 52.96 s | 101.47 s
3 | 1562.00 MW | 1562.00 MW | 1497.00 MW | $2.3598\times 10^{7}$ | $2.3598\times 10^{7}$ | $\mathbf{2.3587\times 10^{7}}$ |  $\mathbf{2.86}$ s | 35.48 s | 154.99 s
IEEE 300-bus | 2 | 16410.95 MW | 16267.45 MW | 16267.45 MW | $6.446\times 10^{5}$ | $6.3456\times 10^{5}$ | $\mathbf{6.3431\times 10^{5}}$ | $\mathbf{7.94s}$ | 7205.2 s |  $92.03$ s
3 | 16571.95 MW | 16571.95 MW | 16267.45 MW | $7.7788\times 10^{5}$ | $7.7847\times 10^{5}$ | $\mathbf{7.5603\times 10^{5}}$ | $\mathbf{35.28s}$ | 7205.9 s | 366.32 s
ActivSg 500-bus | 2 | 6277.30 MW | 6277.30 MW | 6277.30 MW | $3.1670\times 10^{8}$ | $3.1670\times 10^{8}$ | $\mathbf{3.1665\times 10^{8}}$ | 16.34 s |  $\mathbf{10.71}$ s | 130.82 s
3 | 6357.59 MW | 6357.59 MW | 6320.49 MW | $2.0354\times 10^{6}$ | $2.0284\times 10^{6}$ | $\mathbf{2.0152\times 10^{6}}$ |  $\mathbf{118.54}$ s | 7217.2 s | 208.78 s
Polish 2383-bus | 2 | NA | NA | 7511.15 MW | NA | NA | $\mathbf{1.7234\times 10^{6}}$ | NA | NA |  $\mathbf{335.80}$ s
3 | NA | NA | 7579.65 MW | NA | NA | $\mathbf{8.9825\times 10^{6}}$ | NA | NA |  $\mathbf{1328.41}$ s

##  VII Conclusion

In this paper, we investigated the problem of computing a controlled islanding
strategy for large-scale power systems. We formulated the problem by taking
power flow disruption, generator coherency, post-islanding stability, and
blackstart generator allocation into consideration. We translated the
formulated controlled islanding problem into a matroid optimization program.
We presented the concept of hybrid submodularity, and proved that the metrics
considered are hybrid supermodular. Based on this insight, we developed an
efficient local search algorithm with $\frac{1}{2}$-optimality guarantee. We
compared our solution approach with a baseline using mixed-integer linear
program on four test cases including IEEE 118-bus, IEEE 300-bus, ActivSg
500-bus, and Polish 2383-bus systems. The proposed solution approach found an
islanding strategy for each test case that could outperform the baseline
approach. Moreover, our proposed approach scaled well to large cases such as
Polish 2383-bus system, while the baseline did not return a result due to the
problem size.

## References

  * [1] M. R. Lee, J. M. Assante, and T. Conway, “Analysis of the cyber attack on the Ukrainian power grid,” https://www.eisac.com/cartella/Asset/00006542/TLP_WHITE_E-ISAC_SANS_Ukraine_DUC_6_Modular_ICS_Malware\%20Final.pdf?parent=64412.
  * [2] J. W. Busby, K. Baker, M. D. Bazilian, A. Q. Gilbert, E. Grubert, V. Rai, J. D. Rhodes, S. Shidore, C. A. Smith, and M. E. Webber, “Cascading risks: Understanding the 2021 winter blackout in Texas,” _Energy Research & Social Science_, vol. 77, p. 102106, 2021.
  * [3] J. Yardley and G. Harris, “2nd day of power failures cripples wide swath of India,” https://www.nytimes.com/2012/08/01/world/asia/power-outages-hit-600-million-in-india.html.
  * [4] P. Hines, K. Balasubramaniam, and E. C. Sanchez, “Cascading failures in power grids,” _IEEE Potentials_ , vol. 28, no. 5, pp. 24–30, 2009.
  * [5] S. Masood, “Power outage sweeps Pakistan, dropping millions into darkness,” The New York Times, https://www.nytimes.com/2023/01/23/world/asia/pakistan-power-outage-blackouts.html.
  * [6] S. Kamali, T. Amraee, and M. Fotuhi-Firuzabad, “Controlled islanding for enhancing grid resilience against power system blackout,” _IEEE Transactions on Power Delivery_ , vol. 36, no. 4, pp. 2386–2396, 2020.
  * [7] S. Yusof, G. Rogers, and R. Alden, “Slow coherency based network partitioning including load buses,” _IEEE Transactions on Power Systems_ , vol. 8, no. 3, pp. 1375–1382, 1993.
  * [8] H. You, V. Vittal, and X. Wang, “Slow coherency-based islanding,” _IEEE Transactions on Power Systems_ , vol. 19, no. 1, pp. 483–491, 2004.
  * [9] G. Xu and V. Vittal, “Slow coherency based cutset determination algorithm for large power systems,” _IEEE Transactions on Power Systems_ , vol. 25, no. 2, pp. 877–884, 2009.
  * [10] A. Kyriacou, P. Demetriou, C. Panayiotou, and E. Kyriakides, “Controlled islanding solution for large-scale power systems,” _IEEE Transactions on Power Systems_ , vol. 33, no. 2, pp. 1591–1602, 2017.
  * [11] G. Patsakis, D. Rajan, I. Aravena, and S. Oren, “Strong mixed-integer formulations for power system islanding and restoration,” _IEEE Transactions on Power Systems_ , vol. 34, no. 6, pp. 4880–4888, 2019.
  * [12] F. Teymouri and T. Amraee, “An MILP formulation for controlled islanding coordinated with under frequency load shedding plan,” _Electric Power Systems Research_ , vol. 171, pp. 116–126, 2019.
  * [13] P. A. Trodden, W. A. Bukhsh, A. Grothey, and K. I. McKinnon, “Optimization-based islanding of power networks using piecewise linear AC power flow,” _IEEE Transactions on Power Systems_ , vol. 29, no. 3, pp. 1212–1220, 2013.
  * [14] P. Trodden, W. Bukhsh, A. Grothey, and K. McKinnon, “MILP formulation for controlled islanding of power networks,” _International Journal of Electrical Power & Energy Systems_, vol. 45, no. 1, pp. 501–508, 2013.
  * [15] S. Kamali, T. Amraee, and M. Khorsand, “Intentional power system islanding under cascading outages using energy function method,” _IET Generation, Transmission & Distribution_, vol. 14, no. 20, pp. 4553–4562, 2020.
  * [16] M. Esmaili, M. Ghamsari-Yazdel, N. Amjady, and C. Chung, “Convex model for controlled islanding in transmission expansion planning to improve frequency stability,” _IEEE Transactions on Power Systems_ , vol. 36, no. 1, pp. 58–67, 2020.
  * [17] S. A. N. Sarmadi, A. S. Dobakhshari, S. Azizi, and A. M. Ranjbar, “A sectionalizing method in power system restoration based on WAMS,” _IEEE Transactions on Smart Grid_ , vol. 2, no. 1, pp. 190–197, 2011.
  * [18] J. Q. Tortós and V. Terzija, “Controlled islanding strategy considering power system restoration constraints,” in _IEEE Power and Energy Society General Meeting_. IEEE, 2012, pp. 1–8.
  * [19] Q. Zhao, K. Sun, D.-Z. Zheng, J. Ma, and Q. Lu, “A study of system splitting strategies for island operation of power system: A two-phase method based on OBDDs,” _IEEE Transactions on Power Systems_ , vol. 18, no. 4, pp. 1556–1565, 2003.
  * [20] K. Sun, D.-Z. Zheng, and Q. Lu, “Splitting strategies for islanding operation of large-scale power systems using OBDD-based methods,” _IEEE Transactions on Power Systems_ , vol. 18, no. 2, pp. 912–923, 2003.
  * [21] C. G. Wang, B. H. Zhang, Z. G. Hao, J. Shu, P. Li, and Z. Q. Bo, “A novel real-time searching method for power system splitting boundary,” _IEEE Transactions on Power Systems_ , vol. 25, no. 4, pp. 1902–1909, 2010.
  * [22] D. Sahabandu, L. Niu, A. Clark, and R. Poovendran, “A hybrid submodular optimization approach to controlled islanding with heterogeneous loads,” in _IEEE International Conference on Communications, Control, and Computing Technologies for Smart Grids (SmartGridComm)_ , 2022.
  * [23] L. Ding, F. M. Gonzalez-Longatt, P. Wall, and V. Terzija, “Two-step spectral clustering controlled islanding algorithm,” _IEEE Transactions on Power Systems_ , vol. 28, no. 1, pp. 75–84, 2012.
  * [24] K. Pang, C. Wang, N. D. Hatziargyriou, F. Wen, and Y. Xue, “Formulation of radiality constraintsfor optimal microgrid formation,” _IEEE Transactions on Power Systems_ , pp. 1–15, 2022.
  * [25] Z. Liu, A. Clark, L. Bushnell, D. S. Kirschen, and R. Poovendran, “Controlled islanding via weak submodularity,” _IEEE Transactions on Power Systems_ , vol. 34, no. 3, pp. 1858–1868, 2018.
  * [26] D. Sahabandu, L. Niu, A. Clark, and R. Poovendran, “A submodular optimization approach to stable and minimally disruptive controlled islanding in power systems,” in _American Control Conference (ACC)_. IEEE, 2022, pp. 4587–4594.
  * [27] J. Qin, I. Yang, and R. Rajagopal, “Submodularity of energy storage placement in power networks,” in _IEEE 55th Conference on Decision and Control (CDC)_. IEEE, 2016, pp. 686–693.
  * [28] Z. Liu, A. Clark, P. Lee, L. Bushnell, D. Kirschen, and R. Poovendran, “Submodular optimization for voltage control,” _IEEE Transactions on Power Systems_ , vol. 33, no. 1, pp. 502–513, 2017.
  * [29] A. Khodabakhsh, G. Yang, S. Basu, E. Nikolova, M. C. Caramanis, T. Lianeas, and E. Pountourakis, “A submodular approach for electricity distribution network reconfiguration,” _arXiv preprint arXiv:1711.03517_ , 2017.
  * [30] S. Fujishige, _Submodular Functions and Optimization_. Elsevier, 2005.
  * [31] J. G. Oxley, _Matroid Theory_. Oxford University Press, 2006, vol. 3.
  * [32] J. H. Chow, _Time-Scale Modeling of Dynamic Networks with Applications to Power Systems_. Springer, 1982.
  * [33] “IEEE 118-bus system,” http://www.ee.washington.edu/research/pstca/.
  * [34] “IEEE 300-bus system,” http://www.ee.washington.edu/research/pstca/pf300/pg_tca300bus.htm.
  * [35] A. B. Birchfield, T. Xu, K. M. Gegner, K. S. Shetye, and T. J. Overbye, “Grid structural characteristics as validation criteria for synthetic networks,” _IEEE Transactions on Power Systems_ , vol. 32, no. 4, pp. 3258–3265, 2016\.
  * [36] R. D. Zimmerman, C. E. Murillo-Sánchez, and R. J. Thomas, “MATPOWER: Steady-state operations, planning, and analysis tools for power systems research and education,” _IEEE Transactions on Power Systems_ , vol. 26, no. 1, pp. 12–19, 2010.
  * [37] “Mixed-integer linear programming (MILP) - MATLAB intlinprog,” https://www.mathworks.com/help/optim/ug/intlinprog.html.
  * [38] Gurobi Optimization, LLC, “Gurobi Optimizer Reference Manual,” 2023. [Online]. Available: https://www.gurobi.com

In what follows, we introduce generator coherency and the procedure to compute
the coherence matrix $\bar{A}$. Our discussion is based on [32].

We denote the rotor angle of generator $i$ at a steady state operating point
as $\delta_{i}$. Let $\Delta\delta_{i}$ be the rotor angle deviation of
generator $i$ from the steady operating point. Let
$\Delta\delta=[\Delta\delta_{1},\ldots,\Delta\delta_{|\mathcal{G}|}]^{T}$,
where $\Delta\delta_{i}$ is the deviation of rotor angle of generator $i$ from
the operating point. We assume that the generators follow the linearized swing
equation [32] given as

$\Delta\ddot{\delta}=M^{-1}K\Delta\delta,$

where $M$ is the inertia matrix, and $K$ is defined as

$K_{ij}=\begin{cases}-V_{i}V_{j}B_{ij}\cos(\delta_{i}-\delta_{j}),&\mbox{ if
}i\neq j\\\ -\sum_{k\neq i}^{n}K_{ij},&\mbox{ if }i=j\end{cases}.$ (20)

Here $V_{i}$ is the per unit voltage behind transient reactance of generator
$i$ and $B_{ij}$ is the imaginary part of the admittance.

Two generators $i$ and $j$ are said to be $\epsilon$-coherent if the maximum
difference between their rotor angles is bounded by $\epsilon$ for all time.
Let $\sigma$ be the set of eigenvalues of $M^{-1}K$. The $m$ number of
eigenvalues with smallest magnitudes then represent the slowest modes of
$M^{-1}K$. Let $U$ be the eigenbasis corresponding to these $m$ eigenvalues.
Let $U_{1}$ be the matrix obtained by extracting the rows of $U$ that
correspond to the reference generators. Then the coherence matrix $A$ is
computed as $\bar{A}=UU_{1}^{-1}$.

|  Luyao Niu (M’22) is a postdoctoral scholar at the Network Security Lab
(NSL), Department of Electrical and Computer Engineering, at the University of
Washington - Seattle. He received the B.Eng. degree from the School of
Electro-Mechanical Engineering, Xidian University, Xi’an, China, in 2013. He
received the M.Sc. degree and Ph.D. degree from the Department of Electrical
and Computer Engineering, Worcester Polytechnic Institute (WPI) in 2015 and
2022. He is the author of the GameSec Outstanding Paper (2018) and was
finalist for ACM/IEEE International Conference on Cyber-Physical Systems
(ICCPS) 2020 Best Paper Session Award. His research interests include
optimization, game theory, and scalable and verifiable control and security of
cyber physical systems.
---|---
|  Dinuka Sahabandu (M’23) is a postdoctoral scholar at the Network Security
Lab (NSL), Department of Electrical and Computer Engineering, at the
University of Washington - Seattle. He received the B.S. degree and M.S.
degree in Electrical Engineering from the Washington State University -
Pullman in 2013 and 2016, respectively. He received the Ph.D. degree from the
Department of Electrical and Computer Engineering at the University of
Washington - Seattle in 2023. His research interests include game theory for
network security and control of multi-agent systems.
---|---
|  Andrew Clark (M’15) is an Associate Professor of Electrical and Systems
Engineering at Washington University in St. Louis. He received the B.S. degree
in Electrical Engineering and the M.S. degree in Mathematics from the
University of Michigan - Ann Arbor in 2007 and 2008, respectively. He received
the Ph.D. degree from the Network Security Lab (NSL), Department of Electrical
Engineering, at the University of Washington - Seattle in 2014. He is author
or co-author of the IEEE/IFIP William C. Carter award- winning paper (2010),
the WiOpt Best Paper (2012), and the WiOpt Student Best Paper (2014), and was
a finalist for the IEEE CDC 2012 Best Student-Paper Award. He received the
GameSec Outstanding Paper Award (2018) and was finalist for the ACM ICCPS Best
Paper Award (2016, 2018, 2020). He won the General Motors AutoDriving Security
Award at VehicleSec (2023). He received the University of Washington Center
for Information Assurance and Cybersecurity (CIAC) Distinguished Research
Award (2012), Distinguished Dissertation Award (2014), an NSF CAREER award
(2020), and an AFOSR YIP award (2022). His research interests include control
and security of complex networks, submodular optimization, and control-
theoretic modeling of network security threats.
---|---
|  Radha Poovendran (F’15) is a Professor in the Department of Electrical and
Computer Engineering at the University of Washington (UW) - Seattle. He served
as the Chair of the Electrical and Computer Engineering Department at UW for
five years starting January 2015. He is the Director of the Network Security
Lab (NSL) at UW. He is the Associate Director of Research of the UW Center for
Excellence in Information Assurance Research and Education. He received the
B.S. degree in Electrical Engineering and the M.S. degree in Electrical and
Computer Engineering from the Indian Institute of Technology- Bombay and
University of Michigan - Ann Arbor in 1988 and 1992, respectively. He received
the Ph.D. degree in Electrical and Computer Engineering from the University of
Maryland - College Park in 1999. His research interests are in the areas of
wireless and sensor network security, control and security of cyber-physical
systems, adversarial modeling, smart connected communities, control-security,
games-security, and information theoretic security in the context of wireless
mobile networks. He is a Fellow of the IEEE for his contributions to security
in cyberphysical systems. He is a recipient of the NSA LUCITE Rising Star
Award (1999), National Science Foundation CAREER (2001), ARO YIP (2002), ONR
YIP (2004), and PECASE (2005) for his research contributions to multi-user
wireless security. He is also a recipient of the Outstanding Teaching Award
and Outstanding Research Advisor Award from UW EE (2002), Graduate Mentor
Award from Office of the Chancellor at University of California - San Diego
(2006), and the University of Maryland ECE Distinguished Alumni Award (2016).
He was co-author of award-winning papers including IEEE/IFIP William C. Carter
Award Paper (2010) and WiOpt Best Paper Award (2012).
---|---