Applied Koopman Operator Theory for Power Systems Technology

TL;DR

This paper reviews Koopman operator theory for nonlinear systems and demonstrates its application to power systems for stability analysis, diagnostics, and data-driven insights without relying on explicit models.

Contribution

It introduces a data-centric approach using Koopman spectral analysis for power system stability and diagnostics, expanding the application of this theory beyond traditional nonlinear analysis.

Findings

Koopman spectral analysis enables stability assessment without explicit models.

Application of Koopman theory improves detection of system instabilities.

Data-driven methods facilitate power system diagnostics and coherency identification.

Abstract

Koopman operator is a composition operator defined for a dynamical system described by nonlinear differential or difference equation. Although the original system is nonlinear and evolves on a finite-dimensional state space, the Koopman operator itself is linear but infinite-dimensional (evolves on a function space). This linear operator captures the full information of the dynamics described by the original nonlinear system. In particular, spectral properties of the Koopman operator play a crucial role in analyzing the original system. In the first part of this paper, we review the so-called Koopman operator theory for nonlinear dynamical systems, with emphasis on modal decomposition and computation that are direct to wide applications. Then, in the second part, we present a series of applications of the Koopman operator theory to power systems technology. The applications are established as data-centric methods, namely, how to use massive quantities of data obtained numerically and experimentally, through spectral analysis of the Koopman operator: coherency identification of swings in coupled synchronous generators, precursor diagnostic of instabilities in the coupled swing dynamics, and stability assessment of power systems without any use of mathematical models. Future problems of this research direction are identified in the last concluding part of this paper.

Tables

Table 1. Table 1: Numerical results on Koopman eigenvalues and modes obtained with Arnoldi-type algorithm. Reproduced from [ 42 ]

KM pair	Growth Rate	Argument [rad]	Frequency [Hz]	Norm
$j$	$\left|{\tilde{\lambda }}_j\right|$			$\Vert {\widetilde{𝑽}}_j\Vert$
1	0.9986	$\pm$0.1701	1.3533	3.0021
2	0.9986	$\pm$0.1438	1.1447	2.3930
3	0.9985	$\pm$0.1009	0.8028	0.7039
4	0.9985	$\pm$0.1300	1.0343	0.9753
5	0.9984	$\pm$0.0931	0.7405	0.4507
6	0.9984	$\pm$0.1130	0.8990	0.8162
7	0.9983	$\pm$0.1643	1.3078	6.6147
8	0.9983	$\pm$0.1378	1.0962	7.1941
9	0.9983	$\pm$0.0468	0.3727	2.1006
10	0.9982	$\pm$0.1836	1.4612	1.2238

Equations

$$\frac{d}{dt}\mbox{\boldmath$x$}(t)=\mbox{\boldmath$F$}(\mbox{\boldmath$x$}(t)),$$

$$\mbox{\boldmath$x$}[k+1]=\mbox{\boldmath$T$}(\mbox{\boldmath$x$}[k]),$$

$$\mbox{\boldmath$S$}_{t\to t+h}:{\mathbb{X}}\to{\mathbb{X}};~{}~{}\mbox{\boldmath$x$}(t)\mapsto\mbox{\boldmath$x$}(t+h)=\mbox{\boldmath$x$}(t)+\int^{t+h}_{t}\mbox{\boldmath$F$}(\mbox{\boldmath$x$}(\tau))d\tau.$$

$$f:\mathbb{X}\to \mathbb{C}.$$

$${\mathbf{U}}_{t}f$$

$$\mathbf{U}f$$

$$y(t):=f(\mbox{\boldmath$x$}(t))=f(\mbox{\boldmath$S$}_{t}(\mbox{\boldmath$x$}(0)))=f\circ\mbox{\boldmath$S$}_{t}(\mbox{\boldmath$x$}(0))=({\mathbf{U}}_{t}f)(\mbox{\boldmath$x$}(0)).$$

$${\mathbf{U}}_{t}(\alpha_{1}f_{1}+\alpha_{2}f_{2})=(\alpha_{1}f_{1}+\alpha_{2}f_{2})\circ\mbox{\boldmath$S$}_{t}=\alpha_{1}(f_{1}\circ\mbox{\boldmath$S$}_{t})+\alpha_{2}(f_{2}\circ\mbox{\boldmath$S$}_{t})=\alpha_{1}({\mathbf{U}}_{t}f_{1})+\alpha_{2}({\mathbf{U}}_{t}f_{2}).$$

$${\mathbf{U}}_t{\varphi }_j={\mathrm{e}}^{{\nu }_{j}t}{\varphi }_j,\text{makebox}[20.00003pt]j=1,2,\dots$$

$$\mathbf{U}{\varphi }_j={\lambda }_j{\varphi }_j,\text{makebox}[20.00003pt]j=1,2,\dots$$

$$\frac{d}{dt}\theta (t)=\Omega ,\text{makebox}[20.00003pt]\Omega :\text{real constant}.$$

$$({\mathbf{U}}_t{\varphi }_j)(\theta )={\varphi }_j(S_t(\theta ))={\mathrm{e}}^{\mathrm{i}j(\theta +\Omega t)}={\mathrm{e}}^{\mathrm{i}j\Omega t}{\mathrm{e}}^{\mathrm{i}j\theta }={\mathrm{e}}^{\mathrm{i}j\Omega t}{\varphi }_j(\theta ).$$

$$({\mathbf{U}}_{t}f)(\theta )=\sum _{j\in \mathbb{Z}}{c}_j({\mathbf{U}}_t{\varphi }_j)(\theta )=\sum _{j\in \mathbb{Z}}{c}_j{\varphi }_j(\theta ){\mathrm{e}}^{\mathrm{i}j\Omega t},$$

$$\frac{d}{dt}\theta (t)=\omega (t),\text{makebox}[20.00003pt]\frac{d}{dt}\omega (t)=0.$$

$$f(\theta ,\omega )=\omega +\sum _{j\in \mathbb{Z}\setminus \{0\}}{a}_j(\omega ){\mathrm{e}}^{\mathrm{i}j\theta },$$

$$({\mathbf{U}}_{t}f)(\theta ,\omega )=\omega +\sum _{j\in \mathbb{Z}\setminus \{0\}}{a}_j(\omega ){\mathrm{e}}^{\mathrm{i}j(\theta +\omega t)}=\omega +\sum _{j\in \mathbb{Z}\setminus \{0\}}\left\{a_j(\omega ){\mathrm{e}}^{\mathrm{i}j\theta }\right\}{\mathrm{e}}^{\mathrm{i}j\omega t}.$$

$$\mathbf{U}={\mathbf{U}}_{\mathrm{s}}+{\mathbf{U}}_{\mathrm{r}},$$

$${\mathbf{U}}_{\mathrm{s}}=\sum _{j=1}^{\infty }{\lambda }_j{\mathbf{P}}_j,$$

$${\mathbf{U}}_{\mathrm{r}}={\int }_{\mathbb{T}}{\mathrm{e}}^{\mathrm{i}\theta }dE(\theta ).$$

$$y[k]:=y(kT)=({\mathbf{U}}_{kT}f)(\mbox{\boldmath$x$}(0))=(({\mathbf{U}}_{T})^{k}f)(\mbox{\boldmath$x$}(0)).$$

$$f(\mbox{\boldmath$x$}[k])=({\mathbf{U}}^{k}f)(\mbox{\boldmath$x$}[0])=({\mathbf{U}}^{k}_{\rm s}f+{\mathbf{U}}^{k}_{\rm r}f)(\mbox{\boldmath$x$}[0])=\sum_{j=1}^{\infty}\lambda^{k}_{j}({\mathbf{P}}_{j}f)(\mbox{\boldmath$x$}[0])+\int_{\mathbb{T}}{\rm e}^{{\rm i}k\theta}dE(\theta)f(\mbox{\boldmath$x$}[0]),$$

$$f(\mbox{\boldmath$x$}[k])=\sum_{j=1}^{\infty}\lambda^{k}_{j}\phi_{j}(\mbox{\boldmath$x$}[0])V_{j}+\int_{\mathbb{T}}{\rm e}^{{\rm i}k\theta}dE(\theta)f(\mbox{\boldmath$x$}[0]).$$

$$\mbox{\boldmath$f$}(\mbox{\boldmath$x$}[k])=\sum_{j=1}^{\infty}\lambda^{k}_{j}\phi_{j}(\mbox{\boldmath$x$}[0])\mbox{\boldmath$V$}_{j}+\left[\begin{array}[]{c}{\displaystyle\int_{\mathbb{T}}{\rm e}^{{\rm i}k\theta}dE(\theta)f_{1}(\mbox{\boldmath$x$}[0])}\\ \vdots\\ {\displaystyle\int_{\mathbb{T}}{\rm e}^{{\rm i}k\theta}dE(\theta)f_{m}(\mbox{\boldmath$x$}[0])}\end{array}\right].$$

$${\mathbf{P}}^{\nu}f(\mbox{\boldmath$x$}[0])=\lim_{n\rightarrow\infty}\frac{1}{n}\sum^{n-1}_{k=0}{\rm e}^{-{\rm i}2\pi k\nu}f(\mbox{\boldmath$x$}[k]),\makebox[20.00003pt]{}f\in{\mathcal{F}},$$

$$\left[\begin{array}[]{c}{\mathbf{P}}^{\nu_{j}}f_{1}(\mbox{\boldmath$x$}[0])\\ \vdots\\ {\mathbf{P}}^{\nu_{j}}f_{m}(\mbox{\boldmath$x$}[0])\end{array}\right]=\phi_{j}(\mbox{\boldmath$x$}[0])\mbox{\boldmath$V$}_{j},$$

$$\mbox{\boldmath$f$}=\sum_{j=1}^{\infty}\phi_{j}\mbox{\boldmath$V$}_{j},$$

$$\mbox{\boldmath$f$}(\mbox{\boldmath$x$}[k])=\sum^{\infty}_{j=1}\phi_{j}(\mbox{\boldmath$x$}[k])\mbox{\boldmath$V$}_{j}=\sum^{\infty}_{j=1}{{\mathbf{U}}}^{k}\phi_{j}(\mbox{\boldmath$x$}[0])\mbox{\boldmath$V$}_{j}=\sum^{\infty}_{j=1}\lambda^{k}_{j}\phi_{j}(\mbox{\boldmath$x$}[0])\mbox{\boldmath$V$}_{j}.$$

$$\mathcal{K}_{N}:=\mathrm{span}\{\mathbf{\underline{U}}^{k}\mbox{\boldmath$f$};k=0,\ldots,N-1\},\makebox[20.00003pt]{}\mathbf{\underline{U}}^{k}\mbox{\boldmath$f$}:=[{\mathbf{U}}^{k}f_{1},\ldots,{\mathbf{U}}^{k}f_{m}]^{\top}.$$

$${\mathbf{P}}_N\underline{\mathbf{U}}:{\mathcal{K}}_N\to {\mathcal{K}}_N$$

$$\{\mbox{\boldmath$y$}_{0},\mbox{\boldmath$y$}_{1},\ldots,\mbox{\boldmath$y$}_{N-1}\},$$

$$\mbox{\boldmath$r$}=\mbox{\boldmath$y$}_{N-1}-\mathsf{K}_{N-1}\mbox{\boldmath$c$},$$

Full text

1 Introduction

Complex phenomena of the electric power system occur on a wide range of scales in both space and time. The power system has a more than 100 years of history in operation and is well-known as a large-scale, man-made, networked system of heterogeneous components with internal dynamics. The dynamic components include generation units (plants), loads, transmission lines, and—as a current particular interest—distributed renewable resources such as solar, wind, geothermal, and cogeneration plants. Examples of the complex phenomena are synchronization of individual rotating machines, voltage dynamics and collapse, and cascading outages leading to wide-area blackouts. Cascading outages have been recorded in practical power systems [1, 2, 3, 4, 5, 6, 7, 8]. The details of the phenomena have been studied with direct numerical simulations of nonlinear mathematical models: see, e.g., [3, 9, 10, 11, 12, 13]. Analyzing such complex phenomena is of basic importance for realizing the so-called “smart” society in the future.

In this paper, we have a modest but counterintuitive goal: to develop methodology and tools for analysis of such complex nonlinear power systems through linear techniques. Needless to say, it has been demonstrated that a low-dimensional nonlinear dynamical system exhibits very complicated motions such as chaotic oscillation [14]. For power system applications, such low-dimensional models play a crucial role in revealing the key mechanisms of dynamic phenomena and instabilities: see, e.g., [15, 16, 17, 18]. This direction of analyzing the nonlinear models is based on the state description, called the Poincaré picture of dynamics, and has great success in nonlinear dynamical systems [19, 20] and nonlinear control systems [21]. However, it has been recognized that a high-dimensional nonlinear model with a huge number of states can not be easily handled and understood. To resolve the difficulty, as a counter-part of this description, B. O. Koopman is regarded as the pioneer for the use of linear transformations on Hilbert space to analyze Hamiltonian systems by introducing the so-called Koopman operator and studying its spectrum [22]: see [23, 24] in detail. Even if the governing system is nonlinear and finite-dimensional, the Koopman operator is linear but infinite-dimensional and does not rely on linearization: indeed, it captures the full information of the nonlinear dynamical system. In this setting, the elegant set of linear algebra, function space, and operator theory can be used to analyze the nonlinear system while avoiding the handling of a huge number of states. In addition to the linear treatment, the Koopman picture of dynamics has a close connection to data-centric requirement to analysis of complex dynamics [25, 26], where multiple domains—dynamical systems theory, large-scale simulation, sensing, data mining, and so on—can work together. This is currently desirable in power systems technology according to the new development of real-time Phasor Measurement Units (PMUs), which offers an advanced data collection method using phasors of AC (Alternative-Current) voltages: see, e.g., [27, 28, 29].

The purposes of this paper are two-fold. One is to provide a brief review of the Koopman operator theory of nonlinear dynamical systems. There is a large set of literature of theoretical progress in the Koopman operator: see, e.g., [30, 31, 32, 33, 34, 35, 36, 37, 38, 39] and [25, 26] for comprehensive reviews. In particular, it is shown in [31] via spectral analysis of the Koopman operator that single-frequency modes can be embedded in complex spatiotemporal dynamics. These modes are later named the Koopman Modes (KMs) [40] and used for characterizing global behaviors of complex fluid flows: see [40, 41] and references of [26]. Due to high-dimensional, spatiotemporal nature of the complex phenomena in the power system, it is of basic interest for practitioners to identify a small number of dominant components or modes that approximates the phenomena observed practically and numerically. In the first part of this paper, we review a recent progress of the Koopman operator theory with emphasis on modal decomposition and computation, which we believe are direct to wide applications beyond power systems technology addressed in this paper. The first part of this paper is generally described and not specific to any application domain.

The other purpose is to present a series of applications of the Koopman operator theory to power systems technology that we have recently published since 2010. The applications are established as data-centric methods, namely, how to use massive quantities of data obtained numerically and experimentally: coherency identification of swings in coupled synchronous generators [42], precursor diagnostic of instabilities in the coupled swing dynamics [43], and stability assessment without any use of models [44]. It should be noted that other applications in power and energy systems have been reported: for example, network partitioning related to avoidance of cascading failures [45, 46], analysis of PMU data in Japan [47], analysis of building dynamics [48, 49], quantification of smoothing effects in wind energy [50]. Outside the fluid mechanics and power systems technology, the Koopman operator theory has been applied to background/foreground separation in video [51], data fusion [52], dynamic texture [53], and so on. We will try to describe these applications so that readers can easily understand how the data-centric power system analysis is developed via spectral properties of the Koopman operator.

The rest of this paper is organized as follows: In Section 2 we provide a summarized theory of the Koopman operator for nonlinear dynamical systems based on the literature mentioned above. In Sections 3 to 5 we present a series of applications of the Koopman operator theory in power systems technology based on our previous papers [42, 43, 44]. Conclusion of this paper is made in Section 6 with future problems.

Notation   All vectors are viewed as columns. The sets of all real numbers, complex numbers, and integers are denoted by ${\mathbb{R}}$, ${\mathbb{C}}$, and ${\mathbb{Z}}$, respectively. We use $\mathbb{R}_{\geq 0}:=\{x:x\in{\mathbb{R}},x\geq 0\}$. The symbol $\mathbb{T}$ represents the one-torus or circle. The symbol ${\rm i}$ stands for the imaginary number, $|z|$ for the modulus of $z\in{\mathbb{C}}$, $\mathrm{Arg}(z)$ for its argument, and $z^{\ast}$ for the complex conjugate. For a vector \mbox{\boldmathz}\in{\mathbb{C}}^{m} ($m$ times), \mbox{\boldmathz}^{\top} stands for its transpose, \mbox{\boldmathz}^{\ast} for the complex-conjugate transpose, and ||\mbox{\boldmathz}||:=\sqrt{\mbox{\boldmathz}^{\ast}\mbox{\boldmathz}} for the ${\mathbb{C}}^{m}$-norm. We use \mbox{\boldmathz}^{\mathrm{c}}:=(\mbox{\boldmathz}^{\ast})^{\top}. The symbols of transpose and complex-conjugate transpose are also used for matrices. For a matrix $\mathsf{A}$, its rank is denoted as $\mathrm{rank}(\mathsf{A})$, and the linear space spanned by all linear combinations of the column vectors in $\mathsf{A}$ as $\mathrm{span}\{\mathsf{A}\}$. We use $\mathrm{diag}(a_{1},a_{2},\ldots)$ for a diagonal matrix whose diagonal entries starting in the upper left corner are $a_{1},a_{2},...$ The symbol $\bot$ indicates the orthogonality between vector and space.

2 Summarized theory of Koopman operators

This section provides a summarized theory of Koopman operators for nonlinear dynamical systems. Comprehensive reviews of the emerging theory exist in [25, 26]. Current emphasis is on computation as a tool for how the theory is utilized for solving concrete dynamical problems, which we believe is direct to wide applications beyond power systems.

2.1 Definition

Throughout the paper, we consider a dynamical system evolving on a finite $n$-dimensional manifold ${\mathbb{X}}$: for continuous time $t\in{\mathbb{R}}$,

[TABLE]

and for discrete time $k\in{\mathbb{Z}}$,

[TABLE]

where $x$ is called the state belonging to the state space ${\mathbb{X}}$. The function \mbox{\boldmathF}:{\mathbb{X}}\to\mathrm{T}\mathbb{X} (tangent bundle of $\mathbb{X}$) is a nonlinear vector-valued function and is assumed to be tractable in the region we are interested, and \mbox{\boldmathT}:{\mathbb{X}}\to{\mathbb{X}} a nonlinear vector-valued map assumed in the same manner. For the continuous-time system (1), the following finite time-$h$ map \mbox{\boldmathS}_{t\to t+h} ($h\in{\mathbb{R}}$) is defined as

[TABLE]

Since (1) is autonomous, the map does not depend on $t$, simply denoted as \mbox{\boldmathS}_{h}. The one-parameter group of maps, \{\mbox{\boldmathS}_{h};h\in{\mathbb{R}}\}, is called the flow. Below, we will use $t$ for the suffix of the flow instead of $h$ for clarity of the presentation.

Here we introduce the Koopman operators for the continuous- and discrete-time dynamical systems. To do this, the so-called observable $f$ is introduced as a scalar-valued function defined on the state space ${\mathbb{X}}$, namely

[TABLE]

This observable is a mathematical formulation of observation of the state’s dynamics and is well-known as the output function in the control engineering community. Below, we will denote by ${\mathcal{F}}$ a given space of observables (its properties will be mentioned later). For an observable $f\in{\mathcal{F}}$, the Koopman operator ${\mathbf{U}}_{t}$ for the continuous-time system (1) and ${\mathbf{U}}$ for the discrete-time system (2) map $f$ into a new function as follows:

[TABLE]

That is, the Koopman operator ${\mathbf{U}}_{t}$ (or ${\mathbf{U}}$) describes the time $t$ (or $k$) evolution of observation (output of the system) along the state’s dynamics: for the continuous-time case, the output $y(t)$ is described as

[TABLE]

Although the dynamical systems (1) and (2) can be nonlinear and evolve in the finite-dimensional space, the Koopman operators are linear but infinite-dimensional. The linearity is easily proven: for the continuous-time case, by picking up any two observables $f_{1},f_{2}\in{\mathcal{F}}$ and two scalars $\alpha_{1},\alpha_{2}\in{\mathbb{C}}$, we have

[TABLE]

This type of composition operator is defined for a large class of nonlinear dynamical systems [24] and does not rely on linearization: indeed, it captures the full information on the original nonlinear systems. The essential idea outlined in this paper is clear: to analyze dynamics described by the nonlinear systems (1) and (2) through the linear operators.

Here it is noted which class of function space is chosen for ${\mathcal{F}}$ of observables. In the traditional ergodic theory, which has mainly studied measure-preserving dynamical systems on a compact space, the standard $L^{2}$ space is chosen. For this space, the celebrated Birkhoff’s point-wise ergodic theorem and so on are presented: see [54] in detail. For a general, possibly dissipative dynamical system on a non-compact space, there is still room for consideration. Even in the dissipative system, if it has a compact attractor, the $L^{2}$ space is effectively used for investigating asymptotic dynamics on the attractor through the Koopman operator. In [34], for a nonlinear system with a stable equilibrium point, the space of analytic functions in the domain of attraction of the point is chosen as ${\mathcal{F}}$. Also, in [37] generalized Hardy space is introduced as ${\mathcal{F}}$ for analysis of off-attractor dynamics of a nonlinear system with a stable equilibrium point or limit cycle. As a summary, the choice of a space of analytic functions is better for analyzing dissipative dynamics via the Koopman operator with a relatively-simple spectrum.444except for a nonlinear system with a chaotic attractor

2.2 Spectral properties

Spectrum plays an important role in understanding the Koopman operators and, in our purpose, investigating the original nonlinear systems. Since the linear operators are infinite-dimensional, it is necessary to consider both discrete and continuous spectra [55, 56].

First of all, let us introduce the discrete part of spectrum, that is, eigenvalues of the Koopman operators. In this part, the concept in finite-dimensional matrices is naturally applied. For the continuous-time system (1), eigenfunctions and eigenvalues of ${\mathbf{U}}_{t}$ are defined as follows: for non-zero $\phi_{j}\in{\mathcal{F}}$ and $\nu_{j}\in{\mathbb{C}}$ such that

[TABLE]

The functions $\phi_{j}$ are referred to as Koopman eigenfuctions, and the constants $\nu_{j}$ as the associated Koopman Eigenvalues (KEs). The set of all the KEs is called the discrete spectrum of ${\mathbf{U}}_{t}$. Also, for the discrete-time system (2), the definitions of Koopman eigenfunctions and eigenvalues are almost same: for non-zero $\phi_{j}\in{\mathcal{F}}$ and $\lambda_{j}\in{\mathbb{C}}$ such that

[TABLE]

For example, consider the following first-order continuous-time system: for $\theta\in\mathbb{T}$,

[TABLE]

The corresponding flow is $S_{t}(\theta)=\theta+{\mathit{\Omega}}t~{}(\textrm{mod}\,2\pi)$. It is easily checked that the $j$-th Fourier basis $\phi_{j}(\theta)={\rm e}^{{\rm i}j\theta}$ satisfies (4) with KE $\nu_{j}={\rm i}j{\mathit{\Omega}}$:

[TABLE]

For a suitable choice of ${\mathcal{F}}$, we can expand the observable $f(\theta)$ in terms of the Fourier bases and obtain its composition with $S_{t}$:

[TABLE]

where $c_{j}\in{\mathbb{C}}$ is a coefficient for the expansion. The formula (6) suggests that the time evolution of observable under (5) is expressed as an infinite summation over $\mathbb{Z}$ in oscillations with distinct, constant (angular) frequencies $j{\mathit{\Omega}}$ of KEs.

Next, we mention the continuous part of spectrum that is of particular interest in infinite-dimensional linear operators. The mathematical definition of continuous spectrum is available in textbooks (see, e.g., [55, 56]), and we now describe its dynamical motivation with the simple, second-order continuous-time linear system: for $(\theta,\omega)\in\mathbb{T}\times{\mathbb{R}}$,

[TABLE]

The induced flow is \mbox{\boldmathS}_{t}(\theta,\omega)=[\theta+\omega t\,(\textrm{mod}\,2\pi),\omega]^{\top}. Obviously, the value [math] corresponds to an eigenvalue with eigenfunction $\phi_{0}(\theta,\omega)=\omega$, which can be replaced by any function of $\omega$. Here, since the state space is periodic in $\theta$, in the similar manner as above, we may expand the observable $f(\theta,\omega)$ in terms of the Fourier bases like

[TABLE]

where $a_{j}(\omega)$ is a function dependent on $f$. Thus, we have the following composition of $f$ with \mbox{\boldmathS}_{t}:

[TABLE]

The formula is apparently similar to (6), but the angular frequencies depend on the state $\omega$ that can be taken in ${\mathbb{R}}$. This suggests that no decomposition of the time evolution $f(\theta(t),\omega(t))$ based on distinct, constant frequencies is obtained. Rather than this, the time evolution is represented in terms of continuous quantity—integral over ${\mathbb{R}}$. This is a situation in which the continuous spectrum appears in the Koopman operator.

Now, based on [31], we introduce the spectral representation of the Koopman operator ${\mathbf{U}}$ that is unitary on ${\mathcal{F}}=L^{2}({\mathbb{X}})$ with a compact ${\mathbb{X}}$. The unitary case has been traditionally studied in ergodic theory [23] because it results from a measure-preserving map on a compact state space. The following representation holds for the same as for the continuous-time case. For unitary ${\mathbf{U}}$, its spectrum is restricted to the unit circle in the complex domain. Thus, ${\mathbf{U}}$ admits a unique decomposition into its singular and regular parts:

[TABLE]

where there exists a (Hilbert) subspace $\mathcal{H}\subset{\mathcal{F}}$ such that for any $f\in{\mathcal{F}}$, the new function ${\mathbf{U}}_{\rm s}f$ belongs in $\mathcal{H}$, and ${\mathbf{U}}_{\rm r}f$ does in the orthogonal complement of $\mathcal{H}$ . The singular part ${\mathbf{U}}_{\rm s}$ implies the discrete spectrum, determined by the KEs of ${\mathbf{U}}$:

[TABLE]

where ${\mathbf{P}}_{j}$ is the projection operator associated with KE $\lambda_{j}$ of unity modulus and eigenfunction $\phi_{j}$; ${\mathbf{P}}_{i}\circ{\mathbf{P}}_{j}=\delta_{ij}{\mathbf{P}}_{i}$ where $\delta_{ij}$ is the Kronecker delta. The regular part ${\mathbf{U}}_{\rm r}$ implies the continuous spectrum and is represented with a continuous spectral measure $E(\theta)$ on $\mathbb{T}$ as

[TABLE]

Last of all, in terms of data-centric applications, we introduce a connection of discrete spectra between the continuous- and discrete-time systems. For this, consider a uniformly-sampled observation by letting $T(>0)$ be the sampling period. From (3) we have

[TABLE]

This clearly shows that the sampled observation of the dynamics described by the continuous-time system (1) is modeled by the Koopman operator ${\mathbf{U}}={\mathbf{U}}_{T}$ for the discrete-time system (2) with eigenfunctions $\phi_{j}$ and eigenvalues $\lambda_{j}={\rm e}^{\nu_{j}T}$. Note that the Koopman eigenfunctions $\phi_{j}$ are common for the original continuous- and derived discrete-time systems. Since the main idea in this paper is to analyze the dynamics from sampled observational data, we mainly focus on the discrete-time setting in the following.

2.3 Koopman modes: formulation and computation

2.3.1 Mathematical formulation

The idea in [31] is to analyze dynamics governed by the nonlinear system (2) using spectrum of the Koopman operator. First of all, let us consider the time evolution of observable, f(\mbox{\boldmathx}[k]), for the unitary case of ${\mathbf{U}}$. Now, from its spectral decomposition in (7), we have

[TABLE]

where we used the property of the projection operator ${\mathbf{P}}_{j}$. By using the formula ${\mathbf{P}}_{j}f=\phi_{j}V_{j}$, where $V_{j}\in{\mathbb{C}}$ is a constant for the projection depending on $f$, we obtain the following spectral representation of f(\mbox{\boldmathx}[k]) in terms of ${\mathbf{U}}$:

[TABLE]

The first term represents the contribution of KEs to \{\mbox{\boldmathf}(\mbox{\boldmathx}[k])\} and describes the average and quasi-periodic parts in \{\mbox{\boldmathf}(\mbox{\boldmathx}[k])\}. The second term represents the contribution of continuous spectrum and describes the aperiodic part contained in \{\mbox{\boldmathf}(\mbox{\boldmathx}[k])\}.

Furthermore, we can extend this representation for a vector-valued observable \mbox{\boldmathf}:=[f_{1},\ldots,f_{m}]^{\top}:{\mathbb{X}}\to{\mathbb{C}}^{m}, where $f_{i}\in{\mathcal{F}}$. The observable corresponds to a vector of any quantities of interest such as frequencies and voltages metered at various points in a power system, which is direct to the applications in Sections 3 to 5. Then, by the similar argument to the scalar-valued observable, the time evolution \mbox{\boldmathf}(\mbox{\boldmathx}[k])=[f_{1}(\mbox{\boldmathx}[k]),\ldots,f_{m}(\mbox{\boldmathx}[k])]^{\top} is exactly represented as follows:

[TABLE]

Thus, we refer to \mbox{\boldmathV}_{j}\in{\mathbb{C}}^{m} as the Koopman Mode (KM) (named in [40]) associated with KE $\lambda_{j}$, corresponding to the observable $f$. The KE $\lambda_{j}$ characterizes the temporal behavior of the corresponding KM \mbox{\boldmathV}_{j}: the argument of $\lambda_{j}$ determines its frequency. The representation (8) shows via spectral analysis of the Koopman operator that single-frequency modes can be embedded in the complex spatio-temporal dynamics.

It is shown in [30, 31] that the terms \phi_{j}(\mbox{\boldmathx}[0])\mbox{\boldmathV}_{j} are defined and computed with a projection operation associated with ${\mathbf{U}}$ applied to the observable $f$. The projection operation of a non-zero ${\mathbf{P}}^{\nu}$ with KE $\lambda={\rm e}^{{\rm i}2\pi\nu}$ is constructively represented as follows:

[TABLE]

where $\nu\in[-1/2,1/2)$. The projections of the $p$ components $f_{1},\ldots,f_{m}$ of $f$ based on ${\mathbf{P}}_{j}$ are thus obtained:

[TABLE]

where $\nu_{j}=\mathrm{Im}[\ln\lambda_{j}]/2\pi$. This formula (9) associates \phi_{j}(\mbox{\boldmathx}[0])\mbox{\boldmathV}_{j} with the projection operation based on ${\mathbf{P}}^{\nu}$. The left-hand sides of (9) are like the discrete Fourier transform of data \{\mbox{\boldmathf}(\mbox{\boldmathx}[k])\}, and hence \phi_{j}(\mbox{\boldmathx}[0])\mbox{\boldmathV}_{j} can be directly computed from data.

Until now, we have assumed that the Koopman operator is unitary, that is, a case where the dynamics of (2) are on an attractor. Even if this is not the case, that is, we consider dynamics off attractors of (2), each KM oscillates with a single frequency. If each of the $p$ components of $f$ lies within the span of eigenfunctions $\phi_{j}$, then, as in [40], we may expand $f$ in terms of the Koopman eigenfunctions as

[TABLE]

where \mbox{\boldmathV}_{j} are regarded as vector coefficients for the expansion. The time evolution \{\mbox{\boldmathf}(\mbox{\boldmathx}[k])\} is identically given as

[TABLE]

Also, if the dynamics observed here have only a finite number of peaks in the frequency domain (practical experience suggests that this situation is normal in power systems), then we can expect the expansion gives an approximation of the dynamics. For dynamics off attractors, the KE $\lambda_{j}$ still characterizes the temporal behavior of the corresponding KM \mbox{\boldmathV}_{j}: the magnitude of $\lambda_{j}$ determines the growth rate. We call the decomposition (10) based on the discrete spectrum of the Koopman operator the Koopman Mode Decomposition (KMD). We will show its capability in power systems in the latter part of this paper, Sections 3 to 5.

2.3.2 Algorithms for computation

For data-centric applications, KMD is considered directly from time-series data without any use of description of the original system, such as explicit description of $T$ and $F$. A common idea for the existing KMD algorithms is to find a finite-dimensional approximation of the action of the Koopman operator directly from finite-time data. The following description in this paragraph is based on [25]. Now consider a fixed vector-valued observable \mbox{\boldmathf}=[f_{i},\ldots,f_{m}]^{\top}\in{\mathcal{F}}^{m} ($m$ times) and associated finite-time data (snapshots) \{\mbox{\boldmathy}_{k}=\mbox{\boldmathf}(\mbox{\boldmathx}[k]);k=0,\ldots,N-1,N<\infty\}. Then, the following function space $\mathcal{K}_{N}$ is sufficient for investigation of the finite-time data:

[TABLE]

In the following, we will assume that \{\mathbf{\underline{U}}^{k}\mbox{\boldmathf};k=0,\ldots,N-1\} is a linearly independent set so that these functions form a basis for $\mathcal{K}_{N}$. Here, let ${\mathbf{P}}_{N}:{\mathcal{F}}^{m}\to\mathcal{K}_{N}$ be a projection from the space of vector-valued observables onto $\mathcal{K}_{N}$. Thus, the following operator

[TABLE]

is finite-dimensional and has a matrix representation, $\mathsf{C}_{N}:{\mathbb{C}}^{N}\to{\mathbb{C}}^{N}$, in terms of the \{\mathbf{\underline{U}}^{k}\mbox{\boldmathf}\}-basis. Note that the matrix $\mathsf{C}_{N}$ is dependent on (i) vector-valued observable $f$, (ii) the number $N$ of snapshots, and (iii) the projection ${\mathbf{P}}_{N}$. If (\lambda,\mbox{\boldmathV}_{\lambda}) is a pair of eigenvalue and eigenvector of $\mathsf{C}_{N}$, where \mbox{\boldmathV}_{\lambda}=[V_{\lambda,0},\ldots,V_{\lambda,N-1}]^{\top}\in{\mathbb{C}}^{N}, then $\lambda$ is an eigenvalue of ${\mathbf{P}}_{N}\mathbf{\underline{U}}$ and \mbox{\boldmath\phi}_{\lambda}:=\sum^{N-1}_{j=0}V_{\lambda,j}(\mathbf{\underline{U}}^{j}\mbox{\boldmathf}) the associated eigenfunction. By restricting our attention to a fixed observable $f$ and associated space $\mathcal{K}_{N}$, we have reduced the problem of locating eigenvalues and modes of the Koopman operator to locating eigenvalues and eigenvectors for a matrix $\mathsf{C}_{N}$. It is shown in [39] that as the number $N$ of snapshots goes to infinity, the KMs located with ${\mathbf{P}}_{N}\mathbf{\underline{U}}$ approach to exact KMs with respect to the fixed $f$ for the Koopman operator $\mathbf{\underline{U}}$. The existing KMD algorithms provide a procedure of how to compute $\mathsf{C}_{N}$ directly from data. Below, we introduce three existing algorithms: Arnoldi-type algorithm [40], dynamic mode decomposition [57, 58, 59, 60], and vector Prony analysis [61].

(a) Arnoldi-type algorithm

First, the Arnoldi-type algorithm [40] for computation of KMs and KEs is introduced. Consider the finite-time data under uniform sampling, again given by

[TABLE]

where \mbox{\boldmathy}_{k}\in{\mathbb{R}}^{m} is the $k$-th snapshot and $N$ the number of available snapshots. Then, the so-called empirical Ritz values $\tilde{\lambda}_{j}$ and empirical Ritz vectors \mbox{\boldmath\tilde{V}}_{j} of the data are defined with the following algorithm:

• (i)

Define a constant vector \mbox{\boldmathc}=[c_{0},\ldots,c_{N-2}]^{\top} such that for vector \mbox{\boldmathr}\in{\mathbb{R}}^{m} satisfying \mbox{\boldmathr}\bot $\mathrm{span}\{\mathsf{K}_{N-1}\}$,

[TABLE]

with the Krylov subspace defined as

[TABLE]

• (ii)

Define the matrix $\mathsf{C}_{N-1}$ as

[TABLE]

called the companion matrix. Then, locate its $N-1$ eigenvalues as the empirical Ritz values $\tilde{\lambda}_{1},\ldots,\tilde{\lambda}_{N-1}$.

• (iii)

Define the Vandermonde matrix $\mathsf{T}$ using $\tilde{\lambda}_{j}$ as

[TABLE]

• (iv)

Define the empirical Ritz vectors \mbox{\boldmath\tilde{V}}_{j} to be the columns of $\mathsf{V}:=\mathsf{K}_{N-1}\mathsf{T}^{-1}$.

It is shown in [40] that if all of the empirical Ritz values are non-zero and distinct, then the following decompositions of the data are obtained:

[TABLE]

Comparing with (10), the empirical Ritz values $\tilde{\lambda}_{j}$ and vectors \mbox{\boldmath\tilde{V}}_{j} behave precisely in the same manner as the KEs $\lambda_{i}$ and the terms \phi_{i}(\mbox{\boldmathx}[0])\mbox{\boldmathV}_{i} containing the Koopman eigenfunctions and KMs, but for the finite sum (15) instead of the infinite sum (10).

Now, it is necessary to show how to choose the constant vector $c$ from the data (11). It is clearly stated in [25] that this is to minimize the ${\mathbb{C}}^{m}$-norm of $r$ in (12), corresponding to choosing the projection operator $\mathbf{P}_{N-1}$ so that \mathbf{P}_{N-1}\mathbf{\underline{U}}^{N-1}\mbox{\boldmathf} is the least-squares approximation to \mathbf{\underline{U}}^{N-1}\mbox{\boldmathf} at \mbox{\boldmathx}_{0}\in\mathbb{X} satisfying \mbox{\boldmathy}_{0}=\mbox{\boldmathf}(\mbox{\boldmathx}_{0}). Here, because of \mbox{\boldmathr}_{N-1}\bot\,\mathrm{span}\{\mathsf{K}_{N-1}\}, by multiplying $\mathsf{K}_{N-1}^{\top}$ on both sides of (12) from left, we have the following linear equation:

[TABLE]

If $\mathsf{K}_{N-1}^{\top}\mathsf{K}_{N-1}\in{\mathbb{R}}^{(N-1)\times(N-1)}$ is regular, then there exists the unique solution $c$ for (16). If this is not the case, we use the Moore-Penrose pseudo-inverse. Especially, for data with low spatial dimension, that is, $m\ll N$ (this situation normally appears in power system applications), because of $\mathrm{rank}(\mathsf{K}_{N-1})\leq m$, the rank-deficiency of $\mathsf{K}_{N-1}$, namely the singularity of $\mathsf{K}_{N-1}^{\top}\mathsf{K}_{N-1}$ is inevitable when choosing $c$. This implies no uniqueness of $c$, and that the derived decomposition and information (frequencies, modes, etc.) should be cross-checked with other means of analysis such as the Fourier-based formula (9) and standard Fourier analysis. We will use the algorithm in Sections 3 to 5 and show the cross-check in Section 3.

Here, based on [58, 25], let us see that the companion matrix $\mathsf{C}_{N-1}$ is a finite-dimensional approximation of the action of $\mathbf{\underline{U}}$. Using $\mathsf{C}_{N-1}$, although not mathematically rigorous, we have

[TABLE]

where \mbox{\boldmathe}^{\top}:=[0,\ldots,0,1]. This suggests that if the last term \mbox{\boldmathr}\mbox{\boldmathe}^{\top} is less dominant (e.g., |\mbox{\boldmathr}| is close to zero), then $\mathsf{C}_{N-1}$ is thought of as an approximation to the action of the Koopman operator onto the finite-dimensional space $\mathcal{K}_{N-1}$. Based on the fact that the Vandermonde matrix diagonalizes the companion matrix as long as its eigenvalues are distinct, namely $\mathsf{C}_{N-1}=\mathsf{T}^{-1}\mathsf{\tilde{\Lambda}}\mathsf{T}$, where $\mathsf{\tilde{\Lambda}}:=\mathrm{diag}(\lambda_{1},\ldots,\lambda_{N-1})$, we have

[TABLE]

That is, if ||\mbox{\boldmathr}\mbox{\boldmathe}^{\top}\mathsf{T}^{-1}|| is small, then $\mathsf{V}$ and $\mathsf{\tilde{\Lambda}}$ approximate eigenvectors and eigenvalues of $\mathbf{\underline{U}}$, respectively.

(b) Dynamic mode decomposition

Second, the Dynamic Mode Decomposition (DMD) is a well-known algorithm in the fluid mechanics community [57, 59]. Connections with the aforementioned Arnoldi-type algorithm are studied in [58, 60]. DMD is also applied in [62] to characterizing global behaviors of power systems. Also, an interesting extension of DMD for numerically approximating not only KEs and KMs but also Koopman eigenfunctions is proposed in [63].

Since many comprehensive literatures on DMD exist, we here provide a brief introduction of DMD. The author of [57] considers the form (13) of a snapshot sequence. Here, let us consider a singular value decomposition of $\mathsf{K}_{N-1}$ as $\mathsf{Q}_{1}\mathsf{\Sigma}\mathsf{Q}_{2}^{\ast}$. From (2.3.2), we roughly have

[TABLE]

This right-hand side can be computed directly from the data and is hence thought of as an approximation of the action of the Koopman operator in the projected space with $\mathsf{Q}_{1}$, that is, $m$-dimensional space. For data with high spatial dimension in fluid applications, namely $m\gg N$, DMD has scored a great success for characterizing complex fluid fields [57, 58, 60]. For data with low spatial dimension ($m\ll N$), the approximation might be not enough to gaining the dynamical information (frequencies, modes, etc.) of \mbox{\boldmathy}_{k}.

(c) Vector Prony analysis

Finally, the so-called vector Prony analysis [61] is introduced and connected to the Arnoldi-type algorithm. Prony analysis is widely used in applications and is a standard technique to reconstruct a scalar function described by a sparse sum of exponentials from finite function values [64, 65, 66, 67]. Since the family of exponentials is utilized, as clearly stated in [65], the Prony analysis is intended to investigating dynamic data as a solution of linear differential equation. The authors of [68] generalize the Prony analysis for reconstruction of sparse sums of eigenfuntions of linear operators such as the differential operator.

Here, let us introduce the procedure of vector Prony analysis. The current problem is to find the best approximation of the data (11) with the $M$ complex values:

[TABLE]

where we call $\hat{\lambda}_{j}\in{\mathbb{C}}$ the $j$-th Prony value and \hat{\mbox{\boldmathV}}_{j}\in{\mathbb{C}}^{m} the associated $j$-th Prony vector. The classical Prony analysis provides a rigorous way to obtain $M=N$ pairs of the two Prony quantities from $2N$ sampled scalars, namely $m=1$. The authors of [61] extended it and provided the following algorithm for obtaining $M=N$ pairs of the Prony values $\hat{\lambda}_{j}$ and vectors \hat{\mbox{\boldmathV}}_{j} from $2N$ sampled vectors \{\mbox{\boldmathy}_{0},\ldots,\mbox{\boldmathy}_{2N-1}\}:

1. (i)

Determine the $N$ coefficients $p_{k}$ of the Prony polynomial defined as

[TABLE]

Clearly, the $N$ roots of the Prony polynomial coincide with the $N$ Prony values. Thus, for $\ell=0,\ldots,N-1$, we have

[TABLE]

This induces

[TABLE]

Therefore, to determine the coefficients $p_{k}$, we derive the following linear equation that can be overdetermined because of $m\geq 1$:

[TABLE]

where $\mathsf{H}\in{\mathbb{R}}^{(m\cdot N)\times N}$ and we call it vector Hankel matrix given by

[TABLE]

Now, by assuming $\mathsf{H}$ to be of full column rank, the true solution of the target coefficients is explicitly expressed as follows:

[TABLE]

If this is not the case, use the Moore-Penrose pseudo-inverse of $\mathsf{H}$. 2. (ii)

Find the $N$ roots $\hat{\lambda}_{j}$ of the Prony polynomial. This is equivalent to locating the $N$ eigenvalues of the companion matrix $\mathsf{C}^{\prime}_{N}$:

[TABLE] 3. (iii)

and (iv)  Same as (iii) and (iv) of the Arnoldi-type algorithm.

Here we provide a geometrical interpretation of vector Prony analysis [69]. By following the argument in the beginning of Section 2.3.2, define the function space $\mathcal{K}^{\prime}_{N}$ generated by the operation of time-shifting as

[TABLE]

Then, by letting ${\mathbf{P}}^{\prime}_{N}:{\mathcal{F}}^{m}\to\mathcal{K}^{\prime}_{N}$ be a projection from the original space ${\mathcal{F}}^{m}$ onto $\mathcal{K}^{\prime}_{N}$, the following operator

[TABLE]

is still finite-dimensional and has a matrix representation with the same dimension of $\mathsf{C}_{N}$ in the \{[(\mathbf{\underline{U}}^{k}\mbox{\boldmathf})^{\top}, \ldots,(\mathbf{\underline{U}}^{k+N-1}\mbox{\boldmathf})^{\top}]^{\top}\}-basis, which we have denoted by $\mathsf{C}^{\prime}_{N}:{\mathbb{C}}^{N}\to{\mathbb{C}}^{N}$. Needless to say, we have implicitly assumed that the $N$ new functions (“fictive” observables) are linearly independent; this would be more relevant than in $\mathcal{K}_{N}$ with a single \mathbf{\underline{U}}^{k}\mbox{\boldmathf}. Now, denoting the vector Hankel matrix $\mathsf{H}$ by the new Krylov subspace $\mathsf{K}^{\prime}_{N}$ and introducing \mbox{\boldmathr}\in{\mathbb{R}}^{m\cdot N} as a residual with \mbox{\boldmathr}\bot\,\mathrm{span}\{\mathsf{K}^{\prime}_{N}\}, we rewrite (19) as follows:

[TABLE]

In analogy with (12) in Arnoldi-type algorithm, we are able to state that in vector Prony analysis, $p$ is chosen to minimize the ${\mathbb{C}}^{m\cdot N}$-norm of $r$ in (20), corresponding to choosing the projection operator $\mathbf{P}^{\prime}_{N}$ so that \mathbf{P}^{\prime}_{N}\mathbf{\underline{U}}^{N}\mbox{\boldmathf} is the least-squares approximation to \mathbf{\underline{U}}^{N}\mbox{\boldmathf} at \mbox{\boldmathx}_{0}\in\mathbb{X} satisfying \mbox{\boldmathy}_{0}=\mbox{\boldmathf}(\mbox{\boldmathx}_{0}). By comparison with $\mathsf{K}_{N}\in{\mathbb{R}}^{m\times N}$, the rank-deficiency of $\mathsf{K}^{\prime}_{N}\in{\mathbb{R}}^{(m\cdot N)\times N}$ when choosing $p$ can be avoided even in the case of data with low spatial dimension ($m\ll N$). Thus, the unique solution of $p$, in other words, an unique decomposition of \mbox{\boldmathy}_{k} can be derived. The possibility of avoiding the rank-deficiency is a result of time-shifting in which we are able to utilize more data for computing $\mathsf{C}^{\prime}_{N}$. The matrix $\mathsf{C}^{\prime}_{N}$ is thought of as an approximation to the action of the Koopman operator on the fictive finite-dimensional space $\mathcal{K}^{\prime}_{N}$.

In summary, we state mathematically that the difference of the Arnoldi-type algorithm and vector Prony analysis for KMD is the choice of finite-dimensional function spaces on which the original Koopman operator is projected. Both the derived matrices are thought of as the least-squares approximations of the corresponding action of the Koopman operator. Technically, the difference is just the number of snapshots necessary for computation: to obtain $N$ pairs of KMs and KEs, $N+1$ snapshots are needed for Arnoldi-type algorithm, while $2N$ snapshots for vector Prony analysis. A comparison of the two algorithms using realistic power system data is presented in [69].

2.3.3 Coherency notion

We define the notion of coherency in the context of KM. The case of oscillatory KM, in which the KE has a non-zero imaginary part, is now addressed because the study on coherency identification in power systems normally deals with oscillatory responses following a disturbance. For an oscillatory pair of KMs \mbox{\boldmathV}_{i_{1}} and \mbox{\boldmathV}_{i_{2}}, called the KM pair $\{i_{1},i_{2}\}$, with KEs $\lambda_{i_{1}}=|\lambda_{i}|{\rm e}^{{\rm i}2\pi\nu_{i_{1}}}$ and $\lambda_{i_{2}}=\lambda^{\ast}_{i_{1}}=|\lambda_{i_{1}}|{\rm e}^{-{\rm i}2\pi\nu_{i_{1}}}$, the associated modal dynamics are defined as

[TABLE]

where we have used the amplitude coefficient $A_{ij}$ and initial phase $\alpha_{ij}$, defined as

[TABLE]

Thus, we can say that a set of observables $\mathbb{I}\subseteq\{1,\ldots,m\}$ is coherent with respect to the KM pair $\{i_{1},i_{2}\}$ if the amplitude coefficients $A_{i_{1}j}$ are the same for all $j\in\mathbb{I}$, and the initial phases $\alpha_{i_{1}j}$ are also the same. It is noted that the definition is strict compared with the definitions of slow-coherency [70, 71] and near-coherency [72], because it does not admit any finite, constant phase difference of oscillations along the system’s dynamics. A relaxed notion of the coherency is introduced in [25, 49, 45] with a small positive parameter $\epsilon$: if for the two observables $f_{j_{1}}$ and $f_{j_{2}}$ in $f$, the initial phases satisfy $|\alpha_{i_{1}j_{1}}-\alpha_{i_{1}j_{2}}|<\epsilon$, then the two observables are called $\epsilon$-phase-coherent with respect to the KM pair $\{i_{1},i_{2}\}$. Below, we will sometime denote the KM pair $\{i_{1},i_{2}\}$ by $i_{1}$ for simplicity of the notation.

2.4 Model-order reduction

By starting at [34], we look at a connection between spectrum of the Koopman operator and canonical equations of motion, and we then provide an idea on model-order reduction of nonlinear systems that is partly demonstrated in [43] and Section 4.

2.4.1 Hamilton’s canonical form and Koopman eigenfuctions

First of all, let us review the beginning of Hamiltonian mechanics based on [73]. The so-called Hamilton’s equation of $n$-degree-of-freedom motions is the following: for displacement variables \mbox{\boldmathq}=[q_{1},\ldots,q_{n}]^{\top}\in\mathbb{M} (configuration manifold) and momentum variables \mbox{\boldmathp}=[p_{1},\ldots,p_{n}]^{\top}\in\mathrm{T}^{\ast}_{q}\mathbb{M} (cotangent space of $\mathbb{M}$ at $q$),

[TABLE]

where H:\mathrm{T}^{\ast}\mathbb{M}~{}\textrm{(cotangent bundle of \mathbb{M})}\to{\mathbb{R}} is the Hamiltonian that is normally related to energy of the system. If the dynamical system (26) is integrable [73], then the system is written as the famous form of action-angle variable representation: for angle variables \mbox{\boldmath\theta}=[\theta_{1},\ldots,\theta_{n}]^{\top}\in\mathbb{T}^{n} and action variables \mbox{\boldmathI}=[I_{1},\ldots,I_{n}]^{\top}\in{\mathbb{R}}^{n}_{\geq 0},

[TABLE]

where $\mathit{\Omega}_{j}:{\mathbb{R}}^{n}_{\geq 0}\to{\mathbb{R}}$ is the action-dependent function of angular frequency. In this, the $n$ action variables $I_{1},\ldots,I_{n}$ are called $n$ constants of motion because their time differentiations are equal to zero. The discrete-time version of (27) is given by

[TABLE]

Also, the following perturbed form of (28) has been widely analyzed for modal interaction and instabilities in nearly-integrable Hamiltonian systems:

[TABLE]

where $F_{j}$ and $G_{j}:\mathbb{T}^{n}\times{\mathbb{R}}^{n}_{\geq 0}\to{\mathbb{R}}$ are functions that represent (smooth) perturbations to the integrable motion.

Here we derive the same system as (28) for the general system (2) by introducing the Koopman eigenfunctions $\phi_{j}$ and associated eigenvalues $\lambda_{j}$. First, based on [34], let us consider the Koopman eigenfunctions as a transformation of the state $x$ like

[TABLE]

Then, we have

[TABLE]

and by taking $\theta_{j}=\mathrm{Arg}(z_{j})$ and $I_{j}=|z_{j}|$, we have

[TABLE]

This equation is equivalent to the action-angle representation (28) if $|\lambda_{j}|=1$ holds. The condition of KEs implies that the Koopman operator ${\mathbf{U}}$ is unitary, for which the state’s dynamics are on an attractor. Here we should emphasize that (31) is derived for the general (possibly dissipative) dynamical system (2) via the Koopman eigenfunctions. The vector form of the equation (30) is derived as follows: for \mbox{\boldmathz}:=[z_{1},z_{2},\ldots]^{\top},

[TABLE]

Note that the state space should be taken as $\ell_{2}$ or $\ell_{\infty}$ under a certain condition of the infinite-dimensional matrix $\mathsf{\Lambda}$.

Here, if the stable dynamics observed have only a finite number of peaks in the frequency domain, they would be dominated by KEs which magnitudes are less than or equal to unity. When we sort the KEs in the slowest direction, that is,

[TABLE]

the asymptotic component in the dynamics is dominated by several leading KEs, and the level sets of the associated Koopman eigenfunctions form sub-spaces in the state space, on which the asymptotic component is well captured; this concept is developed in [34] as isostable.

2.4.2 Koopman spectrum-preserving model-order reduction

The finding of canonical-like form of equations leads to an idea of model-order reduction of the discrete-time system (2). Now, let us suppose that the first $M$ variables $z_{1},\ldots,z_{M}$ in (32) are sufficient for capturing the state’s dynamics of interest. For this, let us introduce a finite-dimensional truncation of the first $M$ Koopman eigenvalues and eigenfunctions represented by a infinite-dimensional matrix with a finite $M$ rank,

[TABLE]

Thus, we have the following system to represent the state’s dynamics:

[TABLE]

Here, in order to focus on the time-varying part above, by defining \mbox{\boldmathz}^{M}:=[z_{1},\ldots,z_{M}]^{\top}\in{\mathbb{C}}^{M} and $\mathsf{\Lambda}_{M}:=\mathrm{diag}(\lambda_{1},\ldots,\lambda_{M})\in{\mathbb{C}}^{M\times M}$, we have

[TABLE]

For the $n$-dimensional state’s dynamics in (2), if $M\ll n$ holds, then (33) is a relevant reduced-order representation of the original system (2). Since the dominant part of the discrete spectrum is preserved in the procedure, we call it the Koopman spectrum-preserving model-order reduction.

The derivation of the reduced-order model (33) from (2) is possible by conducting both the KMD for finite data on the state’s dynamics and the Petrov-Galerkin projection method [74] applied to the original system (2). The following idea does not require explicit formula of the Koopman eigenfunctions $\phi_{j}$. Now, consider $N$ pairs of KEs and KMs computed from finite-time data on the state’s dynamics \{\mbox{\boldmathy}_{k}=\mbox{\boldmathx}_{k};k=0,1,\ldots\} as

[TABLE]

where $k=0,1,\ldots,N-1$. Suppose that we sort the $N$ pairs of KEs and KMs such that the first $M(\leq N)$ pairs are sufficient for capturing the finite-time dynamics of interest. Then, the following transformation of the state is defined:

[TABLE]

Thus, by substituting the above to the original system (2), we have

[TABLE]

and by assuming $\mathsf{\tilde{V}}^{\ast}\mathsf{\tilde{V}}\in{\mathbb{C}}^{M\times M}$ to be regular and using the standard argument of Petrov-Galerkin projection, we have

[TABLE]

Here we represent the $j$-th component of \mbox{\boldmath\tilde{z}}^{M} as $\tilde{I}_{j}{\rm e}^{{\rm i}\tilde{\theta}_{j}}$. Along the state’s dynamics used for the computation of KMs, without the truncation after the first $M$ KMs in (34), the variable $\tilde{I}_{j}$ does not change as time goes on, and $\tilde{\theta}_{j}$ rotates with a constant speed $\mathrm{Arg}(\tilde{\lambda}_{j})$ given by the KE $\tilde{\lambda}_{j}$. In this way, we can regard the right-hand side of (35) as another formulation of the reduced-order model (33) that is computed by a combination of the model (2) and finite-time data \{\mbox{\boldmathx}_{k}\}. Furthermore, we can re-write (35) as follows:

[TABLE]

where \mbox{\boldmath\tilde{\theta}}=[\tilde{\theta}_{1},\ldots,\tilde{\theta}_{M}]^{\top} and \mbox{\boldmath\tilde{I}}=[\tilde{I}_{1},\ldots,\tilde{I}_{M}]^{\top}, and we have

[TABLE]

or simply

[TABLE]

The right-hand sides of (36) are regarded as terms in an action-angle formulation of the reduced-order model (33) that are analogue to (29).

3 Coherency identification of coupled swing dynamics

This section provides the first application of KMD to power systems technology—identification of coherency of swing dynamics in coupled generator networks [42]. The so-called coherency identification is related to short-term swing or transient stability as an ability of a power system to maintain synchronism when subjected to a large disturbance. The identification is to find a group of synchronous generators swinging together with the in-phase motion and is a basis of development of reduced-order models, which is traditionally used to reduce computational effort and currently employ on-line stability assessment, and dynamical system analysis of loss of transient stability (see [75]). Many groups of researchers have developed methods for coherency identification: see [76, 70, 71, 77, 72, 78] and references of [42]. In this section, we apply the technique in Sections 2.3.2(a) and 2.3.3 to data on short-term swing dynamics in the New England 39-bus test system (NE system) [79] and demonstrate how the technique works for identification of modes and coherency.

3.1 New England test system and coupled swing dynamics

The NE system is a well-known benchmark model for power systems research shown in Fig. 1(a) and contains 10 generation units (equivalently 10 synchronous generators, circled numbers in the figure), 39 buses, and AC transmission lines. Most of the buses have constant active and reactive power loads. See [79] for the details of the system.

Here, we introduce the equations of motion of generators in the NE system. Assume that bus 39 is the infinite bus555A source of voltage constant in phase, magnitude, and frequency, and not affected by the amount of current withdrawn from it. The short-term swing dynamics of generators 2–10 are represented by the following nonlinear differential equations [80]:

[TABLE]

where the integer label $i=2,\ldots,10$ stands for generator $i$. The variable $\delta_{i}$ is the angular position of rotor in generator $i$ with respect to bus 1 and is in radians [rad]. The variable $\omega_{i}$ is the deviation of rotor speed in generator $i$ relative to that of bus 1 and is in radians per second [rad/s]. We set the variable $\delta_{1}$ to a constant, because bus 39 is assumed to be the infinite bus. The parameters $f_{\rm b}$, $H_{i}$, $D_{i}$, $P_{\mathrm{m}i}$, $E_{i}$, $G_{ii}$, $G_{ij}$, and $B_{ij}$ are in per unit system except for $H_{i}$ and $D_{i}$ in seconds [s], and for $f_{\rm b}$ in Hertz [Hz]. The mechanical input power $P_{\mathrm{m}i}$ to generator $i$ and its internal voltage $E_{i}$ are normally constant in the short-term regime [80]. The parameter $H_{i}$ is the per unit time inertia constant of generator $i$, and $D_{i}$ its damping coefficient. The parameter $G_{ii}$ is the internal conductance, and $G_{ij}+{\rm i}B_{ij}$ is the transfer impedance between generators $i$ and $j$. Electrical loads are modeled as constant impedances.

Now, we present an example of data on short-term swing dynamics in the NE system. Fig. 1(b) shows the time responses of rotor speed deviations $\omega_{i}$ under the initial condition:

[TABLE]

where $\delta_{i}^{\ast}$ represents the value of stable equilibrium point of (38). The initial condition physically corresponds to a local disturbance at generator 8. The details of the numerical setting are in [42]. The generators do not show any divergence motion in the figure, that is, they do not show any loss of transient stability for the selected disturbance. Generators 2, 3, 6, and 7 show coherent swings excited by the local disturbance. We call these generators the coherent group. The other generators show incoherent swings in the figure. Generator 9 shows swings similar in frequency and phase to the coherent group, but the amplitude of swings is a little larger. Generators 8 and 10 have swings of larger amplitudes than the others, because the initial condition is localized at generator 8, and the two generators are electrically close.

3.2 Coherency identification via Koopman mode decomposition

We compute the KMs and KEs for the data on coupled swing dynamics shown in Fig. 1(b). The computation is done with the two different algorithms: Fourier-based formula (9) and Arnoldi-type algorithm in Section 2.3.2(a). For computation we need to choose the observable \mbox{\boldmathf}(\mbox{\boldmath\delta},\mbox{\boldmath\omega}) where \mbox{\boldmath\delta}=(\delta_{2},\ldots,\delta_{10})^{\top} and \mbox{\boldmath\omega}=(\omega_{2},\ldots,\omega_{10})^{\top}. In this section, we use the variables of rotor speed deviations, $\omega$, as the observable: \mbox{\boldmathf}(\mbox{\boldmath\delta},\mbox{\boldmath\omega})=\mbox{\boldmath\omega}. This observable has a clear physical meaning in power systems: one measures rotor speeds or frequencies for every generation plant. We use the simulation output shown in Fig. 1(b) that extracts \{\mbox{\boldmath\omega}(nT)\}_{n=0}^{N}, where the uniform sampling period $T=1/(50{\mu\Omega\rm\,Hz})$ and the number of snapshots $N+1=1001$.

A Fourier analysis of the waveforms in Fig. 1(b) shows that generators 2, 3, 6, and 7 have the similar shape of spectrum with peak frequencies 0.35 Hz and 1.10 Hz, and that generators 8 and 10 have the peak frequencies 1.10 Hz and 1.30 Hz. Hence we compute the terms including KMs, \phi_{i}(\mbox{\boldmathx}[0])\mbox{\boldmathV}_{i}, using the projection operator ${\mathbf{P}}^{\nu}$ with $\nu=(0.35{\mu\Omega\rm\,Hz})T$, $(1.10{\mu\Omega\rm\,Hz})T$, and $(1.30{\mu\Omega\rm\,Hz})T$. We use the finite-time approximation of (9) from $k=0$ to $N$. The numerical results are shown in Fig. 2. The amplitude coefficients and initial phases, which are defined in (25), are also shown. For 0.35 Hz, the values of amplitude coefficients are close for each of the generators, and their initial phases are also close except for generator 8. The generators except for 8 hence show in-phase swings with 0.35 Hz. For 1.10 Hz, the values of amplitude coefficients and initial phases are close for generators 2, 3, 6, 7, and 9, and hence they show in-phase swings with 1.10 Hz. The two KMs with 0.35 Hz and 1.10 Hz capture the coherent motion of generators 2, 3, 6, 7, and 9. For 1.10 Hz and 1.30 Hz, the amplitude coefficients for generators 8 and 10 are larger than the others. These two KMs capture the large swings of generators 8 and 10 observed in Fig. 1(b). Thus, we can extract spatial modes oscillating with a single frequency directly from the data on simulation output of the NE system.

Next, we compute the KEs and KMs (empirical Ritz values $\tilde{\lambda}_{j}$ and associated vectors \tilde{\mbox{\boldmathV}}_{j}) using the Arnoldi-type algorithm. Fig. 3(a) shows the empirical Ritz values $\tilde{\lambda}_{j}$. The norm of $r$ in Step (i) was of order $10^{-12}$. Many KMs are obtained and are close to the unit circle $|\tilde{\lambda}_{j}|=1$. Now let us focus on KMs that have both large growth rates $|\tilde{\lambda}_{j}|$ and large norms of \tilde{\mbox{\boldmathV}}_{j}. Such modes represent sustained swing components for the time duration of data and have dominant magnitudes in the time series. Tab. 1 shows numerical results on KEs and KMs, which we call KM pair 1 to 10. The order of KM pairs in Tab. 1 is based on the magnitudes of growth rates. We now pick up KM pair 7 to 9 that have large norms in the table. KM pair 1 and 2 have large norms, too. But, their frequencies are close to KM pair 7 and 8, respectively. Fig. 3(b) shows the KMs \tilde{\mbox{\boldmathV}}_{j} for Mode $j$ ($j=7,8,9$). The amplitude coefficients $A_{ji}$ ($j=7,8,9,i=2,\ldots,10$) and initial phases $\alpha_{ji}$ are also shown. We display the results in order to make it easy to compare this with the Fourier-based results in Fig. 2. For example, the frequency 1.3078 Hz of KM pair 7 is close to one of the dominant frequencies for generators 8 and 10, that is, 1.30 Hz. In fact, since the values of amplitude coefficients $A_{7,i}$ in Fig. 3(b) are large for generators 8 and 10, the corresponding modal dynamics are localized at these generators. Thus, we can decompose the coupled swing dynamics in the NE system into a set of KMs, namely, spatial modes of oscillation with single frequency.

The decomposition into KMs makes it possible to extract coherent generators in the coupled swing dynamics. In fact, the two KM pairs 8 and 9 capture a coherent motion related to the coherent group of generators. These frequencies, 1.0962 Hz and 0.3727 Hz, are close to the frequencies of the coherent group, 1.10 Hz and 0.35 Hz. For KM pair 8, the values of $A_{8,i}$ are close for each of generators 2, 3, 6, 7, and 9, and their initial phases $\alpha_{8,i}$ are also close. The distribution of $(A_{8,i},\alpha_{8,i})$ is plotted in Fig. 4. The points ($\times$) for generators 2, 3, 6, 7, and 9 are clustered around the coordinate $(1.2,3\pi/4)$. Generators 2, 3, 6, 7, and 9 hence show in-phase swings with 1.0962 Hz. For KM pair 9, the values of $A_{9,i}$ and $\alpha_{9,i}$ are close for each of the generators except for generator 8, and hence they show in-phase swings with 0.3727 Hz. The corresponding plots ($\circ$) in Fig. 4 are clustered except for generator 8. The two KMs capture the coherent motion of generators 2, 3, 6, 7, and 9. In this way, we can identify the coherent group observed in Fig. 1(b) by using the decomposition into KMs. The plot as Fig. 4 provides a systematic way to identify coherent swings and generators by using an automatic clustering algorithm.

In this section, we used the Fourier-based formula (9) and Arnoldi-type algorithm in Section 2.3.2(a) for cross-check of the computation of KMs and KEs. The results for amplitude coefficients and initial phases between Fig. 2 and Fig. 3(b) with KM pair 7 to 9 are compared. The results for amplitude coefficients are qualitatively similar. However, the results for initial phases are somewhat different, especially at generator 8. Also, for the amplitude coefficients, there are quantitative differences between the results in Figs. 2 and 3(b). These differences might be due to the fact that the Fourier-based formula (9) assumes the unitary case of the Koopman operator, that is, the dynamics are on an attractor. Indeed, the current analysis is performed for off-attractor dynamics. However, identification of coherent swings and generators is possible using both the algorithms. In fact, the coherent swings of generators 2, 3, 6, 7, and 9 are captured well.

4 Precursor diagnostic of coupled swing instabilities

This section studies a precursor to phenomena of loss of stability in the coupled swing dynamics [43]. Precursor-based monitoring of instabilities has attracted a lot of interest in terms of prediction: see [81, 82, 83] and references of [43]. The precursor in this section is based on the discovery of [75], an emergent transmission path of energy from many oscillatory modes to one oscillatory mode that represents an instability phenomenon of interest. The pathway from high frequency modes to the lowest frequency mode is called the Coherent Swing Instability (CSI). The modes can be extracted from sensor data or data provided by simulation outputs of power system oscillations via KMD. The CSI transmission path is identified by computation of the so-called action transfer operator that is derived by refining the mathematical model (38) using the technique in Section 2.4. This provides a new technique for monitoring the instability by a combination of practical data, mathematical modeling, and computation.

4.1 Derivation of action transfer operator

The authors of [84] gave an insight into the transmission paths of energy in a nearly-integrable Hamiltonian system by linearizing a perturbed action-angle representation. The linearized system provides the time-dependent operator as the system evolves. In this section, we will term this operator the action transfer operator that can quantify the change of action variables. Now, we apply the approach developed in [84] to the discrete-time system (37) and show that by means of KM, the action transfer operator is defined for arbitrary dynamical systems. The linearizd system of (37) around a state (\mbox{\boldmath\tilde{\theta}}^{\ast},\mbox{\boldmath\tilde{I}}^{\ast}) is represented as

[TABLE]

Since the action variables are often associated with amount of energy contained in every KM, the sub-matrix \mathrm{D}_{\tilde{I}}\mbox{\boldmath\tilde{G}}(\mbox{\boldmath\tilde{\theta}}^{\ast},\mbox{\boldmath\tilde{I}}^{\ast}) provides a quantitative index of the model interaction via energy. Thus, we can define the Jacobian matrix \mathsf{J}(\mbox{\boldmathx}[k]) that describes the infinitesimal change of \mbox{\boldmath\tilde{I}}[k] at a state \mbox{\boldmathx}[k]:

[TABLE]

The right-hand side of (40) can be computed numerically: see Appendix of [43] in detail. The Jacobian matrix is a time-dependent operator and quantifies the change of action variables for the $M$ KMs. We term the Jacobian matrix the action transfer operator for the $M$ KMs. In this section, by using a concrete example of power system analysis, we will illustrate that the action transfer operator works for identification of a precursor to the instability.

4.2 Coherent swing instability of New England test system

We present other examples of coupled swing dynamics in the NE system. Fig. 5 shows the time responses of angular positions $\delta_{i}$ of 9 generators under the two initial conditions:

[TABLE]

and

[TABLE]

The initial conditions are close to (39) and physically correspond to local disturbances at generator 8. In Fig. 5(a) the generators do not show divergence motion, that is, they do not show any loss of transient stability for the selected disturbance. The behavior is close to that in Fig. 1(b). In Fig. 5(b) they are bounded during the period from $t=0{\mu\Omega\rm\,s}$ to $7{\mu\Omega\rm\,s}$ and then begin to grow coherently. Every generator loses synchronism with the infinite bus at the same time. This corresponds to the growth of amplitude of inter-area mode oscillation between the NE system and the infinite bus, namely, the outside of the system. This is typical of the CSI phenomenon.

In [75] we showed that CSI involves the divergent motion in the projection of full-system dynamics onto the plane of collective variables. The collective variables correspond to the well-known COA (Center-Of-Angle) coordinates in power systems analysis [85, 79]. For the NE system, the COA $\delta_{\rm COA}$ and its time derivative $\omega_{\rm COA}$ are defined as

[TABLE]

where $H:=\sum^{10}_{i=2}H_{i}$. The variables $\delta_{\rm COA}$ and $\omega_{\rm COA}$ describe the averaged motion of all the generators in the system. Fig. 6 plots the trajectories of (38) shown in Fig. 5 on $\delta_{\rm COA}$–$\omega_{\rm COA}$ plane. The trajectories start near the origin at time $0{\mu\Omega\rm\,s}$ and make a couple of almost periodic loops around the initial point. In Fig. 6(a) the trajectory keeps executing bounded loops. In Fig. 6(b) it escapes and finally diverges.

4.3 Computation of Koopman modes

We compute the KMs for the coupled swings shown in Fig. 5(a). The computation of KMs is done with the Fourier-based formula (9). The state $x$ itself is chosen as the observable, namely \mbox{\boldmathf}(\mbox{\boldmathx})=\mbox{\boldmathx}. We use the simulation output shown in Fig. 5(a) that extracts the data \{\mbox{\boldmathx}(nT)\}_{n=0}^{N}, where the uniform sampling period $T=1/(50{\mu\Omega\rm\,Hz})$ and the number of snapshots $N+1=1001$. Also, in order to use (9) for computing KMs, we need to identify dominant frequencies in the course of coupled swings. A Fourier analysis shows that the time responses in Fig. 5(a) have clean peaks at the five frequencies, $0.40{\mu\Omega\rm\,Hz}$, $0.90{\mu\Omega\rm\,Hz}$, $1.00{\mu\Omega\rm\,Hz}$, $1.25{\mu\Omega\rm\,Hz}$, and $1.45{\mu\Omega\rm\,Hz}$.

Now we compute the terms \phi_{i}(\mbox{\boldmathx}[0])\mbox{\boldmathV}_{i} using the projection operator ${\mathbf{P}}^{\nu}$ with the five frequencies: $i=1$ for 0.40 Hz, 2 for 0.90 Hz, 3 for 1.00 Hz, 4 for 1.25 Hz, and 5 for 1.45 Hz. The numerical results are shown in Fig. 7. The angular positions $\delta_{i}$ only are displayed. The amplitude coefficients and initial phases, which are defined in (25), are also shown. For 0.40 Hz, the values of amplitude coefficients are close for each of the generators, and their initial phases are also close. All the swings of the 9 generators are hence said to be coherent with respect to the KM with 0.40 Hz. We call it the coherent KM in the same reason as in Section 3.

4.4 Precursor diagnostic to coherent swing instability

We compute the action transfer operator ${\sf J}$ for the CSI phenomenon shown in Fig. 5(b). The KMs in Fig. 7 are computed for the bounded swings in Fig. 5(a). The results on Figs. 5(a) and (b) are obtained with different initial conditions. Hence, we consider how ${\sf J}$ behaves under dynamics perturbed by a slight change of initial conditions. We use the tools (model-order reduction and action-angle representation) developed in Section 2.4.

Figure 9 shows several snapshots of the action transfer operator ${\sf J}$ along the COA trajectory shown in Fig. 6(b). In the snapshots, the $i$-th KM ($i=1,\ldots,5$) corresponds to the $2i-1$ and $2i$ rows. The two variables $z_{2i-1}=I_{2i-1}{\rm e}^{{\rm i}\theta_{2i-1}}$ and $z_{2i}=I_{2i}{\rm e}^{{\rm i}\theta_{2i}}$ are assigned to the $i$-th KM. By definition, the action variables $I_{2i-1}$ and $I_{2i}$ have the same behaviors, and the angle variables $\theta_{2i-1}$ and $-\theta_{2i}$ also have the same. Because both the variables $I_{2i-1}$ and $I_{2i}$ affect the infinitesimal change of $I_{j}$ ($j=1,\ldots,10$), in $\mathsf{J}$ we need to plot $\partial\tilde{G}_{j}/\partial I_{2i-1}$ as well as $\partial\tilde{G}_{j}/\partial I_{2i}$.666See Fig. 9. For the 1st KM ($j=1$), in snapshot 2 the variable $I_{9}$ for the 5th KM dominantly affects the time change of $I_{1}$. On the other hand, in snapshot 4 the variable $I_{10}$ dominantly affects $I_{1}$. These observations both are regarded as the effect of the 5th KM to the 1st KM. First, the magnitude of elements in the upper triangular part change as time goes on. On the other hand, the magnitude of elements in the lower triangular part does not change and is almost constant. This salient feature of $\mathsf{J}$ suggests a directed transmission path of energy from high-frequency KMs to low ones. This is consistent with that in [84] that analyzes global dynamics occurring in nearly-integrable Hamiltonian systems. The current result is obtained for a set of KMs embedded in dynamics of a general nonlinear power system model. Second, we see that the transmission path of energy from high-frequency KMs to the lowest one becomes stronger when the COA trajectory largely grows in the neighborhood of the $\delta_{\rm COA}$-axis: see the snapshots at time (4) 7.2 s and (6) 9.0 s in Fig. 9. The transmission path appears several times, and then the magnitude $|z_{1}|$ of the lowest KM begins to diverge as shown in Fig. 9, while the magnitudes of the other KMs are still bounded. Since the lowest KM with 0.40 Hz is the coherent KM, the divergence of magnitude $|z_{1}|$ implies CSI. This provides a precursor diagnostic for CSI by a combination of modeling and data: the emergent transmission paths of energy from high-frequency KMs to the lowest coherent KM.

5 Stability assessment of power systems without models

This section provides stability assessment of power systems without any development of models [44]. Many methods for stability assessment of power systems have been developed [80]. Traditional methods are mainly model-based, that is, the stability assessment is performed by investigating a mathematical model that represents the target dynamics of a power system. It is widely recognized that cascading outages are fairly complicated emergent phenomena in the high-dimensional nonlinear systems. Therefore, it is difficult to obtain a mathematical model that explains all events and time evolution of a cascading outage. Even if we could obtain such a model, it would be not easy to gain a dynamical insight to the cascading outages from the model because of its complexity. Also, the large increase of renewable sources makes it hard to obtain a relevant deterministic model because of their uncertain nature. Thus, in contrast to the model-based approach, it is necessary to develop methods that indicate spatio-temporal structure of instabilities and their precursors from data. Several reports exist on this line of research [86, 87, 88]. In this section, we present an approach to stability assessment based on measured physical power flow data via KMD. In contrast to the existing methods, our approach provides not only dynamic patterns of power flows, which we refer to as the base flow patterns below, but also stability information.

5.1 Mathematical formulation

Consider the finite-time data on dynamics of physical power flows (more precisely, active power flows) under uniform sampling, given by

[TABLE]

where \mbox{\boldmathP}_{k}\in{\mathbb{R}}^{m} is the snapshot of power flows at the discrete time $k$, $m$ the number of measurement sites (for example, generation plants, transmission lines, and substations), and $N$ the number of available snapshots. Generally speaking, the power flows are determined by the internal states of a power system such as rotating frequencies and voltages of AC generators in power plants, bus voltages in substations, power consumptions in loads, and states of controllers in plants and substations [80]. Here we use $x$ to represent all the internal states belonging to a space ${\mathbb{X}}$ and assume that the time evolution of $x$ is represented by the deterministic system as follows:

[TABLE]

where \mbox{\boldmathT}:{\mathbb{X}}\to{\mathbb{X}} is a nonlinear vector-valued map and constructed via power system modeling techniques (see, e.g., [80]); thus the dimension of the map is possibly very large. Now, in order to associate the finite-time data with the internal state’s dynamics, we define a vector-valued observable defined on ${\mathbb{X}}$. For the current analysis, it is reasonable to define the observable as a map \mbox{\boldmathf}_{\mathrm{P}}:{\mathbb{X}}\to{\mathbb{R}}^{m} such that the snapshot \mbox{\boldmathP}_{k} of power flows at time $k$ is written as

[TABLE]

The idea which we now propose for stability assessment is straightforward from KMD. Here we develop the analysis for the case when there are unstable eigenvalues of the Koopman operator, but we have access to data for only a short period of time. If each of the $m$ components of \mbox{\boldmathf}_{\mathrm{P}} lies within the span of eigenfunctions $\phi_{j}$, then vector-valued function \mbox{\boldmathf}_{\mathrm{P}} is expanded in terms of these eigenfuntions as

[TABLE]

where \mbox{\boldmathV}_{j}\in{\mathbb{C}}^{m} are regarded as vector coefficients in the expansion. Thus, the time evolution \mbox{\boldmathP}_{k}=\mbox{\boldmathf}_{\mathrm{P}}(\mbox{\boldmathx}_{k}) starting at \mbox{\boldmathP}_{0}=\mbox{\boldmathf}_{\mathrm{P}}(\mbox{\boldmathx}[0]) is represented as

[TABLE]

where \mbox{\boldmathV}_{j} is the $j$-th KM associated with KE $\lambda_{j}$. An unstable KE is an eigenvalue whose modulus is larger than one, corresponding to the KM that grows exponentially as time increases. Here we propose the following method for detection of power system instability: if for a given set of data on power flows obtained over a time interval, there exists no unstable KE, then we conclude that the associated power system is stable. If this is not the case, then the power system behaves in an unstable manner.

Now we discuss the notion of stability addressed in this section. Our analysis is based on data for only a short period of time, in other words, is intended to describe short-term dynamics exhibited by a nonlinear dynamical system. In this sense, the current notion of stability is different from the conventional one that is intended to long-term, asymptotic dynamics. In a dynamical viewpoint, the current notion deals with a portion of a trajectory and is formulated by examining the unstable eigenvalues of the finite-time Koopman operator corresponding to the time period over which the trajectory is executed. The difference of notions of stability can be seen by considering a periodic orbit that slows down as it passes near a saddle and then speeds up. Portions of such an orbit could be judged stable or unstable by the method of the section depending on which portion is chosen. In this section, we intend to provide a new method for detecting such a short-term trend in dynamics from observational data, thereby obtaining an insight to control measures for (nearly) real-time management.

We use the above approach to detecting instabilities based on KMD to analyze data on power flows sampled from the 2006 system disturbance of the European interconnected grid. Our main idea of the application is to decompose the data on power flows into a set of KMs, implying base flow patterns. A base flow pattern is geographically distributed over the system, exhibits power flows that grow possibly, and is regarded as a coherent spatial unit of power flows. Let \mbox{\boldmath\tilde{V}}_{j}=[\tilde{V}_{j1},\tilde{V}_{j2},\ldots,\tilde{V}_{jm}]^{\top} denote the $j$-th empirical Ritz vector associated with the possibly complex-valued Ritz value $\tilde{\lambda}_{j}=\tilde{r}_{j}{\rm e}^{{\rm i}2\pi\tilde{\nu}_{j}}$ for the finite-time data. For convenience we will call the empirical Ritz vector \mbox{\boldmath\tilde{V}}_{j} the KM, although it differs from it possibly by a multiplicative constant. Since the original data are real-valued, if the KM \mbox{\boldmath\tilde{V}}_{j} is complex-valued, then there exists the conjugate KM \mbox{\boldmath\tilde{V}}_{j}^{\mathrm{c}} with the conjugate Ritz value $\tilde{\lambda}_{j}^{\ast}$ in the expansion (45). Here we suppose that the $(j+1)$-th KM \mbox{\boldmath\tilde{V}}_{j+1} corresponds to the conjugate vector \mbox{\boldmath\tilde{V}}_{j}^{\mathrm{c}} and the $(j+1)$-th Ritz value $\tilde{\lambda}_{j+1}$ to $\tilde{\lambda}^{\ast}_{j}$. In this way, the dynamics of base flow pattern described by the KM pair $\{i,i+1\}$ are given in the same manner as (24):

[TABLE]

where $k=0,1,\ldots,N-2$. The spatial shape of the dynamics is captured by the modulus and argument of the KM. In the following, we investigate the dynamics of power flows by directly computing the KMs from measured data.

5.2 Demonstration in 2006 European grid-wide disturbance

This disturbance was recorded in the night of November 4, 2006 and affected the UCTE’s (Union for the Co-ordination of the Transmission of Electricity’s) synchronously European interconnected grid [7]. It begun in the North German transmission grid, led to a splitting of the UCTE grid into three areas, and caused an interruption of power supply for more than 15 millions European households. This splitting is due to the tripping of interconnected AC lines and implies that the three areas were asynchronously operated. In each of the areas a significant imbalance of power emerged.

The data on dynamics of power flows in the UCTE grid are introduced in Fig. 10A, where we show the dynamics of power exchange deviations from 21:50 till 22:10. A power exchange deviation is defined as the difference between the actual power exchanges and the scheduled power exchanges with other countries. The areas focused in this figure are Belgium, Switzerland, D-RWE, Spain, France, Italy, Netherland, and Portugal. The data in Fig. 10A were obtained by sampling values from the original figure in Appendix 7 of [7] under the uniform sampling (period $T=30\,\mathrm{s}$). The final report [7] states, “These transient deviations are the result of a global shift in physical power flows within the UCTE synchronous grid area to changes in generation programs and exchange programs around 22:00. These curves show a situation which is rather normal and typical at that time.” However, the next single operation (coupling of busbars in the Landesbergen substation in Germany at 22:10) after the transient deviations in Fig. 10A initiated the widespread disturbance in the entire European grid. The operation of coupling the busbars was intended to resolve the heavy power flow on the Landesbergen-Wehrendorf 380 kV line in Germany, but it resulted in its immediate trip and initiated a cascade of line trips, starting with the 220 kV Bielefeld/Ost-Spexard line and continuing within the E.ON Netz grid [7]. Therefore, the N-1 criterion777The N-1 criterion is a basic principle in power system operation [7]. This rule is that any single loss of transmission or generation element should not jeopardize the secure operation of the interconnected network, that is, trigger a cascade of line trippings or the loss of a significant amount of consumption. of secure operation of the grid was not fulfilled before the busbar coupling [7].

We performed the KMD of the data on power exchange deviations presented in Fig. 10A. The data are the forty-one ($N=41$) snapshots \{\mbox{\boldmathP}_{0},\mbox{\boldmathP}_{1},\ldots,\mbox{\boldmathP}_{40}\} at eight different $(m=8)$ countries. The forty pairs of KEs and KMs $\tilde{\lambda}_{j}$ and \mbox{\boldmath\tilde{V}}_{j} were computed and shown in Figs. 10B and 10C. In Fig. 10B, we see that one conjugate pair of $\tilde{\lambda}_{j}$ exists outside of the unit circle, and the corresponding vector represents the unstable mode. The other $\tilde{\lambda}_{j}$’s exist inside of the unit circle, that is, they are associated with the stable modes. The modulus vectors $[|\tilde{V}_{j1}|,|\tilde{V}_{j2}|,\ldots,|\tilde{V}_{j8}|]^{\top}$ and argument vectors $[\mathrm{Arg}(\tilde{V}_{j1}),\mathrm{Arg}(\tilde{V}_{j2}),\ldots,\mathrm{Arg}(\tilde{V}_{j8})]^{\top}$ of the unstable mode are shown in Fig. 10C. Also, the dynamics of the unstable base flow pattern are presented in Fig. 10D. The period of the unstable mode is 37 minutes, and thus its oscillation is relatively slow. The spatial shape of the mode indicates that there are four large magnitude countries involved.

The computation shown in Fig. 10B reflects the KEs obtained during the computation over all 41 snapshots. It is of interest to consider the evolution of the KE computation over time, to understand the persistence of the unstable KM in relatively noisy data. We computed the KEs (and the associated KMs) starting from the 1st sample and computing until $N=8,9,\ldots,40$, namely, for many sets of snapshots with different final number $N$. The result is shown in Fig. 11 as the number density plot, where the color changes from blue to black and depicts the number density of unstable KEs in each cell. In this figure, the black cells are close to the unstable KEs shown in Fig. 10. Therefore, the unstable KM persists under a change of the sample length. As reviewed above, the dynamics of power exchange deviations are closely related to the initiation of the system disturbance. The base flow pattern of the unstable KM presented in Fig. 10D vividly shows how the large-scale unstable dynamics can be extracted from the noisy data in the report—the signature of the instability, in other words, the breaking of the N-1 criterion is clearly present in the dynamics even before time 22:00. It should be recalled that the UCTE report argues that the situation is rather normal and typical at that time.

6 Conclusion and future problems

In this paper, we have outlined our recent efforts on applications of the Koopman operator theory to power systems technology. The first part of this paper (Section 2) was devoted to an introduction to the Koopman operator theory for nonlinear dynamical systems, with emphasis on modal decomposition and computation. The second part (Sections 3 to 5) was devoted to the three topics of the applications: coherency identification of swings in coupled synchronous generators, precursor diagnostic of instabilities in the coupled swing dynamics, and stability assessment of power systems without any development of models. We emphasize that these are established as data-centric methods, namely, how to use massive quantities of data obtained numerically and experimentally, via the KMD (Koopman Mode Decomposition).

Here we note a topic excluded due to space constraints of this paper. In [45, 46] the authors demonstrate the KMD’s capability for partitioning a power network, which determines the points of separation in an islanding strategy. A practical data-driven method incorporating the KMD for network partitioning is proposed. It is shown that the data-driven method can cover network partitions derived from both spectral graph theory [89] and slow-coherency-based method [77]. The KMD-based partitioning pins down the frequency of oscillations as well as information on damping and participation, thus it can be used for monitoring and control purposes, namely the controlled islanding technique: see, e.g., [90].

In closing, we identify two future problems, not addressed in or following the current analysis in Sections 3 to 5:

• •

The applications presented here are specific to swing dynamics governing active power transmission and frequency maintenance. The power system has the other important physical aspect: reactive power and voltage. Loss of the normal level of voltage is often crucial in the course of cascading outages [1, 4, 5]. This poses a problem of characterization of non-oscillatory responses in terms of spectrum of the Koopman operator, which is still a difficult problem in the existing Koopman operator techniques.

• •

Control and operation of the power system in terms of the Koopman operator should be pursued. As stated in the beginning of Section 1, the dynamics of power systems occur in a wide range of scales in space and time. Depending on our purposes of control and operation, we need to extract dynamical information from the data collected practically or predicted in numerical simulations, and then to make a decision or actuation in order to archive a desirable physical response. KMD has two main advantages: (i) clear separation of spatio-temporal scales, even from relatively-noisy data as demonstrated in Section 5, and (ii) direct computation from data without any use of description of the underling system.

Acknowledgements

The first author appreciates Dr. Ryan Mohr, Mr. Hassan Arbabi, and Prof. Alexandre Mauroy for valuable discussions. He also thanks many collaborators and students, especially, Prof. Yutaka Ota, Mr. Yohei Kono, and Mr. Kyoichi Sako for valuable discussion on individual theoretical and applied topics. The work presented in this paper is supported in part by JSPS Postdoctoral Fellowships for Research Abroad, JST-CREST #JPMJCR15K3, and MEXT KAKEN #15H03964.