# Modal Decomposition of the Linear Swing Equation
in Networks with Symmetries

Kshitij Bhatta , Majeed Hayat and Francesco Sorrentino This work was supported
by the National Science Foundation through Grant No. 1727948 and Grant
CRISP-1541148.K. Bhatta, is with the Department of Mechanical Engineering,
Univeristy of New Mexico, Albuquerque, NM 87131 USAM. Hayat, is with the
Department of Electrical and Computer Engineering, Marquette University,
Milwakee, WI 53233 USAF. Sorrentino is with the Department of Mechanical
Engineering, University of New Mexico, Albuquerquem NM 87131 USA (e-mail:
fsorrent@unm.edu).

(August 2020)

###### Abstract

Symmetries are widespread in physical, technological, biological, and social
systems and networks, including power grids. The swing equation is a classic
model for the dynamics of powergrid networks. The main goal of this paper is
to explain how network symmetries affect the swing equation transient and
steady state dynamics. We introduce a modal decomposition that allows us to
study transient effects, such as the presence of overshoots in the system
response. This modal decomposition provides insight into the peak flows within
the network lines and allows a rigorous characterization of the effects of
symmetries in the network topology on the dynamics. Our work applies to both
cases of homogeneous and heterogeneous parameters. Further, the model is used
to show how small perturbations propagate in networks with symmetries.
Finally, we present an application of our approach to a large power grid
network that displays symmetries.

##  I Introduction

Many papers have investigated network models which describe the dynamics of
power grids [27, 13, 14, 49, 38, 39, 1, 21, 43, 10, 17, 16, 20, 33, 34, 26].
Simplified models for the propagation of cascading failures on networks have
been proposed in [27, 13, 14]. More realistic models for the propagation of
cascading failures are based on the swing equation [38] or DC power flows [9,
34]. Symmetries play a significant role in the study of networked systems.
References [32, 18, 36, 8, 15, 44, 30, 51, 12, 5, 45, 46] have proposed tools
based on graph theory and group theory to analyze the dynamics of complex
networks with symmetries. The presence of symmetries in power grid networks
has been documented in [46, 15]. Reference [23] has analyzed how network
symmetries may affect synchronization modes of power grids and suggested that
symmetries may enhance the emergence of complete synchronization.

Despite this previous work, it appears that the effects of the network
symmetries on the dynamics of the swing equation have not been fully
elucidated. In what follows we first provide a definition for network
symmetries in the context of the swing equation, then we show how a reduced
representation of the dynamics based on the so-called ‘quotient network’ can
be achieved both in the cases of homogeneous and heterogeneous parameters.
Only a subset of the modes of the original network are _inherited_ by the
quotient network. Neglecting the remaining modes leads to an approximation,
which nonetheless can be quantified by considering an appropriately defined
error dynamics and applying a modal analysis to the error dynamics.

An important application of the study of network dynamics is network design.
Knowledge of how the dynamics changes in response to structural and dynamical
perturbations is key for designing resilient complex systems. A complex
network approach to study vulnerabilities of power grids has been presented in
[2]. Further, several articles have studied the existence of vulnerabilities
inherent to the power network structure [35, 11, 3, 53, 7, 52, 29] and the
importance of the transient dynamics in the propagation of failures has been
emphasized in [40].

Here we propose a mathematical analysis of the classic swing equation based on
a simplified description in terms of a network of coupled forced second order
systems. Our analysis provides immediate understanding of the swing equation
transient dynamics via a modal decomposition. Further, it allows us to
rigorously address the presence of network symmetries and their effects on the
dynamics.

The rest of this paper is organized as follows. Section II introduces the
swing equation for a generic network and presents its modal decomposition. The
effects of network symmetries are discussed in Sec. III. In the presence of
symmetries, a lower dimensional dynamical representation based on the so-
called quotient network is possible, either in the cases that the powers at
the network nodes respect the symmetries or not. In the latter case, the
quotient description only provides an approximation of the full transient
response. The deviation between the full network response and the quotient
network response can be characterized in terms of a properly defined error
dynamics. Effects of heterogeneity in the damping terms are presented in Sec.
IV. Finally, the conclusions are given in Sec. V.

##  II The Swing Equation and its Modal Decomposition

The swing equation is a classic model for the dynamics of power grid networks.
A power grid can be represented using an undirected graph such that each node
can be either a generator (generating power) or a load (consuming power). Each
node is regarded as a rotating machine (oscillator) and the presence of an
edge between two nodes corresponds to the presence of a transmission line
connecting them. The network connectivity is given by the symmetric matrix
$\tilde{A}=\\{\tilde{A}_{ij}\\}$, $\tilde{A}_{ij}=\tilde{A}_{ji}>0$ if nodes
$i$ and $j$ are connected, $\tilde{A}_{ij}=\tilde{A}_{ji}=0$ otherwise.

######  Definition 1.

An undirected network is defined by the set of nodes
$\mathcal{V}=(1,2,...,n)$, $|\mathcal{V}|=n$ and the set of edges or lines
$\mathcal{E}=\mathcal{V}\times\mathcal{V},$ such that $(i,j)\in\mathcal{E}$,
$i\in\mathcal{V}$, $j\in\mathcal{V}$, if $\tilde{A}_{ij}=\tilde{A}_{ji}>0$.

The state of each oscillator $i=1,...,n$ is characterized by its node
displacement $\theta_{i}$ and the nodal velocity $\omega_{i}=\dot{\theta_{i}}$
relative to a reference frequency $\omega_{R}$. $H_{i}$ is the inertia
constant of oscillator $i$, $D_{i}$ is its damping constant and $b_{i}$ is the
power generated (consumed) by the node. If we assume that
$\omega_{R}\approx\omega$ and lossless transmission lines, we can write the
swing equation [31, 22],

$\frac{2H_{i}}{\omega_{R}}\ddot{\vartheta_{i}}+\frac{D_{i}}{\omega_{R}}\dot{\vartheta_{i}}={b_{i}}-\sum_{{j=1,j\neq
i}}^{n}{\tilde{A}_{ij}}\sin(\vartheta_{i}-\vartheta_{j}),\hskip
8.53581pti=1,2,...n,$ (1)

$b_{i}>0$ for generators and $b_{i}<0$ for loads. Dividing both sides of Eq. 1
by $\frac{2H_{i}}{\omega_{R}}$, we get

$\ddot{\vartheta_{i}}+\frac{D_{i}}{2H_{i}}\dot{\vartheta_{i}}=\frac{\omega_{R}}{2H_{i}}{b_{i}}-\frac{\omega_{R}}{2H_{i}}\sum_{{j=1,j\neq
i}}^{n}{\tilde{A}_{ij}}\sin(\vartheta_{i}-\vartheta_{j}),$ (2)

$i=1,2,...n.$

If we define, $\frac{D_{i}}{2H_{i}}=\gamma_{i}$,
$J_{i}=\frac{\omega_{R}}{2H_{i}}$, ${p}_{i}=J_{i}b_{i}$,
${A}_{ij}=J_{i}\tilde{A}_{ij}$, Eq. (2) becomes

$\ddot{\vartheta_{i}}=-\gamma_{i}\dot{\vartheta_{i}}+p_{i}+\sum_{{j=1,j\neq
i}}^{n}{A}_{ij}\sin(\vartheta_{j}-\vartheta_{i}),$ (3)

$i=1,2,...n$. The flow in an edge/line $(i,j)$, with coupling ${A}_{ij}$ at
time $t$ is given by [40]:

$F_{ij}(t)={A}_{ij}\sin(\vartheta_{j}(t)-\vartheta_{i}(t)).$ (4)

######  Definition 2.

We say that the network is balanced if $\sum_{i=1}^{n}b_{i}=0$.

Unless differently noted, we proceed under the assumption that the network is
balanced.

In any large network, including power-grid systems, there are redundancies in
the form of symmetries [32, 18, 36, 15, 44, 30, 51, 12, 5].

######  Definition 3.

A symmetry for the set of Eqs. (3) is a permutation matrix $P$ such that
$P\boldsymbol{\gamma}=\boldsymbol{\gamma}$ and $P{A}={A}P$. The automorphism
group $\mathcal{G}$ is the set of all symmetries with the operation
composition. The set of all symmetries in the group will only permute certain
subsets of nodes (the _orbits_ or _clusters_) among each other. The set of
nodes $\mathcal{V}$ is partitioned into $q$ disjoint subsets of nodes
$\\{S_{1},S_{2},...S_{q}\\}$, $\cup_{i=1}^{q}S_{i}=\mathcal{V}$, $S_{i}\cap
S_{j}=\emptyset$ for $i\neq j$, with $n=n_{i}\sum_{i=1}^{q}n_{i}$ where
$n_{i}=|S_{i}|$.

The nodes in each subset are mapped into each other by application of one or
more symmetries in $\mathcal{G}$; however, there is no symmetry in
$\mathcal{G}$ that will map into each other nodes in different subsets. We
refer to such subsets of nodes as ‘clusters’ or ‘orbits’ of the automorphism
group. For a review of graph automorphisms, see [46].

######  Remark 1.

Consider a permutation $P$ of the network nodes that satisfies $AP=PA$. Say
$v,w\in\mathcal{V}$ two network nodes, call $v^{\prime}$
$(w^{\prime})\in\mathcal{V}$ the network node $v$ $(w)$ gets mapped to by
application of the permutation $P$. It follows that
$A_{vw}=A_{v^{\prime}w^{\prime}}$ and $A_{wv}=A_{w^{\prime}v^{\prime}}$.

######  Lemma 1.

A flow-invariant ‘synchronous’ solution
$\vartheta_{i}^{*}(t)=\vartheta^{k}(t)$ for all $i\in S_{k}$, $k=1,...,q$, is
induced by the automorphism group $\mathcal{G}$.

###### Proof.

Assume $\vartheta_{i}(0)=\vartheta^{k}(0)$ and
$\dot{\vartheta}_{i}(0)=\dot{\vartheta}^{k}(0)$ for all $i\in S_{k}$,
$k=1,...,q$. It follows that $\ddot{\vartheta}_{i}(0)=\ddot{\vartheta}^{k}(0)$
for all $i\in S_{k}$, $k=1,...,q$, from which the assertion follows. ∎

Stability of the nonlinear swing equation (3) has been studied in [50] by
using a Lyapunov function approach. The conditions of [50] can be directly
applied to ensure convergence of the flow invariant solution on a stable fixed
point $\vartheta_{i}^{*}=\vartheta^{k}$, for all $i\in S_{k}$, $k=1,...,q$.

The linearized swing equation, which models the propagation of small
disturbances (e.g., affecting the initial condition or affecting the power
supplied/demanded at different nodes) [47], is obtained by linearizing Eqs.
(3) about the stable fixed point $\vartheta_{i}^{*}$, $i=1,...,n$,

$\ddot{\theta_{i}}=-\gamma_{i}\dot{\theta_{i}}+{p}_{i}+\sum_{{j=1}}^{n}{L_{ij}}\theta_{j},$
(5)

$i=1,2,...n$, where the Laplacian matrix $L=\\{L_{ij}\\}$ has entries
$L_{ij}={A}_{ij}\cos(\vartheta_{j}^{*}-\vartheta_{i}^{*})-\delta_{ij}\sum_{j}{A}_{ij}\cos(\vartheta_{j}^{*}-\vartheta_{i}^{*})$,
and $\delta_{ij}$ is the Kronecker delta.Each term $p_{i}$ on the right hand
side of Eq. (5) effectively represents a small power deviation. In the rest of
this paper we will often approximate the nonlinear swing equation (3) with the
linearized swing equation (5), under the assumption of small power deviations.

######  Definition 4.

A symmetry for the set of Eqs. (5) is a permutation matrix $P$ such that
$P\boldsymbol{\gamma}=\boldsymbol{\gamma}$ and $P{L}={L}P$.

######  Lemma 2.

The set of Eqs. (3) and the set of Eqs. (5) have the same set of symmetries.

###### Proof.

We break the proof in two parts. We first show that (i) $PA=AP$ implies
$P\hat{A}=\hat{A}P$, where the matrix $\hat{A}=\\{\hat{A}_{vw}\\}$ has entries
$\hat{A}_{vw}=A_{vw}\cos(\vartheta_{v}^{*}-\vartheta_{w}^{*})$. Then show that
(ii) $P\hat{A}=\hat{A}P$ implies $PL=LP$. Assume $P$ is a symmetry for the set
of Eqs. (3). From Remark 1, it follows that for all $v,w\in\mathcal{V}$,
$A_{vw}=A_{v^{\prime}w^{\prime}}$ and $A_{wv}=A_{w^{\prime}v^{\prime}}$ where
$v^{\prime}$ ($w^{\prime}$) is the node $v$ ($w$) is mapped to by $P$. We are
now going to show that also $\hat{A}_{vw}=\hat{A}_{v^{\prime}w^{\prime}}$ and
$\hat{A}_{wv}=\hat{A}_{w^{\prime}v^{\prime}}$. If $v$ gets mapped into
$v^{\prime}$ they belong to the same cluster, and so also $w$ and
$w^{\prime}$, hence $\vartheta^{*}_{v}=\vartheta^{*}_{v^{\prime}}$ and
$\vartheta^{*}_{w}=\vartheta^{*}_{w^{\prime}}$. Hence,
$\cos(\vartheta^{*}_{v}-\vartheta^{*}_{w})=\cos(\vartheta^{*}_{v^{\prime}}-\vartheta^{*}_{w^{\prime}})$,
which proves (i). To prove (ii) we just need to note that for two nodes $k$
and $l$ to be in the same cluster it must be necessarily verified that
$\sum_{j}A_{kj}=\sum_{j}A_{lj}$ and so also that
$\sum_{j}\hat{A}_{kj}=\sum_{j}\hat{A}_{lj}$, see also [36, 46]. ∎

The linear approximation yields this expression for the flow in a line
$(i,j)\in\mathcal{E}$:

$F_{ij}=A_{ij}(\theta_{j}(t)-\theta_{i}(t)).$ (6)

In the rest of this paper, except for Sec. IV, we introduce the following
assumption of homogeneous damping terms:

######  Assumption 1.

All the damping terms are the same $\gamma_{i}=\gamma$, $i=1,...,n.$

The regime in which all the $\gamma_{i}=\gamma$ is of interest because it has
been shown to lead to the best properties in terms of stability of the
synchronous solution [28]. So, unless specified otherwise, in Secs. II and
III, we set $\gamma_{i}=\gamma=0.9$. In Sec. IV we will remove assumption 1
and generalize our results to the case of arbitrary $\gamma$’s.

######  Lemma 3.

If the original matrix $\tilde{A}$ is symmetric, the spectrum of the matrix
$L$ is real.

###### Proof.

We note that
$\hat{A}_{ij}=A_{ij}\cos(\vartheta^{*}_{i}-\vartheta^{*}_{j})=J_{i}\bar{A}_{ij}$,
where the matrix $\bar{A}=\\{\bar{A}_{ij}\\}$ with entries
$\bar{A}_{ij}=\tilde{A}_{ij}\cos(\vartheta^{*}_{i}-\vartheta^{*}_{j})$ is
symmetric. By defining the symmetric matrix $\bar{L}=\\{{\bar{L}}_{ij}\\}$
with entries
${\bar{L}}_{ij}={\bar{A}}_{ij}-\delta_{ij}\sum_{j}{\bar{A}}_{ij}$, we can
write $L=D\bar{L}$ and the matrix $D$ is a diagonal matrix with diagonal
entries $D_{ii}=J_{i}$. By construction, both the matrix $L$ and the matrix
$\bar{L}$ have sums of their rows equal to zero. However, the matrix $L$ is
generically asymmetric, while the matrix $\bar{L}$ is symmetric. We then write
the eigenvalue equation $L\mathbf{v}=D\bar{L}\mathbf{v}=\lambda\mathbf{v}$,
where $\mathbf{v}$ is an eigenvector and $\lambda$ is an eigenvalue for $L$.
Then, we premultiply the above equation by $D^{-1/2}$ and obtain:

$D^{1/2}\bar{L}\mathbf{v}=\lambda D^{-1/2}\mathbf{v}.$

Now, by setting $\mathbf{w}=D^{-1/2}\mathbf{v}$ we can write

$D^{1/2}\bar{L}D^{1/2}\mathbf{w}=\lambda\mathbf{w},$

where the matrix $L^{\prime}=D^{1/2}\bar{L}D^{1/2}$ is symmetric, and thus the
eigenvalues $\lambda$ are real.
To conclude, $D\bar{L}$ is a special case of a generically asymmetric matrix
with real spectrum. However, under generic conditions, the real eigenvalues of
$L$ and $\bar{L}$ are not the same. ∎

Since the sums of the entries along all the rows of the matrix $L$ (of the
matrix $\tilde{L}$) are zero, it follows that this matrix has at least one
eigenvalue equal zero. The multiplicity of the zero eigenvalue is equal to the
number of connected components of the network. We now distinguish the two
cases that the network is (i) connected and (ii) not connected. In case (i),
the Laplacian matrix has a single zero eigenvalue $\lambda_{1}=0$ with
associated eigenvector $[1,1,...,1]$ and the remaining eigenvalues
$\lambda_{2}\geq\lambda_{3}\geq....\lambda_{n}$ are all negative. In case
(ii), assume there are $c$ connected components. Then, after appropriate
labeling of the nodes, the Laplacian matrix can be written as

$L=\begin{bmatrix}L_{1}&&\\\ &\ddots&\\\ &&L_{c}\end{bmatrix},$ (7)

where $L_{i}$, $i=1,...,c$ is the Laplacian matrix associated to connected
component $i$. It follows that the multiplicity of the zero eigenvalue of the
Laplacian matrix $L$ is equal to $c$, $\lambda_{1}=...=\lambda_{c}=0$, and the
corresponding eigenvectors are each one associated with a connected component
$i=1,...,c$; these $c$ eigenvectors are called the component vectors of the
network connected components. The remaining eigenvalues
$\lambda_{c+1}\geq\lambda_{c+2}\geq....\lambda_{n}$ are all negative.

######  Definition 5.

Component vectors are vectors that have entries $j$ equal to $1$ if node $j$
is in connected component $i$ and equal to $0$ otherwise.

######  Lemma 4.

Assume the network is connected. The left eigenvector of the matrix
$L=D\tilde{L}$ associated with the zero eigenvalue is equal to
$[D_{11}^{-1},D_{22}^{-1},...,D_{nn}^{-1}]$.

###### Proof.

By assumption, the matrix $L$ has only one zero eigenvalue. Consider the
eigenvalue equation $\textbf{z}^{T}D\tilde{L}=\textbf{0}$, with left
eigenvector z. The vector $\textbf{z}^{T}D$ is the left eigenvector of the
matrix $\tilde{L}$ associated with the zero eigenvalue and therefore, it is
equal to $[1,1,...,1].$ It follows that the entries of z are the reciprocal of
the entries on the main diagonal of the matrix $D$. ∎

We now consider the vector
${\boldsymbol{\theta}}=[\theta_{1},\theta_{2},...\theta_{n}]$ and the vector
${\mathbf{p}}=[{p}_{1},{p}_{2},...,{p}_{n}]$, and rewrite Eq. (5),

$\ddot{\boldsymbol{\theta}}(t)=-\gamma\dot{\boldsymbol{\theta}}(t)+L\boldsymbol{\theta}(t)+\mathbf{p}.$
(8)

We first assume the network is connected. We diagonalize $L$, $L=V\Lambda
V^{-1}$, where $\Lambda=(\lambda_{1},\lambda_{2},...,\lambda_{n})$ is the
matrix of the eigenvalues of $L$ and $V=\\{V_{ij}\\}$ is the matrix of the
eigenvectors of $L$. We multiply Eq. (3) on the left by $V^{-1}$ and by
calling $\boldsymbol{\eta}(t)=V^{-1}\boldsymbol{\theta}(t)$ and
$\mathbf{q}=V^{-1}\mathbf{p}$, we obtain,

$\ddot{{\boldsymbol{\eta}}}(t)=-\gamma\dot{{\boldsymbol{\eta}}}(t)+\Lambda{\boldsymbol{\eta}}(t)+{\mathbf{q}},$
(9)

which can be broken up into $n$ independent equations or ‘modes’,

$\ddot{{{\eta}}}_{i}(t)=-\gamma\dot{{{\eta}}}_{i}(t)+\lambda_{i}{{\eta}_{i}}(t)+{{q}_{i}},$
(10)

$i=1,...,n$. Modal decompositions similar to Eq. (10) have been also obtained
in [4, 19, 24, 47, 37].Tyloo et al. [47] relax the constant inertia to damping
ratio assumption and show that their derivation is still valid with
heterogeneous dynamical parameters. In what follows we study the effects of
network symmetries on the modal decomposition. We also show how the modes can
be exploited in order to compute maximum flows over lines and in designing
line capacities.

######  Lemma 5.

: Assume a balanced power grid, i.e., $\sum_{i}b_{i}=0$. Then, $q_{1}=0$.

###### Proof.

Consider $p_{i}=J_{i}b_{i}$ and $q_{1}=\mathbf{z^{T}}\mathbf{p}$, where
$\mathbf{z}$ is the left eigenvector of $L$ associated with $\lambda_{1}$.
From Lemma 4, we know that
$\textbf{z}^{T}=[\frac{1}{J_{1}},\frac{1}{J_{2}},...\frac{1}{J_{n}}]$, so we
get

$q_{1}=\sum_{i}\frac{1}{J_{i}}p_{i}\\\ $

and it follows that

$q_{1}=\sum_{i}\frac{1}{J_{i}}J_{i}b_{i}=\sum_{i}b_{i}=0.\\\ $

∎

Hence, for $i=1$, Eq. (10) becomes

$\ddot{{{\eta}}}_{1}(t)+\gamma\dot{{{\eta}}}_{1}(t)=0,$ (11)

from which we can see that for a large $t$, $\eta_{1}(t)$ approaches a
constant, which depends on the initial conditions.

For $i>1$, $\lambda_{i}<0$ ensures convergence of the corresponding
$\eta_{i}(t)$ in Eq. (10) for large $t$. Moreover, for $i=2,...,n$, Eq. (10)
can be written as,

$\ddot{{{\eta}_{i}}}(t)+2\zeta_{i}\omega_{i}\dot{{{\eta}_{i}}}(t)+{\omega_{i}}^{2}{{\eta}_{i}}(t)={{q}_{i}},$
(12)

where $-\lambda_{i}=\omega_{i}^{2}$ and $\zeta_{i}={\gamma}/{(2\omega_{i})}$,
$i=2,...,m$. Eq. (12) is the equation of an underdamped ($0<\zeta_{i}<1$,
$i=2,...,m$) or overdamped ($\zeta_{i}>1$, $i=m+1,...,n$) second order system
forced by a step function of amplitude $q_{i}$. Overdamped systems $i>m$ do
not give rise to overshoots. The solution to (12) can be written as
$\eta_{i}(t)=\bar{\eta}_{i}(t)+\hat{\eta}_{i}(t)$, where $\bar{\eta}_{i}(t)$
is the free evolution and $\hat{\eta}_{i}(t)$ is the forced evolution. The
free evolution decays exponentially in time. Thus, for large $t$, we can
assume that $\eta_{i}(t)\simeq\hat{\eta}_{i}(t)$.
For overdamped systems, the forced solution is equal to,

$\eta_{i}(t)=\frac{q_{i}}{\omega_{i}^{2}}\Bigr{[}1+\frac{\chi_{1}e^{{-\omega_{i}\chi_{2}t}}-\chi_{2}e^{{-\omega_{i}\chi_{1}t}}}{2\sqrt{\zeta_{i}^{2}-1}}\Bigr{]}$
(13)

where
$\quad\Bigr{(}\zeta_{i}-\sqrt{\zeta_{i}^{2}-1}\Bigr{)}=\chi_{1},\quad\Bigr{(}\zeta_{i}+\sqrt{\zeta_{i}^{2}-1}\Bigr{)}=\chi_{2}$.

For underdamped systems, the forced solution is equal to,

${\eta}_{i}(t)=\frac{q_{i}}{\omega_{i}^{2}}\Bigl{[}1-\frac{e^{-\zeta_{i}\omega_{i}t}}{\sqrt{1-\zeta_{i}^{2}}}\sin\Bigl{(}\omega_{i}\sqrt{1-\zeta_{i}^{2}}t-\cos^{-1}\zeta_{i}\Bigr{)}\Bigr{]}.$
(14)

In realistic systems with constant damping and inertia, it has been found that
all modes are underdamped and propagate through the whole system with
$\frac{D}{\lambda_{i}}<\gamma^{-1},\quad\forall i>1$ [47, 48]. Thus, for the
remainder of the paper, we will only focus on underdamped modes.

The underdamped system in (14) converges in steady state to

${\eta}_{i}^{ss}=\frac{q_{i}}{\omega_{i}^{2}},$ (15)

and is upper and lower bounded,

$\eta_{i}^{-}(t)\leq\eta_{i}(t)\leq\eta_{i}^{+}(t),$ (16)

where
$\eta_{i}^{\pm}(t)=({\omega_{i}^{2}}\sqrt{1-\zeta_{i}^{2}})^{-1}({q_{i}}\sqrt{1-\zeta_{i}^{2}}\pm{|q_{i}|}e^{-\zeta_{i}\omega_{i}t})$.

The peak time is equal to

$t_{i}=\frac{\pi}{\omega_{i}\sqrt{1-\zeta_{i}^{2}}},$ (17)

and the peak value is given by

${\eta}_{i}^{\mbox{peak}}=\frac{q_{i}}{\omega_{i}^{2}}\Bigl{[}1+\exp({-\pi\zeta_{i}}/{\sqrt{1-\zeta_{i}^{2}}})\Bigr{]}.$
(18)

The larger $q_{i}$ is, the higher is the peak. The smaller
$\omega_{i}^{2}=-\lambda_{i}$ is, the higher is the peak. Finally, the smaller
is $\zeta_{i}$, the higher is the peak. The eigenvalue $\lambda_{2}$ is
associated with the smallest $\omega_{2}$ and also the smallest $\zeta_{2}$,
and so is generically responsible for the largest peak.

As we know the $\eta_{i}(t)$’s, we can also compute the $\theta_{i}(t)$’s
using the formula $\boldsymbol{\theta}(t)=V\boldsymbol{\eta}(t)$, or
equivalently,

$\theta_{i}(t)=\sum_{j}V_{ij}\eta_{j}(t),$ (19)

It is important to note that the $V_{ij}$ can be either positive or negative,
and so are the $q_{i}$.

In general, the peak of $\theta_{i}(t)$ will be more strongly affected by
$\eta_{2}(t)$, followed by $\eta_{3}(t)$, $\eta_{4}(t)$, etc., but it will
also depend on the terms $V_{ij}q_{j}$. We can also compute the values of
$\theta_{i}$ at steady state,

$\theta_{i}^{ss}=\sum_{j}V_{ij}\eta_{j}^{ss}=\sum_{j}V_{ij}\frac{q_{j}}{\omega_{j}^{2}}.$
(20)

From Eq. (6), the flows are equal to

$F_{ij}(t)=\sum_{k=2}^{n}(V_{jk}-V_{ik})\eta_{k}(t),$ (21)

where the summation starts from $k=2$ since for a connected network
$V_{i1}=V_{j1}$ (the first column of the matrix $V$ is the eigenvector
associated with the eigenvalue $0$.) Each of the terms in the summation on the
right hand side of Eq. (21) can be either positive or negative. We rewrite the
right hand side of (21) as
$\sum_{k+}(V_{jk}-V_{ik})\eta_{k}(t)+\sum_{k-}(V_{jk}-V_{ik})\eta_{k}(t)$,
where the first summation is only over positive terms and the second summation
is only over negative terms. The absolute flow in a line is then upper bounded
by

$|F_{ij}(t)|\leq\sum_{k+}(V_{jk}-V_{ik})\eta_{k}^{+}(t)+\sum_{k-}|(V_{jk}-V_{ik})\eta_{k}^{-}(t)|.$
(22)

###  II-A Linear combinations of modes

Equation (19) (equation (21)) shows that the individual nodal displacements
$\theta_{i}(t)$ (the flows $F_{ij}(t)$) can be written as linear combinations
of the modes $\eta_{i}(t)$, $i=1,..,n$. This has immediate practical
implications. For example, because each mode is both upper bounded and lower
bounded, upper bounds on the absolute flows can be computed, see e.g., Eq.
(22). However, an open question is how peaks and peak times can be computed
for linear combinations of second order modes, such as those in Eqs. (19) and
(21). The reason why this is challenging is that peak times for different
$\eta$s are different, thus a linear combination of the peak $\eta$ values
will not return the peak flow. In this subsection we present a simple
numerical technique to approximate the peak flows, based only on minimal
knowledge about the individual modes.

Our original assumption is that the flows $F_{ij}$ are small, which is the
condition for approximating Eq. (3) with the linear swing equation (5).
However, flows may increase as a result of a variety of perturbations,
including internal power surges and outside attacks. It is also possible that
the absolute flow $|F_{ij}|$ over a line may increase above the so-called line
capacity, causing failure of that line and possibly lead to other line
failures in a cascade. A similar event would lead to violation of the small
flow assumption.

Reference [40] has pointed out the importance of the transient dynamic as it
is possible that an absolute flow may exceed a line capacity transiently. Our
analysis allows us to study the transient response of the swing equation
dynamics under the assumption of small flows. Under this condition, we can
compute the peak flows over each line, which can be used to design line
capacities that are robust against the effects of transient perturbations of
small entity.

Next we present a numerical technique to find the maximum absolute flow over a
line. For large $\gamma$, there are two possible ways in which a flow in a
line may evolve over time; _(1) The absolute flow reaches its maximum value at
the first peak time. (2) The absolute flow reaches its maximum value at steady
state_. Thus, to find the maximum flow accurately, one needs to calculate both
the first peak and the steady state value and pick the one which is largest in
magnitude.

In order to compute the peak times, we note that from Eq. (21) we have that

$\dot{F}_{ij}=C_{2}\dot{\eta}_{2}+C_{3}\dot{\eta}_{3}+...C_{n}\dot{\eta}_{n}$
(23)

where for underdamped modes

$\dot{\eta}_{k}=\frac{q_{k}}{\omega_{k}\sqrt{1-{\zeta_{k}}^{2}}}e^{-\zeta_{k}\omega_{k}t}\sin\omega_{k}\sqrt{1-{\zeta_{k}}^{2}}t,\hskip
8.53581ptk=m+1,...,n$ (24)

Since $\zeta_{k}\omega_{k}=\frac{\gamma}{2}$, the exponents in Eq. (24) are
the same for $k=m+1,...,n$. Further, say
$\varsigma_{k}=\omega_{k}\sqrt{1-{\zeta_{k}}^{2}}$. A necessary condition for
the flow $F_{ij}$ to achieve either a maximum or a minimum is

${\sum_{k=m+1}^{n}}\frac{C_{k}q_{k}}{\varsigma_{k}}\sin{\varsigma_{k}t}=0.$
(25)

Since the values of $\varsigma_{k}$ are different for different $k$, we are
unable to compute a closed form solution for (25). However, a root-finding
algorithm can be used to calculate the peak time and the maximum absolute flow
through the line. For example, that can be done by computing an initial guess
$\tau_{ij}$ and by iterating Newton’s method towards convergence on the peak
time.

For large enough $\gamma$ we expect the peak time to be the root of (25)
closest to the origin. Then a good choice of the initial guess for the root-
finding algorithm can be obtained by Taylor expanding (25) up to third order
about the origin and setting the third order expansion equal to zero. This
leads to the following initial guess,

$\tilde{\tau}_{ij}=\sqrt{\frac{6\sum_{k=m+1}^{n}C_{k}q_{k}}{\sum_{k=m+1}^{n}C_{k}q_{k}\varsigma_{k}^{2}}}.$
(26)

To demonstrate this, we consider the ‘bottleneck network’ shown in Fig. 1 with
$\gamma=1.5$ and $\textbf{p}=[0.3,-0.1,-0.1,0.1,-0.1,-0.2,0.1]^{T}$. We use
Eq. (26) to compute the initial guesses which are then used to approximate the
first peak times and subsequently the first peak values that are shown in
Table I. Table I also shows that for this example our approach based on the
linearized swing equation (5) well approximates the maximum flow obtained by
integration of the full nonlinear swing equation (3).

For non large $\gamma$s however, a slightly more cumbersome approach has to be
taken to find the initial guess. Peak $\eta$ values are first multiplied with
their respective flow coefficients and then these values are cumulatively
added in an ascending order of peak time. The peak time corresponding to the
largest of these cumulatively added values can be used as the initial guess
for the flow. It is important to note that the closer the peak times are, the
more accurate this initial guess is and the farther apart they are, the more
iterations will be required to converge to the actual solution.

We wish to emphasize that the approach described in this section can be
applied to any linear combination of second order modes. For example,we can
use this approach to approximate peak flows in quotient networks, Eq. (27) and
in error networks, Eq. (33).

TABLE I:  Steady state, first peak and maximum flow through lines for Fig. 1. The values for the second column are computed using Eq. (20) and the values for the third column are computed using Eq. (26), Eq. (14) and Eq. (21). The fourth and fifth columns are calculated by numerically solving the linear swing equation (5) and the non-linear swing equation (3), respectively. |  |  | Linear | Non-Linear
---|---|---|---|---
Line | Steady | First | Max | Max
| State | Peak | Flow | Flow
$\theta_{2}-\theta_{1}$ | -0.1750 | -0.1536 | -0.1750 | -0.1757
$\theta_{3}-\theta_{1}$ | -0.1250 | -0.1466 | -0.1466 | -0.1468
$\theta_{3}-\theta_{2}$ | 0.0500 | 0.0000 | 0.0500 | 0.0502
$\theta_{4}-\theta_{2}$ | 0.0750 | 0.0527 | 0.0750 | 0.0751
$\theta_{5}-\theta_{2}$ | -0.2000 | 0.0000 | -0.2000 | -0.2014
$\theta_{4}-\theta_{3}$ | 0.0250 | 0.0509 | 0.0509 | 0.0509
$\theta_{6}-\theta_{5}$ | -0.1000 | 0.0000 | -0.1000 | -0.1002
$\theta_{7}-\theta_{5}$ | 0.0000 | 0.0753 | 0.0753 | 0.0754
$\theta_{7}-\theta_{6}$ | 0.1000 | 0.1221 | 0.1221 | 0.1222
Figure 1: Bottleneck network with 7 nodes. Red circles represent motor nodes
while green squares represent generator nodes.

##  III Use of symmetry in the transient analysis of the linear swing equation

Knowledge of the symmetries may be used to reduce the computational burden
associated with the modal decomposition. In fact, we can apply the modal
decomposition to a reduced ‘quotient network’ where all the redundancies are
eliminated, thereby making the analysis easier and faster.

We introduce the $n\times q$ indicator matrix $E$, where $n$ is the number of
nodes and $q$ is the number of orbits or clusters. Each entry $E(i,j)=1$ if
node $i$ belongs to the cluster $j$ and 0 otherwise. Another way to write the
indicator matrix is the following,
$E=[\mathbf{e}_{1},\mathbf{e}_{2},...,\mathbf{e}_{q}]$ where each vector
$\mathbf{e_{i}}$ is an indicator vector, that is entry $j$ of vector
$\mathbf{e}_{i}$ is 1 if node $j$ is in cluster $i$ and is 0 otherwise.

######  Lemma 6.

The set of states such that $\theta_{l}(t)=\theta_{m}(t)$ and
$\dot{\theta}_{l}(t)=\dot{\theta}_{m}(t)$ for all $l,m\in S_{i}$, $i=1,...,q$,
define an invariant manifold [36].

The full network with Laplacian matrix $L$ can be transformed into its
corresponding quotient network with Laplacian matrix $\tilde{L}$ by the
transformation [41],

$\tilde{L}=((E^{T}E)^{-1}E^{T})LE.\\\ $

By pre-multiplying Eq. (8) by $((E^{T}E)^{-1}E^{T})$, we obtain,

$\ddot{\tilde{\boldsymbol{\theta}}}(t)=-\gamma\dot{\tilde{\boldsymbol{\theta}}}(t)+\tilde{L}\tilde{\boldsymbol{\theta}}(t)+\tilde{\mathbf{p}},$
(27)

where the $q$-dimensional vectors
${\tilde{\boldsymbol{\theta}}}(t)=((E^{T}E)^{-1}E^{T}){\boldsymbol{\theta}}(t)$
and $\tilde{\mathbf{p}}=[\tilde{p}_{1},\tilde{p}_{2},...,\tilde{p}_{q}]$ is
equal to $\tilde{\mathbf{p}}=((E^{T}E)^{-1}E^{T}){\textbf{p}}$, .

######  Definition 6.

We say that the power vector p respects the symmetries if for all the nodes
$i$ in the same cluster $S_{j}$, $p_{i}=\tilde{p}_{j}$, $j=1,...,q$.

In what follows we first consider the case that the power vector p respects
the symmetries and then consider the more general case that it does not.

###  III-A The case that the power vector $\mathbf{p}$ respects the symmetries

When the power vector respects the symmetries, the forced evolution for the
full network and the quotient network will be _the same_ , as we explain in
what follows.

######  Definition 7.

A redundant eigenvector $\boldsymbol{\rho}$ of the Laplacian matrix $L$ has
sum of entries corresponding to each cluster equal to zero and a non redundant
eigenvector $\boldsymbol{\nu}$ of the Laplacian matrix have entries
corresponding to each cluster that are all the same [46]. All the redundant
eigenvectors are orthogonal to $\mathbf{e}_{i}$, $i=1,..,q$. Contrarily, the
non-redundant eigenvectors are not. That means
$\mathbf{e}_{i}^{T}\boldsymbol{\rho}=0$ and
$\mathbf{e}_{i}^{T}\boldsymbol{\nu}\neq 0$, $i=1,..,q$, where
$\boldsymbol{\rho}$ is any redundant vector, and $\boldsymbol{\nu}$ any non
redundant vector.

From Definition 7 we see that the matrix $L$ has $q$ non-redundant
eigenvectors $\boldsymbol{\nu}$’s and $(n-q)$ redundant eigenvectors
$\boldsymbol{\rho}$’s [46].

######  Remark 2.

From Definition 7 it follows that the modes (10) can be either redundant or
nonredundant. The numer of redundant modes is equal to $(n-q)$ and the numer
of nonredundant modes is equal to $q$.

All of the eigenvalues of the matrix $\tilde{L}$ are also eigenvalues of the
matrix $L$ [46]. However, the matrix $L$ has additional eigenvalues that are
not eigenvalues of $\tilde{L}$. When an eigenvalue of $L$ is also an
eigenvalue of $\tilde{L}$, we say that it is ‘inherited’ by the quotient
network.

######  Lemma 7.

Only the eigenvalues of $L$ associated with non-redundant eigenvectors are
inherited by the quotient network.

###### Proof.

Consider the equation:

$\tilde{L}=(E^{T}E)^{-1}E^{T}LE=(E^{T}E)^{-1}E^{T}(VDV^{T})E,$ (28)

where $V$ is the matrix of the eigenvectors of $L$ and $D$ is the matrix of
eigenvalues of $L$.

The matrix $V=[V_{\rho},V_{\nu}]$ where $V_{\rho}$ has all the redundant
eigenvectors and $V_{\nu}$ has all the non-reduntant eigenvectors. Further,
$D$ can be written as $D_{\rho}\oplus D_{\nu}$, where $D_{\rho}$ is a diagonal
matrix with all the $(n-q)$ redundant eigenvalues and $D_{\nu}$ is a diagonal
matrix with all the $q$ nonredundant eigenvalues. Now we see that
$E^{T}V_{\rho}=0$ due to the orthogonality. Therefore,

$\tilde{L}=(E^{T}E)^{-1}E^{T}(V_{\nu}D_{\nu}V_{\nu}^{T})E=(E^{T}E)^{-1}(HD_{\nu}H^{T}),$
(29)

where $H=E^{T}V_{\nu}$ is a square $q$-dimensional matrix.

∎

From Eq. (29), we can see that only non-redundant eigenvalues and eigenvectors
are inherited by the quotient network. Moreover since the power vector
$\mathbf{p}$ respects the symmetries, the entries of the vector
$\tilde{{p}}_{j}=p_{i}$ for any node $i$ in cluster $S_{j}$ (see Definition
6).

From $\sum_{j}{p}_{j}=0$ it follows that
$\sum_{j=1}^{m}=n_{j}{\tilde{p}}_{j}=0$, where $n_{j}$ is the number of nodes
in each cluster.

######  Remark 3.

In the case that the network is connected, the eigenvector of the Laplacian
matrix $[11...1]$, associated with the zero eigenvalue, is non-redundant. So,
the zero eigenvalue is inherited by the quotient network. Further, it can also
be proved under generic conditions that in the case of a disconnected network,
the quotient network inherits $c$ zero eigenvalues from the full network,
where $c$ is the number of components.

######  Theorem 1.

Assume the network is connected. The entry of the vector $\tilde{\mathbf{q}}$
corresponding to the zero eigenvalue is 0.

###### Proof.

Consider the equation:

$\mathbf{w^{T}}\tilde{\mathbf{p}}=0$ (30)

Let $\mathbf{v}$ be a right eigenvector and $\mathbf{w^{T}}$ be a left
eigenvector of $\tilde{L}$. We know the right eigenvector associated with the
zero eigenvalue is $\mathbf{v}=[11...1]$. Which implies,

$\tilde{L}\mathbf{v}=0$

$(E^{T}E)^{-1}E^{T}LE\mathbf{v}=0$

By left multiplying both sides of the above equation by $(E^{T}E)$, we obtain,

$E^{T}LE\mathbf{v}=0$

Indicating that the matrix $E^{T}LE$ must have $\mathbf{v}=[11...1]$ as its
right eigenvector associated with the zero eigenvalue. A left eigenvector
$\mathbf{w}^{T}$ associated with zero eigenvalue has to satisfy the equation:

$\mathbf{w^{T}}(E^{T}E)^{-1}E^{T}LE=0$ (31)

This also means $\mathbf{w^{T}}(E^{T}E)^{-1}$ is a left eigenvector with
associated eigenvalue 0 for $(E^{T}LE)$ which is a symmetric matrix. Because
$E^{T}LE$ is symmetric, its left eigenvectors and right eigenvectors,
corresponding to the same eigenvalues, are the same.
Ergo,

$\mathbf{w^{T}}(E^{T}E)^{-1}=\mathbf{v^{T}}$ (32)

Recall that $\mathbf{v^{T}}=[11....1]$. Moreover, the matrix $(E^{T}E)$ is a
diagonal matrix that has the number of nodes in each cluster $n_{1},n_{2},$
etc on the main diagonal. It follows that $\mathbf{w^{T}}=[n_{1}n_{2}...]$.
Since the entries of the vector $\mathbf{w^{T}}$ are the number of nodes in
each cluster, the product $\mathbf{w^{T}}\tilde{\mathbf{p}}$ is equal to the
sum of all the powers in the full network which is zero by the assumption that
the full network is balanced. ∎

In order to demonstrate this quotient transformation, we will use the $n=7$
node graph with $6$ edges shown in Fig. 2, for which

Figure 2: 7-node network. Generator nodes are indicated with squares and motor
nodes with circles. All nodes corresponding to the same cluster or orbit of
the symmetry group are colored the same.

$L$=$\begin{bmatrix}-3&1&1&1&0&0&0\\\ 1&-2&0&0&1&0&0\\\ 1&0&-2&0&0&1&0\\\
1&0&0&-2&0&0&1\\\ 0&1&0&0&-1&0&0\\\ 0&0&1&0&0&-1&0\\\
0&0&0&1&0&0&-1\end{bmatrix}$

with eigenvalues $0,-0.38,-0.38,-1.58,-2.61,-2.61,-4.41$. The power vector
$\mathbf{p}$ and the transformed power vector $\mathbf{q}$ are also shown,

$\mathbf{p}$= $\begin{bmatrix}-0.3\\\ -0.3\\\ -0.3\\\ -0.3\\\ 0.4\\\ 0.4\\\ 0.4\end{bmatrix}$ |  $\mathbf{q}$= $\begin{bmatrix}0.00\\\ 0.00\\\ 0.00\\\ -0.8895\\\ 0.00\\\ 0.00\\\ 0.2208\end{bmatrix}$
---|---

As can be seen the power vector $\mathbf{p}$ respects the symmetries.

Figure 3 shows the flows vs. time over all the network lines (forced
evolution). Due to the network symmetries, several flows are superimposed on
top of each other (and thus indistinguishable).

Figure 3: Flows vs. time for network in Fig. 2. The network has a total of 6
lines. The forced evolution of the flows over each one of these lines is
plotted. However, due to the symmetries, triplets of flows are superimposed on
top of each other.

For the network in Fig. 2, $S_{1}=\\{1\\}$, $S_{2}=\\{2,3,4\\}$,
$S_{3}=\\{5,6,7\\}$. This corresponds to the following indicator matrix:

E= $\begin{bmatrix}1&0&0\\\ 0&1&0\\\ 0&1&0\\\ 0&1&0\\\ 0&0&1\\\ 0&0&1\\\ 0&0&1\end{bmatrix}$ |
---|---

The Laplacian matrix $\tilde{L}$ for the quotient network is,

Figure 4: Quotient Network corresponding to the network in Fig. 2.

$\tilde{L}$= $\begin{bmatrix}-3&3&0\\\ 1&-2&1\\\ 0&1&-1\end{bmatrix}$.
---

with eigenvalues $0,-1.58,-4.41$. Moreover,

$\tilde{V}_{r}$= $\begin{bmatrix}0.898&0.731&0.577\\\ -0.423&0.345&0.577\\\
0.124&-0.589&0.577\end{bmatrix}$
---
$\tilde{V}_{l}$= $\begin{bmatrix}0.562&0.337&0.229\\\ -0.794&0.476&0.688\\\
0.233&-0.813&0.688\end{bmatrix}$
---

$\tilde{\mathbf{p}}$= $\begin{bmatrix}-0.3\\\ -0.3\\\ 0.4\end{bmatrix}$ |  $\tilde{\mathbf{q}}$= $\begin{bmatrix}0.1873\\\ -0.6402\\\ 0\end{bmatrix}$,
---|---

where $\tilde{V}_{r}$ and $\tilde{V}_{l}$ represent the matrix of right and
the left eigenvectors of $\tilde{L}$ respectively. The vectors
$\tilde{\mathbf{p}}$ and $\tilde{\mathbf{q}}$ are the power vector and the
transformed power vector of the quotient network respectively.
As seen in the Fig. 4, the quotient network has three nodes (each representing
a cluster) and two lines. Unlike the full network, this network is directed
and weighted. Under the assumption that the vector p respects the symmetries,
a study of the forced response of the quotient network provides complete
information about the forced response of the full network. If we plot the
flows of the quotient network, we still get Fig. 3 but with only one curve
representing all the nodes in a cluster.

###  III-B The case that the power vector $\mathbf{p}$ does not respect the
symmetries

When the power vector does not respect the symmetries, we can still reduce the
dynamics in the quotient network form (27). In this case the power for each
node of the quotient network will be the average of the powers of the nodes in
each cluster and the forced evolution of the quotient network is the average
of the forced evolutions of the full network within each cluster.

As an example, consider the network in Fig. 2 studied previously but with
power vector,

$\mathbf{p}$= $\begin{bmatrix}-0.3\\\ -0.3\\\ -0.4\\\ -0.3\\\ 0.5\\\ 0.4\\\
0.4\end{bmatrix}$

which does not respect the symmetries. The forced evolutions for the full
network and for the quotient network are given in Figs. 5(a) and 5(b). The
color code consistent with color of nodes in Figs. 2 and 4. All the green
(red) curves in Fig. 5(a) average to the one green (red) curve in Fig. 5(b).

Figure 5: a:$\theta$ vs. time for the network in Fig. 2, b: $\theta$ vs. time
for quotient network in Fig. 4. Color code consistent with color of nodes in
Figs. 2 and 4.

An important observation is that the quotient network removes redundant modes
but these modes might still produce overshoots in the full network which will
be ignored by the quotient network.

Since the forced evolution obtained from the quotient network is not the same
as that of the full network, a relevant problem is to characterize deviations
between the time evolutions of the nodes in the same cluster and the time
evolution of the associated quotient node.

######  Definition 8.

Define the deviations as
$\epsilon_{i}(t)=(\theta_{i}(t)-\tilde{\theta}_{{i*}}(t))$, $i=1,2,...,n$,
where ${i*}$ is the cluster node $i$ belongs to, $\sum_{i\in
S_{j}}\epsilon_{i}(t)=0$.

The dynamics of the deviations is provided by the redundant modes only.
Without loss of generality assume the redundant modes (10) are labeled
$i=1,...,(n-q).$ Then

$\epsilon_{i}(t)=\sum_{j=1}^{n-q}V_{ij}\eta_{j}(t).$ (33)

The modal decomposition of the forced response Eq. (33) can be used to compute
the transient deviation overshoot dynamics, peak times, upper bounds, etc. We
can also approximate the peaks using the numerical approach of Sec. IIA.

######  Lemma 8.

In the case the power vector p respects the symmetries, the redundant modes
(10) have $q_{i}=0$, hence their forced response is zero.

###### Proof.

This follows trivially from Definitions 6 and 7. ∎

###  III-C Large Network Example

We now consider a large network, shown in Figure 7, based off the Chilean
power-grid topology with $n=218$ nodes, out of which $94$ are motors and $124$
are generators. This is the reduced version achieved by the Star-Mesh
transformation [6] of the original Chilean power-grid [25].

####  III-C1 Powers that respect the symmetries

We first choose the power vector such that it respects the symmetries; for
simplicity we set all generator nodes to have a power of 0.0081 and all motor
nodes to have a power of -0.0106. For this network, we have 29 non-trivial
clusters; clusters with more than one node, and 96 trivial clusters; clusters
with only one node. Using this information, we can transform the full network
into the quotient network shown in Figure 8. After the transformation, each
cluster becomes a node in the quotient network with 157 nodes, out of which 74
are generator nodes and 83 are motor nodes. As already discussed previously,
if the vector p respects the symmetries, the forced response of the quotient
network coincides the forced response of the full network (not shown here for
the sake of brevity.)

####  III-C2 Powers that do not respect the symmetries

If the power vector does not respect the symmetries, the quotient network time
evolution provides the average forced evolution of the nodes in each cluster
of the full network. The deviation between the dynamics of the full network
and the quotient network can be characterized using Eq. (33). To demonstrate
this we first choose a cluster where the nodes do not respect the symmetries.
The cluster is chosen containing 7 nodes: 178, 181, 182, 182, 183, 184, 185,
207, with associated powers [0.0181, 0.0081, 0.0081, 0.0081, 0.0081, 0.0081,
0.0081]. Now, using analogous equations to Eq. (25) and (26), we can obtain
the maximum displacement deviation of each node $\epsilon_{i}$. Since we are
finding the maximum of displacements rather than flows, we replace the
coefficients $C_{i}$ with $V_{i}$ and get,

${\sum_{k=1}^{n-q}}\frac{V_{k}q_{k}}{\varsigma_{k}}\sin{\varsigma_{k}t}=0.$
(34)
$\tilde{\tau}_{ij}=\sqrt{\frac{6\sum_{k=1}^{n-q}V_{k}q_{k}}{\sum_{k=1}^{n-q}V_{k}q_{k}\varsigma_{k}^{2}}}.$
(35)

Using the technique mentioned in Section IIA and Eq. (34) and (35), we can
easily calculate maximum displacement deviation. The time evolution of this
deviation for our selected cluster is shown in Fig. 6.

Figure 6: Evolution of the displacement deviation of the quotient network node
with respect to the original nodal displacements. The blue curve is for node
178, the red curve is for nodes 181, 182, 183, 184, 185, 207. Figure 7:
Chilean power-grid Network. Squares indicate generator nodes and circles
indicate motor nodes. Black nodes are in trivial clusters whereas nodes with
the same color (different from black) are in non-trivial clusters. Figure 8:
Quotient network for Chilean Power-grid topology. Squares indicate generator
nodes and circles indicate motor nodes. Black nodes are in trivial clusters
whereas nodes with the same color (different from black) are in non-trivial
clusters. Thickness of the arrows are indicative of the weight of the edges.

##  IV Effects of heterogeneity in the $\gamma_{i}$ terms

In this section we remove the assumption of homogeneous damping terms
$\gamma_{i}=\gamma$. Note that the symmetries of the swing equation are
defined in the generic case of heterogeneous (arbitrary) $\gamma$’s, see
Definition 4. Here we explain how our previous derivations need to be modified
in the case of heterogeneous $\gamma$’s.

Equation (8) is rewritten,

$\ddot{\boldsymbol{\theta}}(t)=-\Gamma\dot{\boldsymbol{\theta}}(t)+L\boldsymbol{\theta}(t)+\mathbf{p},$
(36)

where the $n$-dimensional diagonal matrix
$\Gamma=diag({\gamma}_{1},{\gamma}_{2},...,{\gamma}_{n})$. Analogously, the
quotient network Equation (27) becomes,

$\ddot{\tilde{\boldsymbol{\theta}}}(t)=-\tilde{\Gamma}\dot{\tilde{\boldsymbol{\theta}}}(t)+\tilde{L}\tilde{\boldsymbol{\theta}}(t)+\tilde{\mathbf{p}},$
(37)

where the $q$-dimensional diagonal matrix
$\tilde{\Gamma}=diag(\tilde{\gamma}_{1},\tilde{\gamma}_{2},...,\tilde{\gamma}_{q})$
and $\tilde{\gamma}_{i}$ is the damping ratio of all the nodes in cluster
$S_{i}$.

A separate consideration is needed for the study of the deviations in the case
that (i) the $\gamma$’s are heterogeneous and (ii) the power vector p does not
respect the symmetries. To deal with such a case, we introduce the irreducible
representations of the symmetry group $\mathcal{G}$ [36].

From knowledge of the group of symmetries $\mathcal{G}$, we can compute the
irreducible representations (IRRs) of the symmetry group of the network. This
defines a transformation $T$ into the so-called IRR coordinate system (see
Ref. [36]). The transformation matrix $T$ is orthogonal. Each one of the rows
of the matrix $T$ is associated with a specific cluster. If a row of the
matrix $T$ is associated with cluster $k$, it means all the $i$ entries of
that row are zero for $i$ not in cluster $S_{k}$. The first $q$ rows of the
matrix $T$ are parallel to the $q$ nonredundant eigenvectors
$\boldsymbol{\nu}$.

Therefore, the matrix $T$ can be cast in the following block diagonal form,

$T=\bigoplus_{k=1,..,q}T_{k},$ (38)

where each block $T_{k}$ is an $n_{k}$-dimensional square matrix associated
with cluster $k$. Premultiplying Eq. (36) by $T$, we obtain,

$\ddot{\boldsymbol{\xi}}(t)=-\Gamma\dot{\boldsymbol{\xi}}(t)+B\boldsymbol{\xi}(t)+\mathbf{r},$
(39)

where $\boldsymbol{\xi}=T\boldsymbol{\theta}$, $\mathbf{r}=T\mathbf{p}$, the
block -diagonal matrix $B=TLT^{T}$ and the matrix $\Gamma$ is unaltered by
application of the matrix $T$. The latter property follows from the block-
diagonal structure of the matrix $T$ and the observation that $\Gamma$ can be
recast in a similar form

$\Gamma=\bigoplus_{k=1,...,q}I_{n_{k}}\tilde{\gamma}_{k},$ (40)

where $I_{n_{k}}$ is the $n_{k}$-dimensional identity matrix.

The transformed $n\times n$ block diagonal matrix $B=TAT^{T}$ is a direct sum
$\oplus_{u=1}^{U}\hat{B}_{u}$, where $\hat{B}_{u}$ is a (generally complex)
$p_{u}\times p_{u}$ matrix with $p_{u}$ the multiplicity of the $u$th IRR in
the permutation representation, $U$ the number of IRRs present and $d_{u}$ the
dimension of the $u$th IRR, so that $\sum_{u}d_{u}p_{u}=n$. The matrix $T$
contains information on which perturbations affecting different clusters get
mapped to different IRRs [42].

Due to the block-diagonal structure of the matrix $B$, Eq. (39) can be
decoupled into a number of lower dimensional equations, where each equation
corresponds to a block of the matrix $B$.

There is one representation (labeled $u=1$) which we call trivial and has
dimension $d_{1}=q$. The associated block of the matrix $B$ corresponds to the
dynamics of the quotient network. Hence, the trivial representation is
associated with all the clusters. However, it is possible that other IRR
representations are only associated with some of the clusters (not all of
them.) Each one of these other representations $u>1$ describes the deviation
dynamics of either an isolated cluster or a group of intertwined clusters
[36]. A simple interpretation of isolated vs. intertwined clusters is the
following. If a cluster is isolated a perturbation affecting the power of any
one of its nodes will not affect the deviation dynamics of other clusters. On
the contrary, when a set of two or more clusters are intertwined, a
perturbation affecting the power of any of the nodes in a cluster will affect
the deviation dynamics of the remaining clusters in the set.

### Example A

For this example, we consider the dynamics of Eq. (5), with adjacency matrix
corresponding to the network in Fig. 2, the power vector
$\mathbf{p}=[-0.3,-0.3,-0.3,-0.4,0.4,0.4,0.5]^{T}$ and the matrix
$\Gamma=diag(3,2,1,1,2,1,1)$. By the definition of symmetries (Definition 3),
nodes belonging to the same cluster have the same $\gamma$ value. In this
case, due to the presence of heterogeneity in $\gamma$, there are two non-
trivial clusters, each containing 2 nodes: $\\{3,4\\}$ and $\\{6,7\\}$. All
other nodes are in trivial clusters. Moreover, the power vector $\mathbf{p}$
does not respect the symmetries.

We compute the $IRRs$ of the symmetry group. The transformation matrix $T$ is
equal to,

$T$= $\begin{bmatrix}1&0&0&0&0&0&0\\\ 0&1&0&0&0&0&0\\\
0&0&\frac{1}{\sqrt{2}}&\frac{1}{\sqrt{2}}&0&0&0\\\ 0&0&0&0&1&0&0\\\
0&0&0&0&0&\frac{1}{\sqrt{2}}&\frac{1}{\sqrt{2}}\\\
0&0&\frac{1}{\sqrt{2}}&-\frac{1}{\sqrt{2}}&0&0&0\\\
0&0&0&0&0&\frac{1}{\sqrt{2}}&-\frac{1}{\sqrt{2}}&\end{bmatrix}$
---

With this information, we can calculate the block-diagonal matrix $B=TLT^{T}$,
${\textbf{r}}=T{\textbf{p}}$, and $\Gamma_{orth}=T\Gamma T^{T}$,

$B=\left[\begin{array}[]{@{}ccccc|cc@{}}-3&1&\sqrt{2}&0&0&0&0\\\
1&-2&0&1&0&0&0\\\ \sqrt{2}&0&-2&0&1&0&0\\\ 0&1&0&-1&0&0&0\\\ 0&0&1&0&-1&0&0\\\
\hline\cr 0&0&0&0&0&-{2}&{1}\\\ 0&0&0&0&0&{1}&-{1}\end{array}.\right]$

$\mathbf{r}$=
$\begin{bmatrix}-0.3&-0.3&-0.495&0.4&0.636&\vline{1/10\sqrt{2}}&{-1/10\sqrt{2}}\end{bmatrix}^{T}$.

$\Gamma_{orth}=\left[\begin{array}[]{@{}ccccc|cc@{}}3&0&0&0&0&0&0\\\
0&2&0&0&0&0&0\\\ 0&0&1&0&0&0&0\\\ 0&0&0&2&0&0&0\\\ 0&0&0&0&1&0&0\\\ \hline\cr
0&0&0&0&0&{1}&0\\\ 0&0&0&0&0&0&{1}\\\ \end{array}.\right]$

The transverse (bottom right) block of the matrix $B$ represents the deviation
dynamics for the two non-trivial clusters, which are intertwined. Note that
the elements of $\mathbf{r}$ relative to this block are given by
$r_{6}=\frac{(p_{3}-r_{4})}{\sqrt{2}}$ and
$r_{7}=\frac{(p_{6}-r_{7})}{\sqrt{2}}$. We also note that $\Gamma_{orth}$ for
the transverse block is the same as $\Gamma$ as explained in the theory. After
diagonalizing the transverse block, its dynamics decouples into two
independent modes,

$\ddot{\boldsymbol{\kappa}}=-I\dot{\boldsymbol{\kappa}}+\begin{bmatrix}-0.3820&0\\\
0&-2.6180\end{bmatrix}{\boldsymbol{\kappa}}+\begin{bmatrix}0.02298\\\
-0.09732\end{bmatrix}$ (41)

Now, using the technique derived in Section IIA, we can use Eq. (41) to
calculate the maximum deviations from the quotient network. We can find the
initial guess analogously to Eq. (35). $\tilde{\tau}$ can then be used to find
the time of maximum error by solving the following equation for $t$ until
convergence:

${\sum_{k=m+1}^{n}}\frac{V_{k}r_{k}}{\varsigma_{k}}\sin{\varsigma_{k}t}=0.$
(42)

The obtained time $t$ can then be used to solve (19) which gives us the
deviations in each non-trivial cluster as seen in Table II. Table II also
shows that for this example our approach based on the linearized swing
equation (5) well approximates the maximum error obtained by integration of
the full nonlinear swing equation (3).

TABLE II:  Steady state, first peak and maximum errors for the two non-trivial clusters. The values for the second column are obtained using Eq. (20) and the values for the third column are obtained using Eq. (35), Eq. (42) and Eq. (19). The fourth and fifth columns are calculated by numerically solving the linear swing equation (5) and the non-linear swing equation (3), respectively. |  |  | Linear | Non-Linear
---|---|---|---|---
Deviation | Steady | First | Max | Max
| State | Peak | Error | Error
$\frac{\theta_{4}-\theta_{3}}{\sqrt{2}}$ | 0 | -0.0314 | -0.0314 | -0.0323
$\frac{\theta_{7}-\theta_{6}}{\sqrt{2}}$ | 0.0707 | 0.0707 | 0.0713 | 0.0792

### Example B

For this example, we consider the Chilean power-grid described in Section IIIC
and pictured in Fig. 7. The power vector is chosen so that it does not respect
the symmetries. The coefficients $\gamma_{i}$ are chosen such that all nodes
have $\gamma_{i}=1$, except for $\gamma_{116}=2$. For this network we computed
the $IRR$ representations and the block-diagonal matrix $B$. This allowed us
to study the deviation from the quotient network dynamics for each one of the
$29$ nontrivial clusters. We did not find any intertwined clusters. Next, we
focus on a cluster comprised of the three nodes $\\{77,79,81\\}$, represented
as three cyan squares on the right side of Fig. 7. The transverse block
corresponding to that cluster is reduced in a similar form to Eq. (41),

$\ddot{\boldsymbol{\kappa}}=-I\dot{\boldsymbol{\kappa}}+\begin{bmatrix}-9.123&0\\\
0.0000&-0.876\end{bmatrix}{\boldsymbol{\kappa}}+\begin{bmatrix}-0.00502\\\
0.00643\end{bmatrix}$ (43)

Similarly to Example A, Eq. (43) can be mapped to the deviation of the
individual nodal displacement from the quotient displacement. This can be done
by computing a linear combination of the uncoupled modes (43) and applying the
technique developed in Sec. IIA. Results of this computation are presented in
Table III, which shows our ability to accurately predict the maximum
deviations for all three nodes $\\{77,79,81\\}$ in the cluster.

TABLE III:  Steady state, first peak and maximum errors for the cluster containing nodes $\\{77,79,81\\}$. The values for the second column are obtained using Eq. (20) and the values for the third column are obtained using Eq. (35), Eq. (42) and Eq. (19). The fourth and fifth columns are calculated by numerically integrating the linear and non linear ODE respectively. $\theta_{q}$ refers to the displacement of the quotient node which is given by $\theta_{q}=\frac{\theta_{77}+\theta_{79}+\theta_{81}}{3}$ |  |  | Linear | Non-Linear
---|---|---|---|---
Deviation | Steady | First | Max | Max
| State | Peak | Error | Error
$\theta_{77}-\theta_{q}$ | 0.00083 | 0.001181 | 0.01175 | 0.001182
$\theta_{79}-\theta_{q}$ | -0.000416 | -0.000591 | -0.000588 | -0.000577
$\theta_{81}-\theta_{q}$ | -0.000416 | -0.000591 | -0.000588 | -0.000577

##  V Conclusions

In this paper we have studied how the presence of network symmetries affects
the swing equation dynamics. We have first introduced the nonlinear swing
equation and then modeled the propagation of small perturbations via the
linearized swing equation. We have then shown that the nonlinear and the
linear swing equation have the same symmetries. These symmetries can be
reduced to provide an essential description of the dynamics in terms of a
’quotient network’.
We have introduced a decomposition of the linear swing equation dynamics into
independent modes, each of which corresponds to a forced second order system.
This allowed us to characterize several transient effects such as overshoots,
peaks, peak times, upper bounds etc. We have classified the symmetries into
two cases based on whether the power vector respects the symmetries or not and
presented techniques to obtain a full characterization of the transient
dynamics and/or error quantification.
In the case in which the power demanded and supplied by different nodes
respects the symmetries, the quotient network completely describes the forced
evolution of the full network, otherwise it provides the forced evolution
averaged over the nodes in each cluster. We have also introduced the error
dynamics which describes how the quotient network time evolutions deviate from
the time evolutions of the individual nodes inside the clusters. This error
dynamics can be written as a linear combination of the ‘redundant modes’.
Finally, we have presented how the symmetry analysis can be applied to
networks with heterogeneous $\gamma$. In order to study the deviation dynamics
in the case that the $\gamma$’s are arbitrary and the power vector does not
respect the symmetries, we have introduced the irreducible representation of
the symmetry group. Furthermore, though not discussed in the paper, it is
possible to apply our techniques to the case of different types of forcing
terms, such as sinusoidal forcing; they can also be generalized to other
network models, such as the effective network model (EN) or the structure
preserving model (SP) [31].

## Acknowledgement

This work was supported by the National Science Foundation through grant
(Grant No. 1727948.) The authors thank Cesar Ornelas for collaborating on an
earlier version of this paper, Aranya Chakraborrti for insightful discussions
and Matteo Lodi for helping with the calculations of the irreducible
representations.

|  Kshitij Bhatta received a bachelor’s degree in Mechanical Engineering in
2020 from the University of New Mexico (US). He is currently a Master’s
student in Mechanical Engineering at the University of New Mexico. His
research interests include non-linear dynamic systems, analysis of complex
networks and control systems with feedback in robotic and automotive
applications.
---|---
|  Majeed M. Hayat received the M.S. and the Ph.D. degrees in Electrical and
Computer Engineering from the University of Wisconsin-Madison in 1988 and
1992, respectively. He is currently a Professor and Department Chair of
Electrical and Computer Engineering at Marquette University, Milwaukee, WI,
USA. He was Associate Editor of Optics Express (Photodetectors and Image
Processing) from 2004 to 2010 and Associate Editor and member of the
Conference Editorial Board for the IEEE Control Systems Society. From 2010 to
2013 he was the Chair of the topical committee on Photodetectors, Sensors,
Systems, and Imaging of the IEEE Photonics Society. From 2014 to 2018 he
served as Associate Editor for the IEEE Transactions on Parallel and
Distributed Systems. His research activities cover a broad range of topics
including resilience and reliability of interdependent cyber-physical systems,
dynamical modeling of cascading phenomena with applications to power systems,
avalanche photodiodes, statistical communication theory, signal and image
processing, algorithms for spectral and radar sensing and imaging, optical
communication, and networked computing. He is a recipient of the National
Science Foundation Early Faculty Career Award (1998). Dr. Hayat has authored
or co-authored over 108 peer-reviewed journal articles and 134 conference
papers (over 5,500 citations, H-Index: 37), and has fourteen issued patents,
six of which have been licensed.
---|---
|  Francesco Sorrentino received a master’s degree in Industrial Engineering
in 2003 and a Ph.D. in Control Engineering in 2007 both from the University of
Naples Federico II (Italy). His expertise is in dynamical systems and
controls, with particular emphasis on nonlinear dynamics and adaptive
decentralized control. His work includes studies on dynamics and control of
complex dynamical networks and hypernetworks, adaptation in complex systems,
sensor adaptive networks, coordinated autonomous vehicles operating in a
dynamically changing environment, and identification of nonlinear systems. He
is interested in applying the theory of dynamical systems to model, analyze,
and control the dynamics of complex distributed energy systems such as power
networks and smart grids. Subjects of current investigation are evolutionary
game theory on networks (evolutionary graph theory), the dynamics of large
networks of coupled neurons, and the use of adaptive techniques for dynamical
identification of communication delays between coupled mobile platforms. He
has published more than 50 papers in international scientific peer reviewed
journals.
---|---

## References

  * [1] Réka Albert, István Albert, and Gary L Nakarado.  Structural vulnerability of the north american power grid.  Physical review E, 69(2):025103, 2004.
  * [2] Sergio Arianos, Ettore Bompard, Anna Carbone, and Fei Xue.  Power grid vulnerability: A complex network approach.  Chaos: An Interdisciplinary Journal of Nonlinear Science, 19(1):013119, 2009.
  * [3] Ross Baldick, Badrul Chowdhury, Ian Dobson, DOng Zhaoyang, Bei Gou, and David Hawkins.  Vulnerability assessment for cascading failures in electric power systems.  PES Power Systems Conference and Exposition, 2009.
  * [4] Bassam Bamieh and Dennice F Gayme.  The price of synchrony: Resistive losses due to phase synchronization in power networks.  In 2013 American Control Conference, pages 5815–5820. IEEE, 2013\.
  * [5] Wayne Barrett, Amanda Francis, and Benjamin Webb.  Equitable decompositions of graphs with symmetries.  Linear Algebra and its Applications, 513:409–434, 2017.
  * [6] S Bedrosian.  Converse of the star-mesh transformation.  IRE Transactions on Circuit Theory, 1961.
  * [7] Andrey Bernstein, Daniel Bienstock, David Hay, Meric Uzunoglu, and Gil Zussman.  Sensitivity analysis of the power grid vulnerability to large-scale cascading failures.  ACM SIGMETRICS Performance Evaluation Review, pages 33–37, 2012\.
  * [8] Karen A Blaha, Ke Huang, Fabio Della Rossa, Louis Pecora, Mani Hossein-Zadeh, and Francesco Sorrentino.  Cluster synchronization in multilayer networks: A fully analog experiment with l c oscillators with physically dissimilar coupling.  Physical Review Letters, 122(1):014101, 2019.
  * [9] Benjamin A Carreras, Vickie E Lynch, Ian Dobson, and David E Newman.  Critical points and transitions in an electric power transmission model for cascading failure blackouts.  Chaos: An interdisciplinary journal of nonlinear science, 12(4):985–994, 2002.
  * [10] David P Chassin and Christian Posse.  Evaluating north american electric grid reliability using the barabási–albert network model.  Physica A: Statistical Mechanics and its Applications, 355(2-4):667–677, 2005.
  * [11] Quiming Chen, Chuanwen Jian, Wenzheng Qiu, and J.D. McCalley.  Probability models for estimating the probabilities of cascading outages in high-voltage transmission network.  IEEE Transactions on Power Systems, pages 1423–1431, 2006.
  * [12] Young Sul Cho, Takashi Nishikawa, and Adilson E Motter.  Stable chimeras and independently synchronizable clusters.  Physical review letters, 119(8):084101, 2017.
  * [13] Paolo Crucitti, Vito Latora, and Massimo Marchiori.  Model for cascading failures in complex networks.  Physical Review E, 69(4):045104, 2004.
  * [14] Paolo Crucitti, Vito Latora, and Massimo Marchiori.  A topological analysis of the italian electric power grid.  Physica A: Statistical mechanics and its applications, 338(1-2):92–97, 2004.
  * [15] Fabio Della Rossa, Louis Pecora, Karen Blaha, Afroza Shirin, Isaac Klickstein, and Francesco Sorrentino.  Symmetries and cluster synchronization in multilayer networks.  Nature communications, 11(1):1–17, 2020.
  * [16] Ajendra Dwivedi, Xinghuo Yu, and Peter Sokolowski.  Identifying vulnerable lines in a power network using complex network theory.  In Industrial Electronics, 2009. ISIE 2009. IEEE International Symposium on, pages 18–23. IEEE, 2009.
  * [17] Ajendra Dwivedi, Xinghuo Yu, and Peter Sokolowski.  Analyzing power network vulnerability with maximum flow based centrality approach.  In Industrial Informatics (INDIN), 2010 8th IEEE International Conference on, pages 336–341. IEEE, 2010.
  * [18] Martin Golubitsky and Ian Stewart.  Symmetry and pattern formation in coupled cell networks.  In Pattern formation in continuous and coupled systems, pages 65–82. Springer, 1999.
  * [19] Linqi Guo, Changhong Zhao, and Steven H Low.  Graph laplacian spectrum and primary frequency regulation.  In 2018 IEEE Conference on Decision and Control (CDC), pages 158–165. IEEE, 2018.
  * [20] Zhang Guohua, Wang Ce, Zhang Jianhua, Yang Jingyan, Zhang Yin, and Duan Manyin.  Vulnerability assessment of bulk power grid based on complex network theory.  In Electric Utility Deregulation and Restructuring and Power Technologies, 2008. DRPT 2008. Third International Conference on, pages 1554–1558. IEEE, 2008.
  * [21] Åke J Holmgren.  Using graph models to analyze the vulnerability of electric power networks.  Risk analysis, 26(4):955–969, 2006.
  * [22] Yubo Huang, Junguo Lu, and Weidong Zhang.  Small disturbances can trigger cascading failures in power grids, 2019\.
  * [23] Takayuki Ishizaki, Aranya Chakraborty, and Jun-ichi Imura.  Graph-theoretic analysis of power systems.  Proceedings of the IEEE, 106, 2018.
  * [24] Stefan Kettemann.  Delocalization of disturbances and the stability of ac electricity grids.  Physical Review E, 94(6):062311, 2016.
  * [25] Heetae Kim, David Olave-Rojas, Eduardo Alvarez-Miranda, and Seung Woo Son.  In-depth data on the network structure and hourly activity of the central chilean power grid.  Scientific Data, 2018.
  * [26] Petar Mlinarić, Takayuki Ishizaki, Aranya Chakrabortty, Sara Grundel, Peter Benner, and Jun-ichi Imura.  Synchronization and aggregation of nonlinear power systems with consideration of bus network structures.  2018 European Control Conference (ECC), 2018.
  * [27] Adilson E Motter and Ying-Cheng Lai.  Cascade-based attacks on complex networks.  Physical Review E, 66(6):065102, 2002.
  * [28] Adilson E Motter, Seth A Myers, Marian Anghel, and Takashi Nishikawa.  Spontaneous synchrony in power-grid networks.  Nature Physics.
  * [29] Tommaso Nesti, Alessandro Zocca, and Bert Zwart.  Emergent failures and cascades in power grids: a statistical physics perspective.  Physical review letters, 2018.
  * [30] Vincenzo Nicosia, Miguel Valencia, Mario Chavez, Albert Díaz-Guilera, and Vito Latora.  Remote synchronization reveals network symmetries and functional modules.  Physical review letters, 110(17):174102, 2013.
  * [31] Takashi Nishikawa and Adilson E Motter.  Comparative analysis of existing models for power-grid synchronization.  New Journal of Physics, 17(1), 2015.
  * [32] Koji Okuda and Yoshiki Kuramoto.  Mutual entrainment between populations of coupled oscillators.  Progress of theoretical physics, 86(6):1159–1176, 1991.
  * [33] Sakshi Pahwa, Amelia Hodges, Caterina Scoglio, and Sean Wood.  Topological analysis of the power grid and mitigation strategies against cascading failures.  In 2010 IEEE International Systems Conference, pages 272–276. IEEE, 2010.
  * [34] Sakshi Pahwa, Caterina Scoglio, and Antonio Scala.  Abruptness of cascade failures in power grids.  Scientific reports, 4:3694, 2014.
  * [35] Fabio Pasqualetti, Antonio Bicchi, and Francesco Bullo.  Graph-theoretical characterization of power network vulnerabilities.  Proceedings of the American Control Conference, 2011.
  * [36] Louis M Pecora, Francesco Sorrentino, Aaron M Hagerstrom, Thomas E Murphy, and Rajarshi Roy.  Cluster synchronization and isolated desynchronization in complex networks with symmetries.  Nature Communications, 5, 2014.
  * [37] Graham Rogers.  Power System Oscillations.  Springer Science Business Media.
  * [38] Martin Rohden, Daniel Jung, Samyak Tamrakar, and Stefan Kettemann.  Cascading failures in ac electricity grids.  Physical Review E, 94(3):032209, 2016.
  * [39] Marti Rosas-Casals, Sergi Valverde, and Ricard V Solé.  Topological vulnerability of the european power grid under errors and attacks.  International Journal of Bifurcation and Chaos, 17(07):2465–2475, 2007.
  * [40] Benjamin Schäfer, Dirk Witthaut, Marc Timme, and Vito Latora.  Dynamically induced cascading failures in power grids.  Nature communications, 9(1):1–13, 2018.
  * [41] Michael T Schaub, Neave O’Clery, Yazan N Billeh, Jean-Charles Delvenne, Renaud Lambiotte, and Mauricio Barahona.  Graph partitions and cluster synchronization in networks of oscillators.  Chaos: An Interdisciplinary Journal of Nonlinear Science, 26(9):094821, 2016.
  * [42] Abu Bakar Siddique, Louis Pecora, Joseph D Hart, and Francesco Sorrentino.  Symmetry-and input-cluster synchronization in networks.  Physical Review E, 97(4):042217, 2018.
  * [43] Ricard V Solé, Martí Rosas-Casals, Bernat Corominas-Murtra, and Sergi Valverde.  Robustness of the european power grids under intentional attack.  Physical Review E, 77(2):026102, 2008.
  * [44] Francesco Sorrentino and Louis M Pecora.  Approximate cluster synchronization in networks with symmetries.  Chaos, 26(094823).
  * [45] Francesco Sorrentino, Abu Bakar Siddique, and Louis M Pecora.  Symmetries in the time-averaged dynamics of networks: Reducing unnecessary complexity through minimal network models.  Chaos: An Interdisciplinary Journal of Nonlinear Science, 29(1):011101, 2019.
  * [46] Rubén J. Sánchez-García.  Exploiting symmetry in network analysis.  Communications Physics, 3(1), May 2020.
  * [47] Melvyn Tyloo and Philippe Jacquod.  Primary control effort under fluctuating power generation in realistic high-voltage power networks.  IEEE Control Systems Letters, 5(3):929–934, 2020.
  * [48] Melvyn Tyloo, L. Pagnier, and Philippe Jacquod.  The key player problem in complex oscillator networks and electric power grids: Resistance centralities identify local vulnerabilities.  Science Advances, 2019.
  * [49] Thierry Van Cutsem and Costas Vournas.  Voltage stability of electric power systems.  Springer Science & Business Media, 2007.
  * [50] Thanh Long Vu and Konstantin Turitsyn.  Lyapunov functions family approach to transient stability assessment.  IEEE Transactions on Power Systems, 31(2):1269–1277, 2015.
  * [51] Andrew J Whalen, Sean N Brennan, Timothy D Sauer, and Steven J Schiff.  Observability and controllability of nonlinear networks: The role of symmetry.  Physical Review X, 5(1):011005, 2015.
  * [52] Yang Yang, Takashi Nishikawa, and Adilson Motter.  Small vulnerable sets determine large network cascades in power grids.  Science, 2017.
  * [53] Honghao Zheng and Christopher DeMarco.  A bi-stable branch model for energy-based cascading failure analysis in power systems.  North American Power Symposium 2010, 2010.