# Characterizing the load profile in power grids by Koopman mode decomposition
of interconnected dynamics

Ali Tavasoli, Heman Shakeri
School of Data Science, University of Virginia, Charlottesville, VA, USA,
22904

###### Abstract

Electricity load forecasting, which is essential for operating, monitoring,
and optimizing power grids, has seen different stages of development during
recent decades by various emerging statistical and deep learning approaches.
In this paper, we consider an interpretable machine learning approach by a
data-driven identification of the load dynamics in the operator-theoretic
framework. In particular, we represent the load data by the Koopman operator,
which is intrinsic to the underlying dynamics, and by computing the
corresponding eigenfunctions we decompose the load dynamics into a set of
coherent spatiotemporal patterns that are the most robust features of the
dynamics and each pattern evolves independently according to its single
frequency, making its predictability based on a linear dynamics. We stress on
the construction of the load dynamics based on coherent spatiotemporal
patterns intrinsic to the dynamics by showing that they are able to encode
rich dynamical features at multiple time scales regarding complex interactions
over interconnected power grids and different exogenous effects. For a more
efficient implementation of the Koopman operator approach, we cluster the load
data based on a modern kernel-based clustering approach, namely PHATE, and
discover the power stations with similar load patterns, particularly those
with synchronized dynamics. We examine our approach for a large-scale dataset
for a renewable electric power system within the continental European
electricity system and show the approach suggested based on Koopman
outperforms, in terms of accuracy and computational efficiency, a deep
learning (LSTM) architecture.

_Keywords_ : Power grids, load forecasting, Koopman operator, spectral
analysis, spatiotemporal patterns, clustering.

##  1 Introduction

Load forecasting approaches can be grouped into categories: physics-based and
data-driven modeling. The physics-based models are used to understand
different aspects that co-determine the spatiotemporal multilayer dynamics of
power systems, subject to complex behavior appearing as a result of the
interplay between dynamics of the nodes and topology, such as synchronization
and cascade effects [Filatrella2008physics; Motter2013nature;
Nardelli2014models]. On the other side, since the available physics-based
approaches are less reliable due to the assumed properties and environmental
parameters, and encouraged by the rapid development of sensor technologies and
the emergence of data mining techniques, the data-driven approaches are
increasingly drawing attention to reduce uncertainty and provide more robust
forecasts [Stankovic2021SINDy; Wang2022; Hu2022urban]. While long-term (over
three years) forecasts of the peak electricity demand are needed for capacity
planning and maintenance scheduling [Hyndman2010], very short-term (below one
day), short-term (a day to two weeks), and medium-term (two weeks to three
years) forecasts have recently gained more attention for power systems’ secure
operation and planning, control and scheduling, and commercial issues in
energy markets [Amjady2001; Song2005; Taylor2007; Taylor2008weather; Shu2012;
Charlton2014; Goude2014].

We can broadly classify the data-driven load forecasting methods into two
general categories, i.e., the conventional statistical methods and the
artificial intelligence methods. By conventional statistical models, time
series analysis is used to extract statistical information about the data and
potentially make predictions about future data points by assuming future data
has some relation with historically observed data. Statistical models mainly
include ordinary regression models [Weron2006book], auto-regressive moving
average model (ARMA) [Pappas2010ARMA], Autoregressive integrated moving
average model (ARIMA) [Chodakowska2021ARIMA], and exponential smoothing
[Taylor2007]. Since these regression models are generally for linear data,
methods such as random forest(RF) [Dudek2015RF], support Vector
Regression(SVR) [Chen2017SVR], gradient boost [Patnaik2021XGBoost], Gaussian
process [Shepero2018GP], or fuzzy systems [Efendi2015fuzzy] are used for
nonlinear problems in energy prediction and evaluation. However, because
machine learning approaches yield greater flexibility in constructing
forecasting models in the presence of features such as nonlinearity, high
volatility, and uncertainty of load profiles artificial intelligence models
have better forecasting accuracy than conventional statistical models in most
cases [Deb2017review; Bianchi2017book; Nti2020review].

The initial machine learning approaches for forecasting load were largely
based on Support Vector Machines (SVM) [Chen2004SVM] and basic (shallow)
Neural Networks (NNs) [Tsekouras2006NN; Sozen2007ANN; Geem2009ANN;
Ardakani2014NN], and later they developed into more advanced approaches such
as Deep NNs (DNNs) [Mohi2017DNN], Convolutional NNs (CNNs) [Marcus2018CNN],
and Long Short-Term Memory (LSTM) [Marino2016LSTM; Kong2019LSTM]. An issue
with deep architectures is overfitting, which limits the forecasting
performance, and rises by a casual selection of network parameters and layers.
Therefore, successful deep networks are implemented in more elaborated
structures, with the details and units that are used hierarchically for
different data processing levels [Heng2018pooling; Zhang2018hybrid;
Mousa2020hybrid; Dudek2021pattern; Yin2023photovoltaic].

While in modern machine learning frameworks, physical systems are studied by
the architectures that promote interpretable and generalizable models
[Brunton2016SINDy; Brunton_Kutz_2019], most neural network approaches are
opaque, defying straightforward analysis and interpretations
[Brunton&Kutz2022]. Moreover, while these networks are typically useful for
interpolation tasks, and memorizing the training data, they lack features in
extrapolation tasks and generalization beyond the training data. This limits
the application of NNs in physical and engineering systems, where, oftentimes,
a model extracted from limited training data must be valid outside the
training range to be used for prediction, optimization, and control.
Furthermore, for applications in power grids, where nonlinear interactions in
large-scale interconnected systems come with arbitrary size and structure, we
need approaches that are computationally tractable and can carefully capture
the underlying dynamical features.

We meet these challenges by studying the load data in an operator theoretic
framework where the data is understood by exploring and decomposing the
underlying dynamics manifold. In particular, working in the Koopman operator
setting allows us to represent the underlying dynamics in terms of its
intrinsic coordinates. These coordinates are the most predictable features and
efficiently capture the information of the load dynamics, and are approximated
from the spectral properties of the underlying Koopman operator. In the view
of the Koopman operator, the nonlinear dynamics becomes linear in the
$\infty$-dimensional space of observables. By studying the spectra of the
Koopman operator, we attain a spatiotemporal decomposition of nonlinear load
dynamics modes that determine a set of independent temporal and spatial
patterns in the evolution of the load dynamics. These modes, which unfold the
entanglement of the spatial and temporal energy demand information,
effectively encode the complex interactions in interconnected power grids, and
their richness allows a finite subset of them to approximately span the
evolving dynamics manifold. Hence, we achieve a minimum-order realization of
the large-scale power grids in terms of the intrinsic coordinates that provide
a globally linear description over the whole domain of the collected data. For
a better realization of Koopman modes for power stations that share similar
properties, we use a cluster-based Koopman forecasting for the load dynamics.
The proposed clustering method is based on a modern kernel approach, namely
PHATE [PHATE2019], which captures both local and global nonlinear data
structure using an information-geometric distance between data points. We
observe that the PHATE effectively reveals the stations with synchronized
dynamics.

The rest of the paper is organized as follows. Section 2 is the literature
review. Section 3 presents the Koopman operator, and Section 4 the data-driven
approach to the identification of the Koopman opeartor. In Section 5 we
describe the clustering approach. In Section 6 we presents our results.
Section 7 is devoted to the concluding remarks.

##  2 Literature review

#### Physics-informed approaches

Physics-informed approaches have been recently considered achieving load
forecasting models that are interpretable and generalizable
[Tartakovsky2019PhysInform; Stock2022PhysInform]. In load forecasting, where
large-scale interconnected subsystems work with complex interactions and under
various stochastic and exogenous factors, sparsity-promoting approaches, such
as SYNDy [Brunton2016SINDy], that preserve the original model’s physical
plausibility are particularly important to realize an effective basis for
modeling and achieve reduced dynamic representations [Stankovic2021SYNDy;
Cai2022SYNDy]. However, this basis often belongs to a pre-selected family, and
identifying it requires significant physics-derived information. In addition,
they lead to nonlinear models that block our path to the wealth of matured
linear approaches for subsequent analysis and design purposes.

#### Koopman

The significance of the Koopman operator is its connection with a different
class of physically-interpretable data-driven approaches that use a natural
basis for the learning process by extracting and decomposing the intrinsic
manifold of the data. Additionally, the decomposition gives a set of intrinsic
coordinates where the representation of the dynamics is linear. In the
geometrical study of dynamical systems, the intrinsic coordinates, which are
also known as coherent patterns of the dynamics, are equivalent to the
eigenfunctions of the Koopman operator [Mezic2005Spectral]. In the context of
signal processing, the Koopman operator may be viewed as an extension of the
Fourier analysis by replacing the predefined frequency basis in Fourier by
spectra that are intrinsic to the underlying dynamical system, as well as
replacing the Fourier temporal patterns by (intrinsic) spatiotemporal patterns
[Brunton2021Fourier]. Furthermore, since each Koopman eigenfunction recovers a
coherent temporal pattern with an associated timescale
[Giannakis2019Spectral], it is well suited for studying multiscale systems
[Giannakis2020Neuroscience; Giannakis2021Climate]. Although several settings
have recently been considered to present the Koopman approach, such as Dynamic
Mode Decomposition (DMD) [Tu2014DMD], Extended DMD (EDMD) [Williams2015EDMD],
graph DMD [Fjii2019GraphDMD], Hankel matrix analysis [Mezic2017Ergodic],
spectral moment estimation [Mezic2020Spectral], Galerkin
[Giannakis2019Galerkin], or deep learning [Lusch2018DL], we extract
spatiotemporal load dynamics patterns in an ergodic theory framework first
presented in [Mezic2005Spectral], and then extended based on kernel operators
and delay coordinates in [Giannakis2019Spectral; Giannakis2019Delay;
Klus2020RKHS; Giannakis2021RKHS]. This approach relies on the property of
integral kernel operators to approximate the Laplace-Beltrami operator
eigenfunctions [Giannakis2012NLSA], which recover the geometry of the
nonlinear manifold comprised by the data [Coifman2006Diffusion]. The Laplace-
Beltrami basis was particularly shown to address the limitation of the
traditional approaches such as Principal Component Analysis (PCA) for
capturing the intermittent features of signals [Giannakis2012NLSA;
Giannakis2013NLSA].

In the context of load dynamics, the DMD [Rowley2009; Schmid2010DMD], which is
the projection of Koopman operator into a linear subspace [Tu2014DMD], has
recently been applied for short-term electric load forecasting [Mohan2018DMD;
Kong2020DMD; Kutz2022GP], where it leads to better accuracy and stability
compared to different conventional approaches. Furthermore, the DMD is
adaptive to multiple exogenous factors including time, weather, seasonal
changes, economic aspects, and social activities, while such factors are not
used or estimated in the forecasting model [Mohan2018DMD; Kong2020DMD;
Kutz2022GP]. The DMD performance is also independent of rigorous training
phases, hand-crafted feature selection, or optimal parameter selection, and
the complexity of the method in terms of computation and speed is much lesser
than single structure and ensemble models for load forecasting [Mohan2018DMD].
However, since DMD inherits the limitations of singular value decomposition
and lacks sufficient delicacy to dissect rich nonlinear phenomena and the
associated transient dynamics [Schmid2022DMD; Brunton2021ModernKT], we use the
Koopman operator to employ nonlinear observables that are more adjustable with
the nonlinear and non-stationary load patterns.

#### Spatiotemporal forecasting

When matching large-scale distributed renewable energy with the urban energy
demand, understanding the spatiotemporal energy profiles will be important
[Mikkola2014spatio]. Specifically, while realizing high-resolution temporal
patterns is critical to predict future changes for demand-side management,
increasing the spatial resolution enables regional energy planning and
capturing the spatial dynamics of drivers of energy demand, and yet, since
regional and local energy grids are interconnected with national and
continental networks, high resolution is required everywhere to capture
multiscale effects [Eggimann2019spatio]. Spatiotemporal load forecasting has
recently been considered based on some LSTM deep learning methods
[Peng2021flexible; Zhu2022LSTM]. We depart from deep learning approaches and
present a physically-interpretable data-driven approach by identifying the
load dynamics in an operator theoretic framework, which, as discussed, is the
main objective of this work. Yet, there are two further features that add to
the advantages of this work: the computational efficiency and the
spatiotemporal patterns from Koopman eigenfunctions that are intrinsic to the
dynamics. The computational efficiency is due to the Galerkin approach
[Giannakis2019Galerkin] that we use to extract our forecasting model and it
does not require the challenges of deep learning approaches for rigorous
training phases and optimal parameter selections as also stated above.
Moreover, we present our model based on extracting multiscale spatiotemporal
modes where each temporal mode scales with a frequency of the load dynamics
and the corresponding spatial mode describes different stations’ energy
patterns at that scale. We observe that, by studying various spatial modes at
multiple scales we will be able to identify different energy consumption
behaviors, which get important insights for decision makers
[Niu2021understanding].

#### Clustering

To understand local demand characteristics on different urban scales, the
analysis of load data based on clustering of demand profiles has been proven
effective [Voulis2018understanding]. Furthermore, to increase the flexibility,
intuition, and interpretability of the data-driven approach, low-dimensional
embedding of different load profiles for coordinated visual exploration of
energy consumption patterns has recently been found useful
[Niu2021understanding; Peng2021flexible]. To this end, nonlinear embedding
techniques that project high dimensional data into a lower space and
simultaneously preserve the global and local data relative structures, such as
diffusion maps [Coifman2006Diffusion] and t-SNE [t-NSE2008], are superior to
conventional linear dimensionality reduction approaches such as the PCA
[Moon2000Book]. We utilize a method—known as PHATE- which is optimized for
two- or three-dimensional visualization, and was recently established in
[PHATE2019] based on potential of heat diffusion for affinity-based transition
embedding. We particularly use the PHATE to achieve a clustering approach
based on the intrinsic structure of the load data, and due to its improved
properties, compared to the suboptimal approaches such as t-SNE and PCA, for
addressing concerns of sensitivity to noise and scrambled global structure of
the data [PHATE2019]. Following the original work in [PHATE2019] that focused
on the biological data, in this paper we apply the PHATE to the extent of a
technological subsystem prototype, and explore its power for detecting
distinct energy patterns. Furthermore, we use a different variant of the
original approach in [PHATE2019] by using the variable bandwidth kernels
[Berry2016Variable] in the PHATE. Therefore, the bandwidth is heterogeneous
(variable) and adapted based on the local data structure, which improves the
geometric description of the data set by constructing the data graph using
more accurate realization of local data points similarities and interactions.

#### Spectral analysis in power grids

Frequency domain approaches are considered to describe the characteristics of
load profile in the time domain [Zhong2012frq] and methods such as Fourier
analysis [Zhong2015DFT], wavelet analysis [Li2017wavelet], or PCA
[Rydin2020PCA] are utilized for this purpose. However, the conventional
methods are required to impose conditions extrinsic to the dynamical system
and fail to capture frequencies in multiple scales, thus risking inaccurate
representations of the underlying system [Brunton2021Fourier]. In addition,
the construction by these conventional methods, and their advanced relatives
such as Dynamic Mode Decomposition (DMD) [Schmid2010DMD] is only accurate in
regions of the underlying state space where the dynamics are linear, and the
approaches are not tailored for the nonlinear data manifolds
[Giannakis2012NLSA]. In particular, they provide no information about the
chaotic parts of the dynamics as well as the continuous part of the system
spectrum, both of which are key to how power grids dynamics evolve
[Yu2003chaos; Halekotte2021chaos]. To meet these challenges, Koopman
eigenfunctions give rise to temporal patterns, each recovering a timescale
$2\pi/\omega$, that are connected with a Nonlinear Laplacian Spectral Analysis
(NLSA) [Giannakis2012NLSA] through eigenfunctions of the Laplace–Beltrami
operator[Coifman2006Diffusion] that recovers the geometry of the data manifold
[Giannakis2019Spectral]. Obtained temporal patterns are specifically highly
coherent with predictable temporal evolution intrinsic to the dynamical
system–markedly different from that of the eigenfunctions of covariance
operators approximated by conventional approaches.

#### Exogenous effects

While many approaches use the exogenous effects data to forecast the load,
such data are not available ahead of the forecasted time interval and
constructing additional forecasting models for exogenous effects seems
unavoidable. For example, in [Taylor2008weather], to improve the short-term
load forecasting based on approaches such as ARIMA and exponential smoothing,
an approach based on weather forecast was suggested. However, using a
secondary forecasting module again has its issues with forecasting exogenous
effects and often arise with significant error, which accordingly influence
the primary (load) forecasting task. In addition, environmental factors such
as environmental temperature have little impact on the electricity load of
industrial enterprises, which indicates it is difficult to select key factors
as input variables to build an electricity load forecasting model
[Hu2022industrial]. Furthermore, a secondary forecasting unit increases the
complexity of the whole forecast architecture, which limits the numerical
efficiency as well as the flexibility of the final model for interpretation.
Existing literature on load forecasting adopts a very casual approach in
considering the seasonality, based on either calendar month or some
meteorological parameter, which is inconsistent and inaccurate, especially
during transition periods, thus leading to high forecasting errors
[Sharma2022seasonal]. The Koopman operator extends to systems with input and
nonautonomous systems in a skew product dynamics setting
[Giannakis2020Neuroscience], where Koopman eigenfunctions encode different
dynamical features of system inputs. While Koopman eigenfunctions recover
various input features in the load forecasting application, as we show in this
paper, the eigenfunctions, and accordingly different encoded features, are
prioritized by computing the corresponding Dirichlet energy values.

##  3 Dynamical systems in operator-theoretic setting

###  3.1 Koopman operator

By the Koopman approach, the available data set is treated as observable of a
dynamical system with state space $\mathcal{X}$ and flow map
$\Phi^{t}:\mathcal{X}\mapsto\mathcal{X}$, $t\in\mathbb{R}$. Here,
$x_{t}=\Phi^{t}(x)$, $x\in\mathcal{X}$, gives the state of the system at time
$t$. In this setting, the recorded signal at $d$ different sensors is viewed
as an observation function $F:\mathcal{X}\mapsto\mathbb{R}^{d}$. The dynamical
system $(\mathcal{X},\Phi^{t})$ possesses ergodic measures that show there
exists a probability measure $\mu$ on $\mathcal{X}$, invariant under the flow
map $\Phi^{t}$, such that for every integrable function
$f:\mathcal{X}\mapsto\mathbb{C}$ the time average $\bar{f}$ of $f$ converges
to the expectation value $\bar{f}=\int_{\mathcal{X}}fd\mu$. Associated with
the triplet, $(\mathcal{X},\Phi^{t},\mu)$ we consider a Hilbert space
$\mathcal{H}=L^{2}(\mathcal{X},\mu)$ of square-integrable observables with
respect to $\mu$.

The group of unitary Koopman operators $U^{t}:\mathcal{H}\mapsto\mathcal{H}$
governs the evolution of observables under $\Phi^{t}$. That means given
$f\in\mathcal{H}$, $g=U^{t}f$ is defined as the observable satisfying
$g(x)=f(\Phi^{t}(x))$ for $x\in\mathcal{X}$. An observable
$\psi_{j}\in\mathcal{H}$ is a Koopman eigenfunction if it satisfies the
eigenvalue equation

$\begin{gathered}U^{t}\psi_{j}=e^{i\omega_{j}t}\psi_{j}\end{gathered}$ (1)

for all $t\in\mathbb{R}$. The variable $\omega_{j}$ (often known as
eigenfrequency) is a frequency with a real value associated with the
eigenfunction $\psi_{j}$. Therefore, in measure-preserving dynamical systems,
the Koopman eigenvalues remain on the unit circle in the complex plane, and
the corresponding eigenfunctions evolve periodically under the dynamics. This
is the key to the predictability of coherent patterns of dynamics. Koopman
eigenvalues and eigenfunctions occur as complex-conjugate pairs, and the
Koopman eigenfunctions corresponding to different eigenfrequencies are
orthogonal on $\mathcal{H}$. Consider the Koopman eigenfunctions space
$\mathcal{D}=\overline{\text{span}\\{\psi_{j}\\}}\subseteq\mathcal{H}$,
invariant under $U^{t}$. The space of observables lying in $\mathcal{D}$ has
an associated Koopman eigenfunction expansion: every $f\in\mathcal{D}$ admits
the decomposition $f=\sum_{j}\hat{f}_{j}\psi_{j}$ where
$\hat{f}_{j}=\left<f,\psi_{j}\right>$. Here, $\left<.,.\right>$ denotes the
inner product in $\mathcal{H}$. We also use $\parallel.\parallel$ to denote
the norm in this space. Moreover, the dynamical evolution of $f$ can be
computed in closed form via

$\begin{gathered}U^{t}f=\sum_{j}\hat{f}_{j}e^{i\omega_{j}}\psi_{j}\end{gathered}$
(2)

For every continuous flow $\Phi^{t}$, the family of operators $U^{t}$ has a
generator $V$, which is a skew-adjoint operator, defined as

$\begin{gathered}Vf:=\lim_{t\mapsto 0}\frac{1}{t}(U^{t}f-f),\ \ \ f\in
D(V)\subset L^{2}(\mathcal{X},\mu)\end{gathered}$ (3)

Operators $U^{t}$ and $V$ share the same eigenfunctions,

$\begin{gathered}V\psi_{j}=i\omega_{j}\psi_{j}\end{gathered}$ (4)

For ergodic systems, all eigenvalues of $V$ are simple.

###  3.2 Kernel operators and delay coordinates

Kernel methods have long tail in machine learning [Hofmann2008Kernel], and
especially, they have an established role in detecting nonlinear data
manifolds and order reduction [Coifman2006Diffusion]. Motivated by this
property to unfold complexity, dynamics-dependent kernels have recently been
studied in the context of dynamical systems [Berry2016Kernel]. In particular,
formulated in delay coordinate space, kernel operators have great capacity to
encode rich nonlinear phenomena [Giannakis2015Kernel].

Data can be a time-ordered set consisting of $N$ samples
$F(x_{1}),...,F(x_{N})$, each $F(x_{n})$ lying in $d$-dimensional data space.
Here, $x_{n}=\Phi^{n\Delta t}(x_{0})$ where $\Delta t$ is the sampling
interval. A delay coordinate map is constructed from $F$ by embedding
$\mathcal{X}$ in a manifold in $\mathbb{R}^{Qd}$ as

$\begin{gathered}F_{Q}(x)=(F(x),F(\Phi^{-\Delta t}x),...,F(\Phi^{-(Q-1)\Delta
t}(x)))\end{gathered}$ (5)

where $Q$ (an integer) is the number of delays. Then, define the kernel
function $k_{Q}:\mathcal{X}\times\mathcal{X}\mapsto\mathbb{R}_{+}$ which
provides a measure of similarity of points in $\mathcal{X}$, based on the
observation function $F_{Q}$. In this work, we use the radial Gaussian kernel,

$\begin{gathered}k_{Q}(x,x^{\prime})=\exp({-\frac{\parallel
F_{Q}(x)-F_{Q}(x^{\prime})\parallel^{2}}{\epsilon}})=\exp(-\frac{d_{Q}^{2}(x,x^{\prime})}{\epsilon})\end{gathered}$
(6)

where $\epsilon$ is the positive bandwidth parameter that can be varying in
$\mathcal{X}$ based on the available data density. Specifically, we use the
class of variable bandwidth kernels, also known as self-tuning kernels,
introduced in [Berry2016Variable]. Associated with the square-integrable
kernel $k_{Q}$ is a compact integral operator,

$\begin{gathered}K_{Q}f(x):=\int_{X}k_{Q}(x,y)f(y)d\mu(y)\end{gathered}$ (7)

Then a Markov kernel is constructed by normalizing $K_{Q}$ through functions
$\rho_{Q}=K_{Q}1_{X}$ and $\sigma_{Q}=K_{Q}(1/\rho_{Q})$ where $1_{S}$ is the
constant function equal to 1 on a set $S$:

$\begin{gathered}P_{Q}f:=\frac{\tilde{K}_{Q}f}{\tilde{K}_{Q}1_{X}}\end{gathered}$
(8)

where $\tilde{K}f=K_{Q}(\frac{f}{K_{Q}1_{X}})$.

In the limit, $Q\rightarrow\infty$ the kernel integral operator $P_{Q}$
commutes with $U^{t}$ [Giannakis2019Delay]. Since commuting operators have
common eigenspace, this allows us to compute eigenfunctions of $U^{t}$ through
expansions in the basis obtained from $P_{Q}$. Since $P_{Q}$ is compact, we
have a wide variety of data-driven tools to approximate integral operators.

###  3.3 Galerkin approximation of Koopman spectra

Solving for Koopman eigenvalues is a numerically ill-posed problem, since the
eigenvalues of its generator, $V$, are dense on the imaginary axis in ergodic
measure-preserving systems. This property adversely affects the stability of
numerical approximations of Koopman eigenvalues and eigenfunctions. To address
this challenge, finding the eigenvalues of $V$ is replaced by approximating
them through a regularized operator $L_{\theta}$ with a small amount $\theta$
of judiciously constructed diffusion added for regularization
[Giannakis2019Spectral; Giannakis2019Delay]. Eigenfunctions of the Koopman
operator are sought in a Sobolev subspace $\mathcal{H}^{2}\subset\mathcal{H}$
[Giannakis2019Delay], on which $V$ is a bounded operator. Based on the
standard Galerkin approach, the regularized Koopman eigenvalue problem is to
find $\gamma\in\mathbb{C}$ and $z\in\mathcal{H}^{2}$, such that for every
$f\in\mathcal{H}$ the following weak formulation holds:

$\begin{gathered}\left<z,L_{\theta}f\right>=\gamma\left<f,z\right>,\ \ \
L_{\theta}=V-\theta\Delta\end{gathered}$ (9)

where the operator $\Delta$ is defined based on the eigenfunctions of the
Laplace-Beltrami operator [Giannakis2019Spectral; Giannakis2019Delay]. In (9),
$\gamma$ and $z$ are weak eigenvalues and eigenfunctions of $L_{\theta}$, and
$f$ is test function. The role of the diffusion term is to penalize the
eigenfunctions of $V$ with large values of a Dirichlet energy functional,
$E(f)=\frac{\left<f,\Delta f\right>}{\parallel f\parallel^{2}}$, providing a
measure of roughness of functions in $\mathcal{H}$. In this way, the dynamics
is projected into the least rough set of Koopman eigenfunctions. The
eigenfunctions $\varphi_{j}$, of the kernel operator $P_{Q}$ defined in
Section 3.2 are used as a basis to solve the Galerkin problem (9), as well as
to approximate the eigenfunctions of the Laplace-Beltrami operator. While the
regularized operator $L_{\theta}$ and the generator $V$ share the same
eigenfunctions, the eigenvalues of $L_{\theta}$ are parameterized as
$\gamma_{\theta}=i\omega-\theta\eta$ where $i\omega$ is an eigenvalue of $V$
and $\eta$ an eigenvalue of $\Delta$.

##  4 The data-driven approach

###  4.1 Data-driven Koopman eigenfunctions on sampled data

To put the above approach in a numerical setting, we order the load data as
$x_{k}=x(k\tau)$, $k=1:N$, with the sampling interval $\tau$. Here, each
$x_{k}\in\mathbb{R}^{d}$ is a vector of load data at $d$ spatial points
measured in time $t=k\tau$. The use of the Markov kernel operator is to
specify an orthonormal basis for the Koopman eigenfunctions based on the data
sampled on the orbit. The action of the kernel operator $K_{Q}$ is
approximated through a matrix $K$ with the elements given as

$\begin{gathered}K_{ij}=k_{Q}(x_{i},x_{j})=\exp(-\frac{d^{2}_{Q}(x_{i},x_{j})}{\epsilon}),\\\
d^{2}_{Q}(x_{i},x_{j})=\frac{1}{Q}\sum_{k=0}^{Q-1}\parallel
x_{i-k}-x_{j-k}\parallel^{2}\end{gathered}$ (10)

The normalized Markov kernel $P_{Q}$ is approximated through the elements of
the matrix $P$ as

$\begin{gathered}P_{ij}=\frac{K_{ij}}{(\sum_{k=1}^{N}K_{ik}q_{k}^{-1/2})q_{j}^{1/2}},\
\ \ q_{i}=\sum_{k=1}^{N}K_{ik}\end{gathered}$ (11)

The approximated eigenvalues, $\lambda_{j}$, and eigenfunctions,
$\varphi_{j}$, of the Markov kernel operator are obtained by solving the
algebraic eigenvalue problem for $P$,

$\begin{gathered}P\varphi_{j}=\lambda_{j}\varphi_{j}\end{gathered}$ (12)

The matrix $P$ approximates the operation of $P_{Q}$ as an ergodic average
that is the base idea in kernel analog forecasting [Giannakis2020Operator].
From a numerical perspective, for large sample numbers $N$, $P$ may be made
sparse by zeroing all but the largest $k_{nn}\ll N$ elements in each row of
$K$, and symmetrizing the resulting sparse matrix.

As a Markov matrix, $P$ has real eigenvalues ordered as
$1=\lambda_{1}>\lambda_{2}\geq\lambda_{3}\geq...$, as well as real
eigenvectors $\varphi_{j}$ that are mutually orthogonal in $\mathbb{R}^{N}$
(we observe that the first eigenvector corresponding to $\lambda_{1}=1$ is the
constant eigenvector $\varphi_{1}=(1,...,1)^{T}\in\mathbb{R}^{N}$ that gets
assigned zero Dirichlet energy). As a key feature, since $P_{Q}$ commutes (and
hence has common invariant subspaces) with the Koopman operator $U^{t}$ as
$Q\rightarrow\infty$, Koopman eigenfunctions can be approximated by a linear
combination of kernel operator eigenfunctions,

$\begin{gathered}\psi_{i}=\sum_{j=1}^{l}c_{ji}\varphi_{j},\ \ \
i=1,...,l\end{gathered}$ (13)

where $\psi_{i}\in\mathbb{R}^{N}$ denotes the $i$-th eigenfunction of $U^{t}$
sampled at $N$ data points. Hence, $l\leq N$ first eigenfunctions of the
kernel operator are used to approximate the Koopman eigenfunctions

The matrix $c_{l\times l}=[c_{ij}]$ is determined by a data-driven Galerkin
approach [Giannakis2019Galerkin], where the weak formulation (9) is
approximated by the following matrix generalized eigenvalue problem:

$\begin{gathered}A\hat{c}=\gamma B\hat{c}\end{gathered}$ (14)

where

$\begin{gathered}A_{ij}=\left<\varphi_{i},V\varphi_{j}^{(2)}\right>-\theta\left<\varphi_{i},\Delta\varphi_{j}^{(2)}\right>\\\
B_{ij}=\left<\varphi_{i},\varphi_{j}^{(2)}\right>\\\
\varphi_{i}^{(2)}=\frac{\varphi_{i}}{\eta_{i}}\ \ \
\eta_{i}=\frac{1}{\epsilon}(\lambda_{i}^{-1}-1)\end{gathered}$ (15)

for $i,j=1,...,l$. Each column of matrix $c$ corresponds to a solution of
$\hat{c}$ in (14). Solutions $(\gamma_{i},\psi_{i})$ are placed in order of
increasing Dirichlet energy, calculated for each eigenfunction $\psi_{i}$ as

$\begin{gathered}E(\psi_{i})=\frac{\left<\psi_{i},\Delta\psi_{i}\right>}{\parallel\psi_{i}\parallel^{2}}\end{gathered}$
(16)

Intuitively, the Dirichlet energy measures the spatial variability (or
roughness) of functions on $\mathcal{X}$. A highly oscillatory eigenfunction
in $\mathcal{X}$ has a large Dirichlet energy, while $E(\psi_{i})$ will be
small if $\psi_{i}$ is spatially smooth.

Therefore, while many existing approaches require a library of functions to
construct the Koopman invariant subspaces, using the kernel and delay
coordinates approaches, a data-driven basis is obtained for the Koopman
eigenfunctions. This basis is dynamics-dependent and optimal since it
approximates the Laplace-Beltrami operator eigenfunctions
[Giannakis2019Spectral; Giannakis2019Delay], which recover the geometry of the
nonlinear manifold on which data are collected [Coifman2006Diffusion]. The
Koopman eigenfunctions corresponding to nontrivial eigenvalues ($\lambda=0$
for continuous systems and $\lambda=1$ for discrete systems) form complex-
conjugate pairs with equal Dirichlet energies. The nontrivial eigenvalue has a
constant eigenfunction over $\mathcal{X}$ with zero Dirichlet energy.

Algorithm 1 presents a pseudocode for the approximation of the Koopman
eigenfunctions. We use the approach established in [Berry2016Variable] to
achieve a variable bandwidth kernel where the bandwidth parameter is adapted
inversely proportional to the sampling density:

$\begin{gathered}\epsilon=\epsilon(x_{i},x_{j})\propto\frac{1}{\rho^{\alpha}(x_{i})\rho^{\alpha}(x_{j})},\
\ \ \alpha\in\mathbb{R}^{+}\end{gathered}$ (17)

where the sampling density function $\rho(.)$ and the constant $\alpha$, along
with an estimation of the manifold dimension $\hat{m}$, are computed based on
solving an optimization problem for the available data. The adaptive mechanism
increases the bandwidth whatever the sampling density is small, and it
decreases the bandwidth at regions where data is dense. Several approaches may
be used to estimate the data density $\rho(.)$ [Giannakis2019Spectral], …

Algorithm 1 Data-driven Koopman eigenfunctions approximation.

  * •

Input

    * –

Time series $\\{x_{k}\\}_{k=1}^{k=N}$, $x_{k}\in\mathbb{R}^{d}$, of the load
sampled at $d$ power stations every $\tau$ time units

    * –

Number of delays $Q$

    * –

Number of retained nearest neighbors $k_{nn}$

    * –

Number of Markov kernel eigenfunctions $l\leq N$ used in Galerkin
approximation

    * –

Koopman regularization parameter $\theta\geq 0$

    * –

Number of Koopman eigenfunctions $l^{\prime}\leq l$ to be computed

  * •

Output

    * –

Koopman eigenvalues $\gamma_{1},...,\gamma_{l^{\prime}}$ with
$\gamma_{k}\in\mathbb{C}$

    * –

Koopman eigenfrequencies $\omega_{1},...,\omega_{l^{\prime}}$ with
$\omega_{k}\in\mathbb{R}$

    * –

Koopman eigenfunctions, $\psi_{1},...,\psi_{l^{\prime}}$ with
$\psi_{k}\in\mathbb{C}^{N}$

  * •

Steps

1:Compute the $N\times N$ pairwise distances $d_{Q}(x_{i},x_{j})$ using (10)

2:For each $i$ only retain the $k_{nn}$ smallest values of
$d_{Q}(x_{i},x_{j})$ and consider the set
$\mathcal{N}_{k_{nn}}(x_{i})=\\{x_{j}|x_{j}\ \text{in the}\ k_{nn}\
\text{neighbourhood of}\ \ x_{i}\\}$

3:Symmetrize $d_{Q}(x_{i},x_{j})$ by augmenting the list of retained distances
for given $i$ if $x_{i}\in\mathcal{N}_{k_{nn}}(x_{j})$ but
$x_{j}\notin\mathcal{N}_{k_{nn}}(x_{i})$

4:Determine the optimal bandwidth function $\epsilon=\epsilon(x_{i},x_{j})$
using the procedure described in [Berry2016Variable]

5:Form the $N\times N$ matrix $K$ in (10), treating all pairwise kernel values
$K_{ij}$ not associated with the retained distance data as zero.

6:Compute the $N\times N$ sparse matrix $P$ in (11)

7:Compute the $l$ largest eigenvalues $\lambda_{k}$ of $P$ and the
corresponding kernel eigenfunctions $\psi_{k}$

8:Form the $l\times l$ Galerkin approximation matrices $A$ and $B$ in (15)

9:Solve the generalized eigenvalue problem (14) and obtain $c_{k}$ as the
$k$-th column of the matrix $c_{l\times l}$ and the corresponding $\gamma_{k}$
for $k=1,...,l$

10:Compute Koopman eigenfunctions $\psi_{k}$ for $k=1,...,l^{\prime}$ using
(13)

11:Compute the Dirichlet energy $E(\psi_{k})$ of each eigenfunction using (16)

12:Order the pairs $(\gamma_{k},\psi_{k})$ in order of increasing Dirichlet
energies

13:Compute eigenfrequencies $\omega_{k}=\text{Im}(\gamma_{k})$

###  4.2 Spatiotemporal load reconstruction

In this paper, our purpose in reconstructing the time series of the data is to
realize how many modes actually contribute to the evolution of the system.
This also allows for a spatiotemporal decomposition of the data and extracting
the chronos space (temporal modes) and topos space (spatial patterns) of the
dynamics. Each Koopman eigenfunction induces a time series
$\psi(x_{1}),\psi(x_{2}),...$ that are the temporal modes. The spatial
patterns are constructed in an NLSA space [Giannakis2012NLSA], which has a
higher dynamical significance than the standard PCA …

First, the projections $A_{k}(t)$ of the observable $F$ onto the Koopman
eigenfunctions $\psi_{k}$ are considered as
$A_{k}(t)=\left<\psi_{k},F\circ\Phi^{t}\right>$, which is numerically
evaluated as [Giannakis2018Turbulent]

$\begin{gathered}A_{k}(q\tau)=\frac{1}{N}\sum_{n=1}^{N}\psi^{*}_{nk}x_{n-q},\
\ \end{gathered}$ (18)

where the Koopman eigenfunctions are considered as
$\psi_{k}=(\psi_{1k},...,\psi_{Nk})^{T}$. The component $A_{k}(0)$ is called
the $k$-th Koopman mode associated with the observable $F$
[Mezic2005Spectral]. We use the elements of $A_{k}$ as the spatial patterns
generated by the Koopman eigenfunctions. The corresponding spatiotemporal
patterns associated with the pair $(\psi_{k},A_{k})$ are given by

$\begin{gathered}x_{n}^{(k)}=\frac{1}{Q^{\prime}}\sum_{q=0}^{Q^{\prime}-1}A_{k}(q\tau)\psi_{(n+q,k)},\
\ \ Q^{\prime}=\text{min}\\{Q,N-1-n\\},\ \ 1\leq n\leq N\end{gathered}$ (19)

where $x_{n}^{(k)}$ is the contribution of the $k$ -th Koopman eigenfunction
to the reconstruction of $x_{n}=x(n\tau)$. Since Koopman eigenfunctions form
complex-conjugate pairs $(\psi_{k},\psi_{k+1})$, the sum
$x_{n}^{(k)}+x_{n}^{(k+1)}$ is required to obtain a real-valued load
reconstruction by each pair. By (19), a reconstruction at time $n\tau$ is a
superposition of $Q$ snapshots. This is a result of the time-lagged embedding
and different from a standard PCA analysis [Giannakis2018Turbulent].

###  4.3 Nystrom extension of Koopman eigenfunctions

By the original approach developed in [Giannakis2019Spectral;
Giannakis2019Delay], we obtain the Koopman eigenfunctions sampled at the data
points $X_{s}=\\{x_{1},...,x_{N}\\}$, as described in Section 4.1. To evaluate
the Koopman eigenfunctions at out-of-sample points, in this work, we use the
Nystrom extension approach to extend the Koopman eigenfunctions to
$x\in\mathcal{X}/X_{s}$. For this purpose, we first apply the Nystrom approach
to extend the eigenfunctions $\varphi_{k}$ of the Markov kernel operator and
then use (13) for out-of-sample evaluation of Koopman eigenfunctions. We
present our approach based on the method established in [Coifman2006Nystrom].

Assume the Markov integral operator $P_{Q}$ described in Section 3.2 has the
kernel $p_{Q}:\mathcal{X}\times\mathcal{X}\mapsto\mathbb{R}_{+}$. The Nystrom
extension method is used to solve the eigenfunction of integral equations in
the form

$\begin{gathered}\int
p_{Q}(x,y)\varphi(y)\nu(y)dy=\lambda\varphi(x)\end{gathered}$ (20)

where $\nu$ indicates the underlying probability density function,
$\varphi(x)$ represents the eigenfunction, and $p_{Q}(x,y)$ denotes the
similarity between $x$ and $y$. We consider the data
$X_{s}=\\{x_{1},x_{2},...,x_{N}\\}$, sampled from the manifold $\mathcal{X}$,
as $N$ landmark points at which the eigenfunction $\varphi$ is sampled. For
any given point $x\in\mathcal{X}$ using the Nystrom method, we have

$\begin{gathered}\frac{1}{N}\sum_{i=1}^{N}p_{Q}(x,x_{i})\hat{\varphi}(x_{i})=\lambda\hat{\varphi}(x)\end{gathered}$
(21)

where $\hat{\varphi}(x)$ is an approximation to the true $\varphi(x)$.
Equation (21) cannot be solved directly, as $\hat{\varphi}(x)$ and $\lambda$
are both unknown. Denote the similarity matrix between the landmark points by
$P$ with $P_{ij}=p_{Q}(x_{i},x_{j})$. Substituting $x$ with $x_{i}$ in
Equation (21) and writing it in matrix form, we have

$\begin{gathered}P\hat{\varphi}_{k}=N\hat{\lambda}_{k}\hat{\varphi}_{k},\ \ \
k=1,...,N\end{gathered}$ (22)

For an unsampled point $x$, the $i$-th eigenfunction at $x$ can be
approximated as

$\begin{gathered}\hat{\varphi}_{i}(x)=\frac{1}{N\hat{\lambda}_{i}}\sum_{j=1}^{N}p_{Q}(x,x_{j})\hat{\varphi}_{i}(x_{j})\end{gathered}$
(23)

The Koopman eigenfunctions are correspondingly evaluated by (13).

###  4.4 Forecasting by Koopman eigenfunctions

While the original work in [Giannakis2019Spectral; Giannakis2019Delay]
suggests (2) for the prediction of observables in the Koopman invariant
subspaces, there are two important problems that should be addressed in real
world applications. First, while the regular part of the dynamics, e.g., the
periodic and quasiperiodic dynamics, is connected with the discrete spectra
and has its corresponding Koopman invariant subspace that completely describes
its evolution, the irregular dynamics connected with the continuous spectra,
e.g. the chaotic and mixing dynamics, have no corresponding Koopman invariant
subspaces [Mezic2005Spectral]. As such, some eigenfunctions recovered by the
kernel and delay coordinates approaches are not exact Koopman eigenfunctions,
and in fact, they are only $\epsilon$-approximate eigenfunctions that satisfy
(1) with a bound $\epsilon$ controlled by the length of the delay-embedding
window [Giannakis2021Coherent]. In particular, if $z$ is an
$\epsilon$-approximate eigenfunction, then there is an $\omega\in\mathbb{R}$,
called $\epsilon$-approximate eigenfrequency, such that

$\begin{gathered}\parallel U^{t}z-e^{i\omega t}z\parallel\leq\epsilon\parallel
z\parallel\end{gathered}$ (24)

where $\epsilon$ can be made arbitrarily small as the embedding window $Q$
increases. Here, to maintain the original concept in [Giannakis2021Coherent],
we do not change the “$\epsilon$” notation, but we do not confuse it with the
kernel bandwidth parameter.

Therefore, while the regular dynamics can be exactly predicted using the
Koopman eigenfunctions, the irregular dynamics is approximated by error. To
reduce this error, in this work, we suggest another modification by assuming
that some eigenfunctions $\psi_{i}(x)$ may not be completely single-frequency
evolving functions as in (2), and their error with respect to a regular
eigenfunction can be (partly) compensated by being coupled with the other
eigenfunctions. To this end, we consider the following linear evolution
through an $l^{\prime}\times l^{\prime}$ matrix $K_{\psi}$:

$\begin{gathered}\psi(x(t+1))=K_{\psi}\psi(x(t))\end{gathered}$ (25)

where $\psi=(\psi_{1},...,\psi_{l^{\prime}})^{T}$ is a column vector including
Koopman eigenfunctions and the matrix $K_{\psi}$ is computed based on data as
follows. First, we consider the following matrices involving evaluations of
eigenfunctions on sampled data using Algorithm 1:
$\Psi^{-}=[\psi(x_{1}),...,\psi(x_{N-1})]$ and
$\Psi^{+}=[\psi(x_{2}),...,\psi(x_{N})]$. Then, $K_{\psi}$ is the solution of
the following optimization:

$\begin{gathered}\min_{K_{\psi}}\parallel\Psi^{+}-K_{\psi}\Psi^{-}\parallel_{2}\end{gathered}$
(26)

which is solved using normal equations. In the ideal case where the underlying
dynamics is regular, e.g., periodic or quasiperiodic dynamics, the
eigenfunctions recovered by the kernel operator and delay coordinates
approaches are single-frequency and $K_{\psi}$ will be block-diagonal (recall
that Koopman eigenfunctions are complex conjugate pairs, so that they evolve
in pair and each block will be a $2\times 2$ skew-symmetric matrix). However,
in the case of irregular dynamics, e.g., chaotic or mixing dynamics, some
eigenfunctions may evolve through coupling with others and the block-diagonal
structure will (slightly) break. To illustrate this discussion, Figure 1 shows
how $k_{\psi}$ gets closer and closer to the ideal block-diagonal matrix by
widening the embedding window. In particular, for, $Q\geq 1000$ the proposed
approach recovers coherent patterns that fall into line with exact Koopman
eigenfunctions (with $\epsilon$ values close to zero in (24)).

(a)  (b)  (c)  (d)

Figure 1: Visualization of $K_{\psi}$ for different embedding windows
indicates that, increasing $Q$, the utilized approach decomposes the dynamics
into coherent patterns with each pattern having decoupled dynamics. (a)
$Q=100$, (b) $Q=500$, (c) $Q=1000$, (d) $Q=10000$.

The second problem to be addressed before any real world application is the
data-driven Koopman mode decomposition $\hat{f}_{j}=\left<f,\psi_{j}\right>$
needed to evaluate (2). In fact, due to $\epsilon$-approximate eigenfunctions,
the orthogonality of the recovered eigenfunctions, which is the necessary
condition to derive (2), may not be satisfied for sufficient accuracy.
Furthermore, an accurate data-driven evaluation of the inner product
$\left<f,\psi_{j}\right>$ requires enough data sampled over $\mathcal{X}$, and
any shortage in this case will be a source of error. To reduce these types of
error, in this work, we consider a linear decoder which maps the Koopman
eigenfunctions to system states. By linear decoder, we assume that the
considered set of Koopman eigenfunctions is large enough to span the system
state: $x\in\text{span}\\{\psi_{k}(x)\\}_{k=1}^{l^{\prime}}$. Therefore, the
system state is expanded in terms of Koopman eigenfunctions through a $d\times
l^{\prime}$ matrix $C$ as follows:

$\begin{gathered}x=C\psi(x)\end{gathered}$ (27)

where $x\in\mathbb{R}^{d}$ represents the system state vector. To determine
$C$ based on available data, we consider another optimization problem as
follows:

$\begin{gathered}\min_{C}\parallel X-C\Psi\parallel_{2}\end{gathered}$ (28)

where $X$ and $\Psi$ are the time-ordered matrices $X=[x_{1},...,x_{N}]$ and
$\Psi=[\psi(x_{1}),...,\psi(x_{N})]$.

##  5 Clustering the data

We use a clustering-based Koopman approach to cluster the stations according
to their temporal power patterns, and then forecast the future energy for each
pattern of stations. To this end, we apply a time series clustering analysis
for electricity pattern detection. To some extent, power generation usually
follows a steady pattern of evolution over time. These patterns can be
detected by mapping the power temporal data to the lower dimensional space.
There are many dimensionality-reduction methods for visualization, of which
the most commonly used are PCA [Moon2000Book], diffusion maps
[Coifman2006Diffusion], and t-SNE [t-NSE2008]. To address some concerns
regarding sensitivity to noise, scrambled global structure of data,
computational scalability, optimizing for two- or three-dimensional
visualization, and intrinsic structure of the data, a method based on
potential of heat diffusion for affinity-based transition embedding—known as
PHATE- was recently established [PHATE2019]. While the original work in
[PHATE2019] was mainly concerned with the biological data, in this paper we
apply the PHATE to the extent of a physical subsystem prototype; that is, for
detecting distinct power generation patterns. We first consider the
normalization of a typical vector $x$ as

$\begin{gathered}\hat{x}=\frac{x-x_{min}}{x_{max}-x_{min}}\end{gathered}$ (29)

where $\hat{x}$ denotes the normalized vector. In (29), $x$ is a time-ordered
vector with the minimum and maximum values $x_{min}$ and $x_{max}$. For
clustering, we consider the load data at $d$ spatial points as time-ordered
vectors in $\mathbb{R}^{N}$ and denote the load recorded at station $i$ as
$y_{i}=[y_{i}(t_{1}),...,y_{i}(t_{N})]^{T}$, $i=1,...,d$. Empirical
similarities between different load patterns are measured for the normalized
vectors $\hat{y}_{i}$ and using the variable bandwidth kernel explained in
Section 4.1. This leads to a kernel matrix, $\hat{P}_{d\times d}$ with the
elements $\hat{P}_{ij}$ indicating the similarity between each pair
$(\hat{y}_{i},\hat{y}_{j})$. We note that it does not make sense to consider
any time-delay $Q$ when evaluating similarities of spatial points. In general,
$\hat{P}$ is quite different from the matrix $P$ used in Section 4.1 to
compute similarities between different time points (look at the difference
between the vectors $y_{i}\in\mathbb{R}^{N}$ and $x_{i}\in\mathbb{R}^{d}$
considered in Section 4.1). After defining an initial diffusion operator like
$\hat{P}$, the PHATE takes some deliberately designed steps. First, the
diffusion operator is raised to its $t^{\prime}$-th power, where the diffusion
time scale $t^{\prime}$ is chosen optimally based on von Neumann entropy. This
provides a robust intrinsic data distance that preserves the global structure
of the data, and additionally low-pass filters the data. To resolve
instabilities and embed the global structure in low (two or three) dimensions,
in the next step, PHATE considers
$\hat{U}^{t^{\prime}}_{i}=-\text{log}(\hat{P}^{t^{\prime}}_{i})$, where
$\hat{P}^{t^{\prime}}_{i}$ is the $i$-th row of $\hat{P}^{t^{\prime}}$, and
computes the potential distances:

$\begin{gathered}\Gamma^{t^{\prime}}_{ij}=\parallel\hat{U}^{t^{\prime}}_{i}-\hat{U}^{t^{\prime}}_{j}\parallel_{2}\end{gathered}$
(30)

The last step involves a metric multidimensional scaling (metric MDS) to embed
the data into a low-dimensional (in our work 3-dimensional) space. This is
done by minimizing the following ‘stress’ function:

$\begin{gathered}\text{stress}(\hat{z}_{1},...,\hat{z}_{d})=\sqrt{\sum_{i,j}(\Gamma^{t^{\prime}}_{ij}-\parallel\hat{z}_{i}-\hat{z}_{j}\parallel)^{2}/\sum_{i,j}(\Gamma^{t^{\prime}}_{ij})^{2}}\end{gathered}$
(31)

over embedded $m$-dimensional coordinates $\hat{z}_{i}\in\mathbb{R}^{m}$ of
data points. While the classical MDS assumes the input distances directly
correspond to low-dimensional Euclidean distances, which is overly restrictive
in the PHATE setting, the metric MDS relaxes this assumption by only requiring
the input distances to be a distance metric. The considered clustering
approach is summarized in Algorithm 2.

Algorithm 2 Low-dimensional embedding of data using PHATE for clustering

  * •

Input

    * –

Time series of each station as $y_{i}=[y_{i}(t_{1}),...,y_{i}(t_{N})]^{T}$,
$y_{i}\in\mathbb{R}^{N}$, $i=1,...,d$

    * –

Number of retained nearest neighbors $k_{nn}$

  * •

Output

    * –

PHATE coordinates PHATE${}_{j}\in\mathbb{R}^{d}$

  * •

Steps

1:Get the normalized data $\hat{y}_{i}$ according to (29)

2:Compute the $d\times d$ pairwise distances $d_{1}(\hat{y}_{i},\hat{y}_{j})$
using (10)

3:For each $i$ only retain the $k_{nn}$ smallest values of
$d_{1}(\hat{y}_{i},\hat{y}_{j})$ and consider the set
$\mathcal{N}_{k_{nn}}(\hat{y}_{i})=\\{\hat{y}_{j}|\hat{y}_{j}\ \text{in the}\
k_{nn}\ \text{neighbourhood of}\ \ \hat{y}_{i}\\}$

4:Symmetrize $d_{1}(\hat{y}_{i},\hat{y}_{j})$ by augmenting the list of
retained distances for given $i$ if
$\hat{y}_{i}\in\mathcal{N}_{k_{nn}}(\hat{y}_{j})$ but
$\hat{y}_{j}\notin\mathcal{N}_{k_{nn}}(\hat{y}_{i})$

5:Determine the optimal bandwidth function
$\epsilon=\epsilon(\hat{y}_{i},\hat{y}_{j})$ using the procedure described in
[Berry2016Variable]

6:For $d_{1}(\hat{y}_{i},\hat{y}_{j})$, form the $d\times d$ matrix $\hat{K}$
in (10), treating all pairwise kernel values $\hat{K}_{ij}$ not associated
with the retained distance data as zero.

7:For $\hat{K}$, compute the $d\times d$ sparse matrix $\hat{P}$ in (11)

8:Compute an optimal value for the diffusion time scale $t^{\prime}$ using the
algorithm given in [PHATE2019]

9:For each row $\hat{P}_{i}^{t^{\prime}}$ of $\hat{P}^{t^{\prime}}$ compute
$\hat{U}^{t^{\prime}}_{i}=-\text{log}(\hat{P}^{t^{\prime}}_{i})$

10:Compute the potential distances $\Gamma^{t^{\prime}}_{ij}$ in (30)

11:Do a metric MDS by minimizing the ‘stress’ function (31) and obtain
$\hat{z}_{i}\in\mathbb{R}^{m}$, $i=1,...,d$.

12:Set PHATE${}_{j}=[\hat{z}_{1}(j),...,\hat{z}_{d}(j)]^{T}$, $j=1,...,m$,
where $\hat{z}_{i}(j)$ denotes the $j$-th element of vector $\hat{z}_{i}$.

##  6 Results

In this section, we present the results of applying the Koopman operator
approach to the data described in Section 6.1. A schematic of the proposed
numerical approach to load forecasting using Koopman eigenfunctions is shown
in Figure 2. We evaluate the performance of our proposed approach according to
the following RMSE (root mean squared error),

$\begin{gathered}\text{RMSE}=\sqrt{\frac{1}{n_{x}}\sum_{i=1}^{n_{x}}(x_{i}-x^{f}_{i})^{2}}\end{gathered}$
(32)

which applies to a vector $x=[x_{1}\ ...\ x_{n_{x}}]^{T}$ and its forecasted
value $x^{f}=[x^{f}_{1}\ ...\ x^{f}_{n_{x}}]^{T}$. All the error and RMSE
values in this paper were computed for the normalized load vectors, so that
the RMSE is limited (scaled) between 0 and 1.

(a)

(b)

Figure 2: (a) Computing Koopman eigenfunctions. (b) Load forecasting using
Koopman.

###  6.1 Data description

The data used in this paper are a large-scale dedicated data set for a
renewable electric power system [Jensen2017Data] that models the continental
European electricity system, including demand and renewable energy inflows for
the period 2012-2014. The resolution is as high as 50 km and 1 hour, and the
extent is as large as that of mainland Europe and 3 years. The dataset
combines a transmission network model, as well as information for generation
and demand. It is collected from the transmission network of 1494 transmission
buses with 2156 lines. Generation includes conventional generators with their
technical and economic characteristics, as well as weather-driven forecasts
and corresponding realizations for renewable energy generation. An overview of
the dataset and the power grid is seen in Figure 3.

(a)  (b)  (c)

Figure 3: Overview of dataset. (a) Distribution of generators. (b)
Transmission grid and the 3-year average power of each generator in MWh/h. (c)
Average load of all power plants.

###  6.2 Reduced-order modeling of the load dynamics

Many high-dimensional systems of interest evolve on a low-dimensional manifold
or attractor that is well-approximated using a low-rank basis
[Brunton2016SINDy]. We show that the Koopman eigenfunctions yield a natural
such basis for the load dynamics. Specifically, the load of all stations
(1494) during January 2014 are reconstructed using $l=100$ Koopman
eigenfunctions, with Figure 4 showing some example load reconstructions. Here,
we apply Algorithm 1 to the load data available for January 2014 with delay
$Q=500$ and the regularizing parameter $\theta=10^{-9}$. After identifying the
Koopman eigenfunctions, the load is reconstructed by (18). The load of other
months can be reconstructed in the same manner, and Figure 4 serves as an
example. Moreover, we can reconstruct the load over larger time horizons
(e.g., for one year), and for that, we need to increase the length of the
delay. However, this will induce a huge computational burden due to the data
that is substantially large in both temporal and spatial dimensions.

###  6.3 Nonlinear spatiotemporal spectral analysis of the load dynamics

Therefore, Koopman eigenfunctions effectively encode the spatiotemporal load
dynamics and represent the dynamics on low-order manifold. The dynamics modes
are constructed in a nonlinear Laplacian spectral analysis (NLSA) space
[Giannakis2012NLSA; Giannakis2013NLSA], which has a higher dynamical
significance than the space obtained by the standard PCA or POD. This leads to
coherent patterns of the dynamics that identify the most predictable and
robust features in evolving dynamics [Giannakis2018Turbulent;
Giannakis2021Climate; Giannakis2021Coherent]. Figure 5 shows these patterns
for the load dynamics. These patters appear with the least Dirichlet energy
values in Koopman eigenfunctions. Similar modes are detected for other months,
and Figure 5 represents January 2014 as an example.

By Figure 5, the most variance of the load dynamics is captured by the
frequencies that scale with zero, a month, a week, a day, and a half day.
Thus, the Koopman framework offers an interpretable machine learning approach
that allows us to detect the composing elements of the dynamics. Specifically,
in the load dynamics, while the zero-frequency mode is corresponding to the
average loads generated in individual power stations (compare the spatial
pattern in the first row of Figure 5 to average loads in Figure 3), the mode
scaling with a month in the second row of Figure 5 is due to the seasonal
conditions (temperature, solar radiation, …) that may change during a month.
The modes in the third row of Figure 5, corresponding to a week frequency, is
owing to the working days and holidays during a week, and the mode scaling
with a day frequency in the fourth row of Figure 5 is due to changes of solar
radiation during 24 hours. The last row in Figure 5 detects a half day
frequency that is made by the working hours during a day and generally the
different power usage patterns during date and night.

Thus, with the Koopman operator, we identify the multiscale dynamics in power
grids and detect the contribution of modes at different scales in the load
dynamics. While the general methods of Fourier analysis and POD have been
conventionally utilized to decompose signals into their frequency contents
[Holmes1996], the spectra of the Koopman operator lead to the decomposition of
the dynamics in an optimal nonlinear basis [Brunton2021Fourier], i.e., in the
space of Laplace-Beltrami operator eigenfunctions [Giannakis2012NLSA] which
recover the geometry of the nonlinear data manifold data
[Coifman2006Diffusion]. While the conventional approaches are mainly proposed
to extract the temporal patterns, the optimal basis discovered from Koopman
spectral analysis determines spatial patterns corresponding to each temporal
mode (see the right column in Figure 5). In addition, we do not need the prior
knowledge of a base frequency for the Koopman modal analysis.

The temporal patterns in the middle column of Figure 5 show the nonlinear load
signal can be reconstructed through a few modes that show high degrees of
predictability (near to periodic behavior). With the left column in Figure 5,
each mode has a narrow frequency band. This type of temporal modes are
intermittent patterns that play the key role in capturing features of strongly
nonlinear dynamics, but unlike periodic modes, they are not detectable by
classical singular spectrum analysis [Giannakis2012NLSA].

(a)  Figure 4: Examples of spatiotemporal load reconstruction. The load data
of 1494 stations during January 2014 are reconstructed using 100 Koopman
eigenfunctions.

(a)  (b)  (c)

(d)  (e)  (f)

(g)  (h)  (i)

(j)  (k)  (l)

(m)  (n)  (o)

Figure 5: Spatial and temporal patterns with multiple scales recovered by the
load data during January 2014. The left, middle, and right columns show,
respectively, the power spectrum, time evolution, and spatial pattern induced
by the Koopman eigenfunction. Ordered from above, the first row is associated
with the trivial Koopman eigenfunction ($\omega=0$), and the subsequent rows
are associated with frequencies scaled with a month ($3.86\times
10^{-7}\text{Hz}$), a week ($1.65\times 10^{-6}\text{Hz}$), a day ($1.15\times
10^{-5}\text{Hz}$), and a half day ($2.31\times 10^{-5}\text{Hz}$).

###  6.4 Power usage patterns detected by spatial modes

The right column in Figure 5 shows that each Koopman eigenfunction induces a
spatial pattern that is specific to that mode. Specifically, while the spatial
modes associated with zero and low frequencies (see the three upper rows in
Figure 5) are marked by the station in France, which also generates the
largest average value in Figure 3, the two last rows of Figure 5 indicate that
high-frequency modes scaled with a day and a half day lead to quite different
patterns. In particular, the mode scaled with a day shows two stations in
Northern Europe (Germany), as well as some other stations located in Southern
Europe (Italy, Greece, and Spain), which appear with the largest amplitude in
this mode. Moreover, we observe, from the last row of Figure 5, these southern
stations, after the one in France, generate the largest values in the mode
scaled with a half day. This indicates different stations generating power for
different purposes. Thus, stations that generate power for industrial purposes
have the largest values in low-frequency modes, while stations generating
power for residential and commercial areas mark the higher-frequency (day and
half day) modes!

Comparing different spatial modes of different months reveals how energy
pattern changes during the year. In Figure 6, we show spatial patterns
associated with the frequency scaled with a day during 12 months of 2014.
According to Figure 6, while the modes during the winter are marked by some
stations in some Southern Europe regions, they contribute less in the spatial
patterns generated from June 2014 to October 2014 (i.e., during the last month
of spring, entire summer, and first month of fall). Thus, these areas hold a
more uniform power usage during a day in warmer months, while their usage
changes significantly during a day in cold months.

(a)  (b)  (c)

(d)  (e)  (f)

(g)  (h)  (i)

(j)  (k)  (l)

Figure 6: Extracted spatial patterns showing the frequency scaled with a day.

###  6.5 Clustering the data

Our motivation for cluster-based analysis is that, a decomposition of the
Koopman operator based on the sparse structure of the underlying dynamical
system allows one to consider the system as a family of subsystems
interconnected by a graph where eigenfunctions for the subsystems induce
eigenfunctions for the whole system [Korda2022Sparsity]. While the original
study in [Korda2022Sparsity] requires a priori knowledge of the underlying
dynamics model to detect the sparsity and a decomposition of a dynamical
system, we devise a method based on PHATE [PHATE2019] to detect the sparsity
and the associated subsystems based on the available data. Moreover, we use a
different variant of the original approach in [PHATE2019] by using the
variable bandwidth kernels [Berry2016Variable] in the PHATE.

By applying the Algorithm 2 to the normalized data during 2012-2014, we obtain
the upper plot in Figure 8, where we decide for 10 clusters based on the
unified data. We identify the clusters using the $k$-means algorithm
[Bishop2006]. By looking into the lower plot in Figure 8, we observe a
clustering correlated with the physical graph of the power plants. Thus,
stations in different geographical areas fall into different clusters.
Moreover, we explore clusters with completely synchronized stations in them.
Specifically, the clusters 2, 3, 4, 6, and 7 contain stations that are
synchronized in their clusters. These stations are corresponding to the groups
that depart and unify distinctly from the main body in the upper plot of
Figure 8. Moreover, we observe more clusters (i.e., clusters 1 and 10) being
separated from the main body in upper Figure 8, so that we can extract more
groups of synchronized stations by either increasing the number of clusters in
$k$-means approach or re-clustering the current clusters. Although increasing
the number of clusters will lead to increased accuracy, we choose for 10
clusters and make a trade-off between accuracy and model complexity.
Therefore, with the proposed clustering approach, we are able to split the
large-scale system into several interconnected subsystems by identifying
different power generation patterns.

(a)

(b)

Figure 7: Visualizing the data on the low-dimensional PHATE coordinates
unveils 10 clusters.

###  6.6 Short-term prediction

The short-term load forecasting is an important basis for daily and weekly
dispatching of power grids [Kong2020DMD]. For the considered short-term
prediction, we perform a week-ahead forecasting based on Koopman modes of the
previous week. Therefore, we use the load data of a week to forecast the next
week load. In this way, we assume the Koopman modes do not change
significantly for two successive weeks. Thus, any considerable change to these
modes, e.g., due to weather conditions or holidays, may potentially a source
of error for the considered approach. However, we show that the proposed
approach, with its underlying assumption, outperforms a load forecasting
method based on deep learning.

For forecasting purpose, we use a hierarchical scheme (Figure 8) which works
based on clustered power stations using Algorithm 2. In this regard, we
forecast the load of the individual stations in each cluster by constructing a
single Koopman model for that cluster. Thus, based on our analysis of the load
data in Section 6.5, we construct 10 Koopman models.

Figure 8:  Cluster-based prediction of load.

The considered numerical evaluation of the Koopman operator requires tuning
several parameters. Specifically, the parameters regarding the number of
delays $Q$ and Koopman eigenfunctions $l$ play central roles in the
performance of the proposed approach. Regarding the number of delays $Q$, it
is the underlying factor to extract highly coherent dynamical behaviors, and
theoretically, in the limit of infinitely many delays we can recover the exact
dynamical system properties [Giannakis2019Delay; Giannakis2021Coherent]. We
also observe from Figure 1 that we will get more regular and decoupled
eigenfunctions by increasing the number of delays in the load dynamics.
However, we note that, the number of delays is limited by the availability of
the data in real systems. Moreover, for the finite-dimensional Koopman
subspaces, the delay embedding window shows some correlation with the
frequency ranges of interest [Giannakis2018Turbulent;
Giannakis2020Neuroscience; Giannakis2021Coherent; Giannakis2021Climate], so
that including uncorrelated delays may ignore the dynamical features in these
ranges. Therefore, in the following, we assess the sensitivity of the proposed
approach to the delay horizon $Q$.

To investigate the role of delay horizon $Q$ in load forecasting, we fix the
number of eigenfunctions and compute the RMSE for each station. Figure 9 shows
the results for 12 months of 2014, with $l=100$. As in Figure 9, the RMSE
values reach a minimum at $Q=200$ for most months, and generally they do not
change considerably by increasing $Q$ over 200, which shows the robustness of
the proposed approach for the important parameter $Q$. Considering that, by
the sampling rate of 1 hour, we have 168 samples over a week, this suggests
that a reasonable value for the delay is the one scaled with forecasting
horizon (which is one week in the considered short-time forecasting). Although
in a typical case where due to some special conditions, such as the weather or
substantially noisy data, we may get an increased accuracy by increasing the
delay horizon (for example November in Figure 9), still we care about the
availability of the data, as well as increased computational and model
complexities. Therefore, we suggest for a delay scaled with (equal to) the
forecasting horizon, and still remain cautious in the case of substantially
noisy data or abnormal weather conditions. In such cases, the proposed
approach gives the freedom to increase the delay horizon to get Koopman modes
that are more robust.

(a)  (b)  (c)  (d)

(e)  (f)  (g)  (h)

(i)  (j)  (k)  (l)

Figure 9: RMSE of different stations at each considered time delay ($l=100$).

The number of Koopman eigenfunctions included in the forecasting model is the
next parameter we investigate. Theoretically, the considered library of
eigenfunctions should be rich enough to capture the dynamical features in wide
spectra. However, this does not necessarily mean using arbitrarily large
libraries of Koopman eigenfunctions for real world data. Specifically, the
evaluation of Koopman eigenfunctions with higher energy levels is accompanied
by computational errors, and including them may bring false information and
instability into the model. In fact, the higher modes are not well-behaved on
the data manifold and may not be accurately approximated from finite datasets
[Giannakis2019Spectral]. Moreover, to improve the model interpretability and
avoid increased complexity, the use of modes that carry the least variance of
the data and are of no physical importance is not justified. Further
consideration regarding the number of Koopman eigenfunctions is overfitting in
the regression problem (28), which may arise when the library $\Psi$ is too
large for the given dataset [tavasoli2021koopman].

To illustrate the discussion in the load dynamics, Figure 10 indicates the
RMSE values for different values of $l$ when the delay horizon is $Q=168$. For
most months In Figure 10, the RMSE achieves a minimum (in terms of RMSE
distribution of individual stations and/or the mean value) for some
intermediate values of $l$ (around $l=50$). This indicates that increasing $l$
over some intermediate values gives a library that is too large so that
overfitting ensues.

Figure 11 shows the RMSE values for different stations and different months of
2014, when $Q=168$ and $l=50$. We thus use the same parameters for different
Koopman models of different clusters and months. Figure 11 reports the least
forecasting errors belong to the six months of winter and summer of 2014. The
error is generally larger in spring and fall, and after the extraordinary
error in April, Figure 11 reports the largest errors for the last month of
spring (June), and the two last months of fall (November and December).
Therefore, as we capture the dynamics to forecast, those features that are not
relevant to the dynamics and may be due to the significantly noisy data or
exogenous effects (such as the weather conditions in spring and fall or public
holidays) will be sources of forecasting error. Specifically, what contributes
to extraordinary errors in April 2014 is due to public holidays.

To shed more light, we show a load example in the two last weeks of April 2014
in a separate figure to investigate this month in more details (see Figure
12). We also show the RMSE for April 2014 for the increased delay window
$Q=500$. The load pattern in the left of Figure 12 indicates two public
holidays on April 18 (Friday) and April 22 (Monday). These two public holidays
significantly reduce the correlation between the training and test datasets,
i.e., between the two last weeks of April 2014. Furthermore, the right of
Figure 12 shows that a more accurate forecasting results, and the errors
effectively decrease, by increasing the delay to $Q=500$. Therefore, the role
of the delay embedding in understanding the complex load dynamics and
achieving more robust forecasting in the presence of exogenous effects is
established.

Figure 13 illustrates more performance of the proposed approach by indicating
examples of load forecasting in January 2014. We use the same parameters,
i.e., $Q=168$, $l=50$, and $\theta=10^{-9}$, for different Koopman models
regarding different clusters and different months in Figures 11 and 13.

(a)  (b)  (c)  (d)

(e)  (f)  (g)  (h)

(i)  (j)  (k)  (l)

Figure 10: RMSE of different stations at each considered number of Koopman
eigenfunctions ($Q=168hr$ (one-week delay)).

(a)  (b)  (c)

(d)  (e)  (f)

(g)  (h)  (i)

(j)  (k)  (l)

Figure 11: RMSE for different months with $Q=168hr$ (one-week delay) and
$l=50$.

(a)  (b)

Figure 12: A load example (left) and the RMSE (right) for April 2014 with
$Q=500$ ($l=50$).

(a)  (b)

(c)  (d)

(e)  (f)

(g)  (h)

(i)  (j)

Figure 13: Examples of forecasted load in last week of January 2014 ($Q=168hr$
(one-week delay), $l=50$.

###  6.7 Comparison with deep learning

We explore the system dynamics by unfolding and decomposing the data manifold
structure. Thus, instead of a straightforward search within the data in an
opaque and amorphous environment, we run a targeted use of the load data. We
have already shown that this significant approach yields a clear image of the
dynamics and large freedom for analysis and interpretations. To illustrate
further features of this approach, we discuss its advantage over an LSTM
setting that is the accepted tool in load forecasting (for examples, see
[Marino2016LSTM; Kong2019LSTM; Peng2021flexible; Zhu2022LSTM]). We show that,
the proposed method in this paper is more accurate and computationally
tractable.

Figure 14 compares the RMSE values for different clusters and different months
of 2014. For Koopman, we used $l=50$ eigenfunctions with a delay of $Q=168$.
Moreover, the LSTM uses two layer structures, with 100 hidden units in each
layer. In addition, in LSTM, the activation function is $\tanh$; Adam
optimizer is used for gradient descending; the loss function is MSE; the batch
size is 64, and the epochs is set to 10000. We used a dropout unit with rate
0.2 after each LSTM layer to avoid over-fitting. Both approaches perform a
short-term forecasting for the last week of each month, and the training data
is one-week ahead. The parameter values for the LSTM were tuned to yield the
least forecasting errors while avoiding overfitting, so that increasing the
number of LSTM layers, hidden units, or number of epochs over the considered
values led to overfitting. Moreover, we conduct a cluster-based forecasting,
shown in Figure 8, so that the PHATE cluster-based LSTM is introduced and
utilized for the first time in this work.

Therefore, while the LSTM uses a substantially complex network (two LSTM layer
structures with 100 hidden units in each layer), the Koopman approach uses as
few as 50 modes which are the most regular and predictable patterns of the
dynamics. Moreover, to forecast in each month, while for the Koopman approach
the operating computer took only some minutes (less than 10 minutes) to run,
the LSTM was executed in more than 6 hours with the same computer.

Overall, Figure 14 is marked by the pattern that reports less forecasting
error by the Koopman approach. The main exception to this rule is in April.
Following our discussion in Section 6.6, we know that some exogenous effects
(public holidays) in the second last week of April 2014 lead to the modes that
are less coherent for the delay horizon $Q=168$, and we may increase the delay
to get more robust forecasting. Thus, Figure 15 (in the left) shows that,
increasing the delay to $Q=500$, we are able to get a forecasting that is
comparable to the LSTM. Figure 15 further illustrates that better forecasting
results for the clusters 3 and 7 on November and December by widening the
embedding window (compare clusters 3 and 7 in the middle and right of Figure
15 with their corresponding plots in Figure 14).

(a)  (b)  (c)  (d)

(e)  (f)  (g)  (h)

(i)  (j)  (k)  (l)

Figure 14: RMSE values due to Koopman and LSTM for different clusters in the
last week of each month in 2014.

(a)  (b)  (c)

Figure 15: Comparing LSTM and the Koopman with $Q=500$.

###  6.8 Modal noise and innovation noise

Expanding the load dynamics in terms of Koopman eigenfunctions (see Equation
(27)) enables us to distinguish different sources of forecasting error, and
identify modal and innovation noises [Mezic2022noise]. Specifically, since we
consider a truncation of the infinite-dimensional Koopman operator, the
finite-dimensional expansion $C\psi(x)$ leads to a residual $r$ with respect
to the real load data $x$; that is $x(t)=C\psi(t)+r(t)$. Consequently we split
the residual $r$ into the modal noise $\eta_{m}$ and innovation noise
$\eta_{i}$ as following:

$\begin{gathered}\eta_{m}(t)=CC^{\dagger}r(t)\\\
\eta_{i}(t)=(I-CC^{\dagger})r(t)\end{gathered}$ (33)

where $I$ is the identity matrix and $C^{\dagger}$ is the Moore-Penrose
pseudo-inverse of $C$. Therefore, $\eta_{m}$ is in the plane of truncated
Koopman eigenfunctions, while $\eta_{i}$ remains perpendicular to this plane.
The in-plane modal noise is useful to get a confidence bound for the
prediction model. Furthermore, while the innovation noise decreases by
increasing the number of eigenfunctions and extending the modal plane, the
modal noise is essentially due to mismatch of the regression matrix $C$
between training and test data and is not tunable by increasing the number of
Koopman eigenfunctions. This sets forth a criterion for the optimal set
(number) of Koopman eigenfunctions [Mezic2022noise].

To expand the discussion in power grids, we examine the modal and innovation
noise for January 2014. Moreover, we consider two cases of with and without
measurement noise. In each case, we compute the residuals for all 168 hours of
the last week of January 2014 and for all stations. Figures 16 shows the case
without measurement noise, so that we considered the original data without
adding noise. While a significant modal noise is held by the upper row of
Figure 16, the innovation noise in the lower row decreases by increasing the
number of Koopman eigenfunctions. This indicates that, only a few number of
Koopman eigenfunctions are sufficient to capture the effective dynamics modes
in power grids, and the residual (error) remains in the space of these modes.
Therefore, any improvement in the forecasting may be possible by modifying the
matrix $C$ and accurately assessing the contribution of Koopman modes in the
load dynamics evolution. For example, a modification of the matrix $C$ might
be obtained based on exogenous effects such as weather conditions and public
holidays. In other words, the contribution of Koopman eigenfunctions in the
evolution of the load dynamics of one week (training data) is modified for the
next week (test data) based on exogenous effects.

Next, Figure 17 shows the case with measurement noise, where we have manually
added 10$\%$ noise into the original data. Here, the innovation noise changes
slightly by increasing the number of Koopman eigenfunctions in the lower row
of Figure 17. Furthermore, by the Koopman approach, we are able to distinguish
the parts of the error due to regression mismatch (upper row of 17) and due to
measurement noise (lower row of 17). In fact, since the added noise holds no
dynamics, it does not fall into the span of the Koopman eigenfunctions. Hence,
by distinguishing the components of the dynamics (with low energy) and the
noisy parts (with high energy), the Koopman model puts no effort to capture
the noise dynamics.

(a)  (b)  (c)  (d)

(e)  (f)  (g)  (h)

Figure 16: Original data (without added noise): Modal (upper row) and
innovation (lower row) noise for load forecasting of all stations during 168
hours of last week of January 2014. Ordered from left, the columns show the
results for $l=5,10,20,50$ Koopman eigenfunctions.

(a)  (b)  (c)  (d)

(e)  (f)  (g)  (h)

Figure 17: 10$\%$ noise included: Modal (upper row) and innovation (lower row)
noise for load forecasting of all stations during 168 hours of last week of
January 2014. Ordered from left, the columns show the results for
$l=5,10,20,50$ Koopman eigenfunctions.

##  7 Conclusion

In this paper, we proposed data-driven identification of power grids load
dynamics in the operator-theoretic framework and, for this purpose, we
considered the Koopman operator approach and leveraged its connection with a
generation of modern dynamics-adapted kernel operators in delay coordinates.
In this framework, we identify a reduced-order model for the load dynamics by
extracting a set of coherent spatiotemporal patterns based on a Koopman mode
decomposition. The coherent patterns are the most robust features of the load
dynamics and each pattern recovers a dynamical feature that evolves linearly
with an associated time scale. Therefore, coherent patterns yield a set of
coordinate transformations that leads to a linear representation of the
dynamics, and the linear model is valid over the whole basin. Moreover, they
lead to a spatiotemporal spectral analysis in an optimal basis that is
intrinsic to the nonlinear load dynamics, and by studying the corresponding
spectra we are able to discover various basic energy patterns in multiple time
scales. We implement a cluster-based Koopman approach for the spatiotemporal
load forecasting over large-scale power grids where by the proposed clustering
approach, based on PHATE, we are able to detect the synchronized load profiles
and construct a single forecasting model for them. The advantages of the
proposed approach are discussed by implementing for a large-scale dataset
regarding a renewable electric power system within the continental European
electricity system, and comparing to a deep learning approach. In particular,
we suggest an interpretable machine learning approach and present a
transparent image of the load dynamics by discovering its basic underlying
components. Yet all dynamical features are extracted directly from the
available data in a nonparametric setting, without any assumption or
restricting framework for analyzing the load data. In this aspect, the
proposed approach creates and resides within a bridge connecting the physics-
based and data-driven statistical and deep learning approaches, with the
former carrying limited practical features that may be undermined due to
requiring all properties and environmental parameters prior to any analysis,
and the later defying straightforward analysis and interpretations. Compared
to deep architectures for load forecasting, the proposed approach is
independent of rigorous training phases, hand-crafted feature selection, or
optimal parameter selection, and, in terms of computational speed and
complexity, it is more suited for the size of the data in large-scale power
grids. The forecasting is funded based on the features that are due to some
dynamics, so that least effort is made for capturing random features, and
embedding the data into a space of delay coordinates leads to coherent
patterns that encode all (hidden) dynamical features of interconnected
subsystems and exogenous effects.

A natural extension of this approach will be improving the extracted coherent
patterns, since they are our key to forecasting and analyzing. The first
direction will be realizing the patterns in the presence of varying dynamical
features, where Koopman modes may change, for example, from one week to
another, due to abnormal (out of dynamics) weather conditions, holidays, or
any other change in the interconnected dynamics and exogenous effects. This
may be done by a time-varying realization of Koopman eigenfunctions and by
extending the existing online DMD approaches [Rowley2019onlineDMD;
Mustaffa2020onlineDMD]. The second direction to improve the proposed approach
is distilling the set of extracted coherent patterns by detecting the spectral
pollution, which is the appearance of spurious eigenvalues due to numerical
errors of discretization and which have nothing to do with the underlying
dynamics [Colbrook2021].