# Spectral Koopman Method for Identifying Stability Boundary

Bhagyashree Umathe, and Umesh Vaidya

###### Abstract

The paper is about characterizing the stability boundary of an autonomous
dynamical system using the Koopman spectrum. For a dynamical system with an
asymptotically stable equilibrium point, the domain of attraction constitutes
a region consisting of all initial conditions attracted to the equilibrium
point. The stability boundary is a separatrix region that separates the domain
of attraction from the rest of the state space. For a large class of dynamical
systems, this stability boundary consists of the union of stable manifolds of
all the unstable equilibrium points on the stability boundary. We characterize
the stable manifold in terms of the zero-level curve of the Koopman
eigenfunction. A path-integral formula is proposed to compute the Koopman
eigenfunction for a saddle-type equilibrium point on the stability boundary.
The algorithm for identifying stability boundary based on the Koopman
eigenfunction is attractive as it does not involve explicit knowledge of
system dynamics. We present simulation results to verify the main results of
the paper.

##  I Introduction

Characterizing stability and identifying stability boundary of nonlinear
dynamical systems is of interest to various problems in engineering and
science. Examples include biological systems, biomedicine, robotics, power
systems, power electronics, and economics. Given the significance of the
problem, there is a long history and a considerable body of literature on the
stability theory of systems. Methods for characterizing stability and
stability boundary of equilibrium dynamics can be broadly classified into two
different classes: Lyapunov and non-Lyapunov-based methods [1].

The Lyapunov-based methods rely on constructing the Lyapunov or energy
function for stability verification and providing an estimate for the
stability boundary in the form of a domain of attraction. Roughly speaking,
the domain of attraction is a set of all initial conditions attracted to
equilibrium dynamics asymptotically [2]. One of the main challenges with the
Lyapunov-based methods is that there are no systematic methods to construct
these functions, and the stability boundary estimates provided are often
conservative. More recent work along these lines includes the development of
the maximal Lyapunov function, optimal estimation of stability region based on
a given Lyapunov function, Linear Matrix Inequality, and Sum of Squares-based
optimization methods for constructing Lyapunov functions [3]. These
optimization-based methods suffer from the curse of dimensionality and entail
high computational efforts with the increase in the dimension of the state
space. With the recent advances in computing and neural networks, there are
also efforts to provide approaches based on Deep Neural Networks for
constructing Lyapunov function [4].

The non-Lyapunov-based approach for stability characterization relies on
exploiting the geometrical structure in the form of stable and unstable
manifolds of the dynamical system for identifying stability boundary. In
particular, for a large class of dynamical systems, the stability boundary can
be characterized using the stable manifold of an unstable equilibrium point on
the stability boundary [5].

More recently, linear operator theoretic methods involving Perron-Frobenius
and Koopman operators have become popular for analyzing and synthesizing
nonlinear dynamical systems [6, 7, 8, 9, 10]. In [11], authors have proposed
uncertainty propagation using Koopman spectrum-based reachability analysis.
These operators’ linearity and spectral properties are used to construct
stability certificates in the form of Lyapunov functions and Lyapunov
measures. Data-driven approaches are proposed for constructing the Lyapunov
function using the Koopman operator. In this paper, we exploit the spectral
properties of the Koopman operator to discover a non-Lyapunov-based approach
for characterizing the stability boundary. This is made possible by using the
fact that the spectral properties of the Koopman operator have an intimate
connection with the state space geometry of the dynamical systems. In
particular, a dynamical system’s stable and unstable manifolds can be obtained
as a zero-level curve of the eigenfunctions of the Koopman operator. We
exploit this insight towards the development of systematic methods for the
identification of stability boundary. There are some distinct advantages of
using the Koopman spectrum to develop a non-Lyapunov approach for identifying
stability boundary. We will highlight these advantages as part of the main
contributions of this paper.

The main contributions of this paper are as follows. We characterize the
stability boundary in terms of the principal eigenfunction of the Koopman
operator. In particular, the stability boundary is characterized by the zero-
level curve of the eigenfunction corresponding to a positive eigenvalue. A
path-integral formula is provided to compute the Koopman principal
eigenfunction corresponding to the positive eigenvalue associated with a
saddle-type unstable equilibrium point on the stability boundary. The path-
integral formula is used to compute the value of eigenfunction at a few sample
data points in the state space. We present convergence analysis results for
the approximation of the Koopman eigenfunction based on the value of the
eigenfunction at a few sample data points. One of the main advantages of the
proposed Koopman-based approach for characterizing stability boundary compared
to using the Koopman operator for the computation of the Lyapunov function for
stability assessment is that we are only required to compute one
eigenfunction. This contrasts the computation of all the dominant
eigenfunctions used to construct the Lyapunov function. This is because the
stability boundary forms a manifold of co-dimension one, characterized as a
zero-level curve of one eigenfunction. Another advantage unique to our
proposed computational approach is that we can compute the Koopman
eigenfunction directly from the data instead of obtaining it from the
approximation of the Koopman operator. This is a tremendous computational
advantage for systems with large dimension.

##  II Preliminaries and Notations

Notations: ${\mathbb{R}}^{n}$ denotes the $n$ dimensional Euclidean space. We
denote ${\cal L}_{\infty}({\mathbf{X}}),{\cal C}^{1}({\mathbf{X}})$ as the
space of all essentially bounded real-valued functions and continuously
differentiable functions on ${\mathbf{X}}\subseteq{\mathbb{R}}^{n}$
respectively. Here, ${\mathbf{s}}_{t}({\mathbf{x}})$ is the flow of
$\dot{\mathbf{x}}={\mathbf{f}}({\mathbf{x}})$, at time $t$ with initial
condition ${\mathbf{x}}$.

In this paper, we are interested in characterizing the stability boundary
using the Koopman eigenfunctions for an autonomous dynamical system of the
form.

$\displaystyle\dot{\mathbf{x}}={\mathbf{f}}({\mathbf{x}}),\quad{\mathbf{x}}\in{\mathbf{X}}\subseteq{\mathbb{R}}^{n}$
(1)

where, ${\mathbf{X}}$ is assumed to a smooth manifold without boundary or an
open subset of ${\mathbb{R}}^{n}$. It is assumed that the vector field is
smooth enough to ensure the existence and uniqueness of the solution of the
differential equation. We say an equilibrium point ${\mathbf{x}}^{\star}$ is
hyperbolic if, ${\mathbf{A}}:=\frac{\partial{\bf
f}}{\partial{\mathbf{x}}}({\mathbf{x}}^{\star})$ has no eigenvalues on the
$j\omega$ axis. The following assumption is made on the vector field (1) in
the rest of the paper.

###### Assumption 1

We assume that in system (1), ${\mathbf{f}}({\mathbf{x}})$ is atleast ${\cal
C}^{1}$ function of ${\mathbf{x}}$ with an asymptotically stable equilibrium
point at ${\mathbf{x}}_{s}$, and all the other equilibrium points hyperbolic.

###  II-A Characterization of the Stability Boundary

We start with definitions related to the stability analysis.

######  Definition 1 (Stability Region and Boundary)

For the system (1), with the flow ${\mathbf{s}}_{t}({\mathbf{x}})$, let
${\mathbf{x}}_{s}$ be an asymptotically stable equilibrium point. We define
the region ${\cal A}({\mathbf{x}}_{s})$ as the stability region or the region
of attraction of the stable equilibrium point ${\mathbf{x}}_{s}$ as

$\displaystyle{\cal
A}({\mathbf{x}}_{s}):=\left\\{{\mathbf{x}}\in{\mathbb{R}}^{n}:\quad\lim_{t\rightarrow\infty}{\mathbf{s}}_{t}({\mathbf{x}})={\mathbf{x}}_{s}\right\\}.$
(2)

The boundary of ${\cal A}({\mathbf{x}}_{s})$ is called the stability boundary
$\partial{\cal A}({\mathbf{x}}_{s})$, which is a separatrix of
${\mathbf{x}}_{s}$.

###### Definition 2

[Local Stable and Local Unstable Manifold] We consider ${\mathbf{x}}^{\star}$
as a hyperbolic equilibrium point, and $D\subset{\mathbb{R}}^{n}$ is a
neighborhood of ${\mathbf{x}}^{\star}$. If the flow of (1) is
${\mathbf{s}}_{t}({\mathbf{x}})$ then the local stable manifold of
${\mathbf{x}}^{\star}$ is given by:

$\displaystyle W^{s}_{l}({\mathbf{x}}^{\star}):=\\{{\mathbf{x}}\in
D:{\mathbf{s}}_{t}({\mathbf{x}})\rightarrow{\mathbf{x}}^{\star}\quad\text{as}\quad
t\rightarrow\infty\\}$ (3)

Similarly, the local unstable manifold is given by:

$\displaystyle W^{u}_{l}({\mathbf{x}}^{\star}):=\\{{\mathbf{x}}\in
D:{\mathbf{s}}_{t}({\mathbf{x}})\rightarrow{\mathbf{x}}^{\star}\quad\text{as}\quad
t\rightarrow-\infty\\}$ (4)

###### Definition 3

[Stable and Unstable Manifold] The stable manifold
$W^{s}({\mathbf{x}}^{\star})$ and Unstable manifold
$W^{u}({\mathbf{x}}^{\star})$ are derived by the local manifolds with the flow
$s_{t}({\mathbf{x}})$ in backward and forward in time respectively, as:

$\displaystyle W^{s}({\mathbf{x}}^{\star}):=\bigcup_{t\leq
0}{\mathbf{s}}_{t}\left(W^{s}_{l}({\mathbf{x}}^{\star})\right),\;W^{u}({\mathbf{x}}^{\star}):=\bigcup_{t\geq
0}{\mathbf{s}}_{t}\left(W^{u}_{l}({\mathbf{x}}^{\star})\right).$

For the characterization of the stability boundary, the following assumptions
are made on the vector field (1).

###### Assumption 2

  1. A1.

All the equilibrium points on the stability boundary are hyperbolic.

  2. A2.

The stable and unstable manifolds of the equilibrium points on the stability
boundary satisfy the transversality assumption (The intersection of two
manifolds at point ${\mathbf{x}}$ is said to be transversal if the tangent
space of the two manifolds at point ${\mathbf{x}}$ span the entire space).

  3. A3.

Every trajectory on the stability boundary approaches one of the equilibrium
points at $t\to\infty$.

The following theorem from [1, Theorem 4.8] characterizes the stability
boundary under the above assumption.

###### Theorem 1

For the nonlinear dynamical system (1), satisfying Assumptions 1 and 2, let
${\mathbf{x}}_{i},\;\;i=1,\ldots N$ be the equilibrium points on the stability
boundary $\partial{\cal A}({\mathbf{x}}_{s})$ of the stable equilibrium point
${\mathbf{x}}_{s}$. Then

  1. a)

${\mathbf{x}}_{i}\in\partial{\cal A}({\mathbf{x}}_{s})$ if and only if
$W^{u}({\mathbf{x}}_{i})\cap\partial{\cal A}({\mathbf{x}}_{s})\neq\emptyset$.

  2. b)

$\partial{\cal A}({\mathbf{x}}_{s})=\cup_{i}W^{s}({\mathbf{x}}_{i})$.

Assumptions A1 and A2 are generic. Roughly speaking, we say a property is
generic for a class of systems if that property is true for almost all systems
in the class. Assumption A3 is not a generic property; thus, it needs to be
verified. The existence of the Lyapunov function provides a sufficient
condition for the Assumption A3 to hold. In this paper, we pursue a more
direct approach for computing the stability boundary for the class of system
satisfying Assumption 2. The proposed approach relies on directly computing
the stability boundary based on constructing stable manifolds of the
equilibrium points on the stability boundary. So, we study the class of system
for which the Lyapunov-based approach is used to determine the domain of
attraction as the existence of the Lyapunov function is sufficient for the
Assumption 2 to be satisfied.

###### Definition 4

An equilibrium is said to be type-one if exactly one eigenvalue of the
linearization of the system at that equilibrium point has positive real part
and called source if all the eigenvalues have positive real part.

The following theorem from [1, Theorem 4.10] characterizes the structure of
the equilibrium on the stability boundary.

###### Theorem 2

For the nonlinear dynamical system (1) containing two or more stable
equilibrium points, if the system satisfies Assumption 2, then the stability
boundary $\partial{\cal A}({\mathbf{x}}_{s})$ of the stable equilibrium point
${\mathbf{x}}_{s}$ must contain at least one type-one equilibrium point. If,
furthermore, the stability region ${\cal A}({\mathbf{x}}_{s})$ is bounded,
then $\partial{\cal A}({\mathbf{x}}_{s})$ must contain at least one type-one
equilibrium point and one source.

###  II-B Koopman Operator and its Spectrum

In this section, we provide a brief overview of existing results on the
spectral theory of the Koopman operator. For more details on this topic, refer
to [12, Chapter 7].

######  Definition 5 (Koopman Operator)

$\mathbb{U}_{t}:{\cal L}_{\infty}({\mathbf{X}})\to{\cal
L}_{\infty}({\mathbf{X}})$ for dynamical system (1) is defined as

$\displaystyle[\mathbb{U}_{t}\psi]({\mathbf{x}})=\psi({\mathbf{s}}_{t}({\mathbf{x}})),\;\;\;\psi({\mathbf{x}})\in{\cal
L}_{\infty}({\mathbf{X}}).$ (5)

The infinitesimal generator for the Koopman operator is

$\displaystyle\lim_{t\to
0}\frac{(\mathbb{U}_{t}-I)\psi}{t}=\frac{\partial\psi}{\partial{\mathbf{x}}}{{\mathbf{f}}}({\mathbf{x}})=:{\cal
K}_{{\mathbf{f}}}\psi,\;\;t\geq 0.$ (6)

###### Definition 6

[Eigenvalues and Eigenfunctions of Koopman] A function
$\phi_{\lambda}({\mathbf{x}})\in{\cal C}^{1}({\mathbf{X}})$ is said to be
eigenfunction of the Koopman operator associated with eigenvalue $\lambda$ if

$\displaystyle[{\mathbb{U}}_{t}\phi_{\lambda}]({\mathbf{x}})=e^{\lambda
t}\phi_{\lambda}({\mathbf{x}}).$ (7)

Using the Koopman generator, (7) can be written as

$\displaystyle\frac{\partial\phi_{\lambda}}{\partial{\mathbf{x}}}{{\mathbf{f}}}({\mathbf{x}})=\lambda\phi_{\lambda}({\mathbf{x}}).$
(8)

where
$\frac{\partial\phi_{\lambda}}{\partial{\mathbf{x}}}=(\frac{\partial\phi_{\lambda}}{\partial
x_{1}},\ldots,\frac{\partial\phi_{\lambda}}{\partial x_{n}})$ is a row vector.
The spectrum of the Koopman operator, in general, is very complex and could
consist of discrete and continuous parts. In this paper, we are interested in
approximating the eigenfunctions of the Koopman operator with associated
eigenvalues, same as that of the linearization of the nonlinear system at the
equilibrium point. With the hyperbolicity assumption on the equilibrium point
of (1), this part of the spectrum of interest is known to be discrete and
well-defined.

In the following, we briefly summarize the results from [13]. Equations (7)
and (8) provide a general definition of the Koopman spectrum. However, the
spectrum can be defined over finite time or a subset of the state space.

######  Definition 7 (Open Eigenfunction [13])

Let $\phi_{\lambda}:{\boldsymbol{\mathcal{C}}}\to{\mathbb{C}}$, where
${\boldsymbol{\mathcal{C}}}\subset{\mathbf{X}}$ is not an invariant set. Let
${\mathbf{x}}\in{\boldsymbol{\mathcal{C}}}$, and
$\tau\in(\tau^{-}({\mathbf{x}}),\tau^{+}({\mathbf{x}}))=I_{\mathbf{x}}$, a
connected open interval such that
$\tau({\mathbf{x}})\in{\boldsymbol{\mathcal{C}}}$ for all $\tau\in
I_{\mathbf{x}}$. If

$[{\mathbb{U}}_{\tau}\phi_{\lambda}]({\mathbf{x}})=\phi_{\lambda}({\mathbf{s}}_{\tau}({\mathbf{x}}))=e^{\lambda\tau}\phi_{\lambda}({\mathbf{x}}),\;\;\;\;\forall\tau\in
I_{\mathbf{x}},$

then $\phi_{\lambda}({\mathbf{x}})$ is called the open eigenfunction of the
Koopman operator family ${\mathbb{U}}_{t}$, for $t\in{\mathbb{R}}$ with
eigenvalue $\lambda$.

Based on the above definition, next, we define the concepts of subdomain
eigenfunction, principal eigenfunctions, and Koopman spectrum. If
${\boldsymbol{\mathcal{C}}}$ is a proper invariant subset of ${\mathbf{X}}$,
in which case $I_{\mathbf{x}}={\mathbb{R}}$ for every
${\mathbf{x}}\in{\boldsymbol{\mathcal{C}}}$, then $\phi_{\lambda}$ is called
the subdomain eigenfunction. If ${\boldsymbol{\mathcal{C}}}={\mathbf{X}}$ then
$\phi_{\lambda}$ will be the ordinary eigenfunction associated with eigenvalue
$\lambda$ as defined in (7). The open eigenfunctions, as defined above, can be
extended from ${\boldsymbol{\mathcal{C}}}$ to a larger reachable set when
${\boldsymbol{\mathcal{C}}}$ is open based on the construction procedure
outlined in [Definition 5.2, Lemma 5.1 [13]]. Let ${\cal P}$ be that domain.
The eigenvalues of the linearization of the system dynamics at the origin,
i.e., ${\mathbf{A}}$, will form the eigenvalues of the Koopman operator [13,
Proposition 5.8]. Our interest will be in constructing the corresponding
eigenfunctions defined over the domain ${\cal P}$. We will call these as
principal eigenfunctions. The principal eigenfunctions will be defined over a
proper subset ${\cal P}$ of the state space ${\mathbf{X}}$ (called subdomain
eigenfunctions) or over the entire ${\mathbf{X}}$ [13, Lemma 5.1, Corollary
5.1, 5.2, and 5.8].

The spectrum of the Koopman operator reveals important information about the
state space geometry of the dynamical system [13, 14]. In particular, we have
the following results.

###### Corollary 1

[13, Corollary 5.10] Let ${\mathbf{x}}^{\star}$ be the hyperbolic equilibrium
point of the system (1) with $\lambda_{1},\ldots,\lambda_{p}$ the eigenvalues
of the linearization of system (1) at ${\mathbf{x}}^{\star}$. Let
$\lambda_{1},\ldots,\lambda_{u}$ be eigenvalues with positive real part with
associated open eigenfunctions $\phi_{\lambda_{1}},\ldots,\phi_{\lambda_{u}}$
and $\lambda_{u+1},\ldots,\lambda_{p}$ be eigenvalues with negative real part
with associated open eigenfunctions
$\phi_{\lambda_{u+1}},\ldots,\phi_{\lambda_{p}}$ defined over the set ${\cal
P}$. Then, the joint level set of the eigenfunctions

$\displaystyle W_{{\cal
P}}^{s}=\\{{\mathbf{x}}\in{\mathbf{X}}:\phi_{\lambda_{1}}({\mathbf{x}})=\ldots=\phi_{\lambda_{u}}({\mathbf{x}})=0\\},$
(9)

forms the stable manifold on ${\cal P}$ and the joint level set of the
eigenfunctions

$\displaystyle W_{\cal
P}^{u}=\\{{\mathbf{x}}\in{\mathbf{X}}:\phi_{\lambda_{u+1}}({\mathbf{x}})=\ldots=\phi_{\lambda_{p}}({\mathbf{x}})=0\\},$
(10)

is the unstable manifold on ${\cal P}$ of origin equilibrium point.

The following remark on the eigenfunctions of the Koopman operator and stable,
unstable subspaces corresponding to a linear system is easy to prove.

###### Remark 1

The Koopman eigenfunctions corresponding to the eigenvalue $\lambda_{j}$ of
the linear dynamics, $\dot{\mathbf{x}}={\mathbf{A}}{\mathbf{x}}$, are given by
$\phi_{\lambda_{j}}({\mathbf{x}})={\mathbf{v}}_{j}^{\top}{\mathbf{x}}$, where
${\mathbf{v}}_{j}^{\top}$ is left eigenvector of ${\mathbf{A}}$ with
eigenvalue $\lambda_{j}$ i.e.,
${\mathbf{v}}^{\top}_{j}{\mathbf{A}}=\lambda_{j}{\mathbf{v}}^{\top}_{j}$.

##  III Stability Boundary computation using Koopman Spectrum

This section presents the main results on the computation of stability
boundary using the Koopman spectrum.

###### Theorem 3

Consider the dynamical system (1) satisfying Assumption 2. Let
${\mathbf{x}}_{i}$ for $i=1,\ldots,N$ be the type-one equilibrium point on the
stability boundary of the stable equilibrium point ${\mathbf{x}}_{s}$. Let
$\phi_{u}^{i}({\mathbf{x}})$ be the principal eigenfunctions corresponding to
the positive eigenvalues of the linearization
${\mathbf{A}}_{i}:=\frac{\partial{\mathbf{f}}}{\partial{\mathbf{x}}}({\mathbf{x}}_{i})$
with positive real part and defined in the domain ${\cal P}$. Then, the
stability boundary can be characterized in ${\cal P}$ using the Koopman
eigenfunction as the joint zero-level set as follows

$\displaystyle{\cal
S}_{b}:=\bigcup_{i=1}^{N}\\{{\mathbf{x}}:\phi^{i}_{u}({\mathbf{x}})=0\\}$ (11)

###### Proof:

The proof follows by applying results of Theorem 1, where it is shown that the
stability boundary is formed as the union of stable manifold of unstable
equilibrium points on the boundary i.e., $\partial{\cal
A}({\mathbf{x}}_{s})=\cup_{i=1}^{N}W^{s}({\mathbf{x}}_{i})$. From results of
Corollary 1, we know that stable manifold restricted to $\cal P$ is given by
$W^{s}_{\cal
P}({\mathbf{x}}_{i})=\\{{\mathbf{x}}:\phi_{u}^{i}({\mathbf{x}})=0\\}$. The
desired result is then followed by combining these two results. ∎

The above theorem provides the characterization of stability boundary in terms
of Koopman eigenfunction. However, computing the Koopman eigenfunction is a
challenge. We now present our main results for the computation of the Koopman
eigenfunction corresponding to the eigenvalue with a positive real part. With
no loss of generality, we will assume that the unstable equilibrium point on
the stability boundary is at the origin. Since the equilibrium point is
hyperbolic (Assumption 1 and 2), the system (1) admits the following
decomposition into linear and nonlinear parts.

$\displaystyle\dot{\mathbf{x}}={\mathbf{f}}({\mathbf{x}})={\mathbf{A}}{\mathbf{x}}+{\mathbf{F}}_{n}({\mathbf{x}})$
(12)

where ${\mathbf{A}}:=\frac{\partial{\mathbf{f}}}{\partial{\mathbf{x}}}(0)$ and
${\mathbf{F}}_{n}({\mathbf{x}}):={\mathbf{f}}({\mathbf{x}})-{\mathbf{A}}{\mathbf{x}}$.
Our goal is approximate the eigenfunction corresponding to eigenvalue
$\lambda$ of ${\mathbf{A}}$ with ${\rm Re}(\lambda)>0$. Similar to the
decomposition of the system dynamics, the principal eigenfunction also admits
a decomposition into linear and nonlinear parts. Let
$\phi_{\lambda}({\mathbf{x}})$ be the principal eigenfunction for eigenvalue
$\lambda$, we have

$\displaystyle\phi_{\lambda}({\mathbf{x}})={\mathbf{w}}_{\lambda}^{\top}{\mathbf{x}}+h_{\lambda}({\mathbf{x}}).$
(13)

where ${\mathbf{w}}_{\lambda}^{\top}{\mathbf{x}}$ is the linear part and
$h_{\lambda}({\mathbf{x}})$ is the nonlinear part. Substituting the above
expression of eigenfunction in (8) and using (12), we obtain
${\mathbf{w}}_{\lambda}^{\top}{\mathbf{A}}=\lambda{\mathbf{w}}_{\lambda}^{\top}$
i.e., ${\mathbf{w}}$ is the left eigenvector of ${\mathbf{A}}$ with eigenvalue
$\lambda$. We also obtain the linear partial differential equation to be
satisfied by $h({\mathbf{x}})$ as,

$\displaystyle\frac{\partial
h_{\lambda}}{\partial{\mathbf{x}}}{\mathbf{f}}({\mathbf{x}})-\lambda
h_{\lambda}({\mathbf{x}})+{\mathbf{w}}^{\top}{\mathbf{F}}_{n}({\mathbf{x}})=0$
(14)

In [15], we provided a path-integral formula for the computation of
$h_{\lambda}({\mathbf{x}})$ and hence the Koopman eigenfunction. The results
in [15] were developed to approximate eigenfunction, assuming the system has a
stable equilibrium point. The results from [15] do not apply to our case, as
we are dealing with the unstable or saddle-type equilibrium point on the
stability boundary. We start with the following results on the general
solution formula for the linear PDE (14).

###### Theorem 4

The solution formula for the first order linear PDE (14) can be written as

$\displaystyle h_{\lambda}({\mathbf{x}})=e^{-\lambda
t}h_{\lambda}({\mathbf{s}}_{t}({\mathbf{x}}))+\int_{0}^{t}e^{-\lambda
t}{\mathbf{w}}^{\top}_{\lambda}{\mathbf{F}}_{n}({\mathbf{s}}_{\tau}({\mathbf{x}}))d\tau$
(15)

###### Proof:

The PDE (14) can be written as

$\displaystyle\frac{dh_{\lambda}({\mathbf{s}}_{t}({\mathbf{x}}))}{dt}-\lambda
h_{\lambda}({\mathbf{s}}_{t}({\mathbf{x}}))+{\mathbf{w}}^{\top}_{\lambda}{\mathbf{F}}_{n}({\mathbf{s}}_{t}({\mathbf{x}}))=0.$
(16)

Multiplying throughout by $e^{-\lambda t}$, we obtain $\frac{d(e^{-\lambda
t}h_{\lambda}({\mathbf{s}}_{t}({\mathbf{x}})))}{dt}+e^{-\lambda
t}{\mathbf{w}}^{\top}_{\lambda}{\mathbf{f}}_{n}({\mathbf{s}}_{t}({\mathbf{x}}))=0$.
Next, by integrating the above from $0$ to $t$, we obtain the desired result.

∎

Before proving the main results we have following remark.

###### Remark 2

We are not interested in computing the eigenfunction per se but the zero-level
curve of the eigenfunction characterizing the stability boundary. Hence, our
interest will be in approximating the Koopman eigenfunction in the region
containing the stable manifold. Following the Definitions 2 and 3 of local and
global stable manifolds, we notice that the region containing the stable
manifold can be obtained by backward propagating a small set around the origin
containing the local stable manifold. This observation is crucial to
approximate eigenfunction, as the following theorem proves.

###### Theorem 5

Let $h_{\lambda}({\mathbf{x}})$ be continuous function of ${\mathbf{x}}$ and
$\cal U$ be the set such that for any point ${\mathbf{x}}\in\cal U$ there
exists a time $t({\mathbf{x}})$ for which
${\mathbf{s}}_{t}({\mathbf{x}})\in{\cal U}_{\epsilon}$, the $\epsilon$
neighborhood of the origin equilibrium point. Let $h_{\lambda}({\mathbf{x}})$
be the value of eigenfunction at point ${\mathbf{x}}\in\cal U$, then

$\displaystyle\left|h_{\lambda}({\mathbf{x}})-\int_{0}^{t({\mathbf{x}})}e^{-\lambda\tau}{\mathbf{w}}_{\lambda}^{\top}{\mathbf{F}}_{n}({\mathbf{s}}_{t}({\mathbf{x}}))d\tau\right|\leq\delta(\epsilon)$
(17)

for some $\delta(\epsilon)>0$ and all ${\mathbf{x}}\in\cal U$.

###### Proof:

Using result of Theorem 4, we know $h_{\lambda}({\mathbf{x}})$ satisfies (15).
Consider any point ${\mathbf{x}}\in\cal U$, and let
${\mathbf{y}}={\mathbf{s}}_{t({\mathbf{x}})}({\mathbf{x}})\in{\cal
U}_{\epsilon}$. Subsuiting
${\mathbf{y}}={\mathbf{s}}_{t({\mathbf{x}})}({\mathbf{x}})$ in (15), we obtain

$h_{\lambda}({\mathbf{x}})-\int_{0}^{t({\mathbf{x}})}e^{-\lambda\tau}{\mathbf{w}}_{\lambda}^{\top}{\mathbf{F}}_{n}({\mathbf{s}}_{\tau}({\mathbf{x}}))d\tau=e^{-\lambda
t({\mathbf{x}})}h_{\lambda}({\mathbf{y}})$

Since $\lambda>0$, we have

$\left|h_{\lambda}({\mathbf{x}})-\int_{0}^{t({\mathbf{x}})}e^{-\lambda\tau}{\mathbf{w}}_{\lambda}^{\top}{\mathbf{F}}_{n}({\mathbf{s}}_{\tau}({\mathbf{x}}))d\tau\right|\leq|h_{\lambda}({\mathbf{y}})|$

The proof then follows by using the contuinity property of
$h_{\lambda}({\mathbf{x}})$ as $\|{\mathbf{y}}\|\leq\epsilon$. ∎

Using the results of the above theorem, we can approximate the eigenfunction
for any point ${\mathbf{x}}\in\cal U$ containing the stable manifold as

$\displaystyle\hat{\phi}_{\lambda}({\mathbf{x}})={\mathbf{w}}^{\top}_{\lambda}{\mathbf{x}}+\int_{0}^{t({\mathbf{x}})}e^{-\lambda\tau}{\mathbf{w}}^{\top}_{\lambda}{\mathbf{F}}_{n}({\mathbf{s}}_{\tau}({\mathbf{x}}))d\tau.$
(18)

###  III-A Computational Framework for Stable Manifold

We first outline the algorithm to compute a stable manifold of the unstable
equilibrium point and then use these steps to obtain the simulation results.

Algorithm to determine stability boundary
1) Find all the equilibrium points and determine the equilibrium points,
${\mathbf{x}}_{i}^{\star}$, that are unstable on the stability boundary.
2) Let
${\mathbf{A}}_{i}=\frac{\partial{\mathbf{f}}}{\partial{\mathbf{x}}}({\mathbf{x}}_{i}^{\star})$
and ${\mathbf{w}}_{\lambda}$ be the left eigenvector corresponding to the
eigenvalue $\lambda$ with ${\rm Re}(\lambda)>0$. Construct,

${\cal
U}_{0}=\\{{\mathbf{x}}\in{\mathbf{X}}:{\mathbf{x}}^{\top}{\mathbf{P}}{\mathbf{x}}\leq\epsilon_{1}\\}$

for some positive matrix ${\mathbf{P}}$ whose $n-1$ eigenvectors span the
subspace $\\{{\mathbf{x}}:{\mathbf{w}}_{\lambda}^{\top}{\mathbf{x}}=0\\}$. We
will comment on the choice of $\epsilon_{1}$ and other parameters used in the
algorithm later.
3) Propagate the set ${\cal U}_{0}$ backward in time under the flow
${\mathbf{s}}_{-t}({\mathbf{x}})$ over time interval $[0,T]$ i.e.,

${\cal U}=\bigcup_{t=0}^{T}{\mathbf{s}}_{-t}({\cal U}_{0})$

4) Let $\\{{\mathbf{x}}_{k}\\}_{k=1}^{L}$ be the sample data points uniformly
distributed in the set ${\cal U}$. Use the path integral formula (18) to
compute the value of the eigenfunction corresponding to $\lambda$ at these
points.

5) Let
$\Psi({\mathbf{x}})=(\psi_{1}({\mathbf{x}}),\ldots,\psi_{N}({\mathbf{x}}))^{\top}$
be the finite choice of basis functions used for the approximation of
eigenfunction in the domain $\cal U$. Formulate the following least square
optimization problem to determine the approximation of the eigenfunction in
domain $\cal U$.

$\displaystyle\min_{{\mathbf{u}}\in{\mathbb{R}}^{N}}\|{\mathbf{G}}{\mathbf{u}}-{\mathbf{c}}\|$
(19)
$\displaystyle{\mathbf{c}}=(\hat{\phi}_{\lambda}({\mathbf{x}}_{1}),\ldots,\hat{\phi}_{\lambda}({\mathbf{x}}_{L}))^{\top},\;{\mathbf{G}}=(\Psi({\mathbf{x}}_{1}),\ldots,\Psi({\mathbf{x}}_{L}))^{\top}.$

The finite-dimensional approximation of the eigenfunction in the domain $\cal
U$ is given by

$\displaystyle\tilde{\phi}_{\lambda}({\mathbf{x}})=\Psi({\mathbf{x}})^{\top}{\mathbf{u}}^{\star},\;\;\;{\mathbf{u}}^{\star}={\mathbf{G}}^{\dagger}{\mathbf{c}}$
(20)

6) The approximation of the stability boundary can then be identified as the
zero-level curve of $\hat{\phi}_{\lambda}({\mathbf{x}})$. In particular, we
can determine the approximate stability boundary as $\hat{\cal
S}_{b}=\\{{\mathbf{x}}:|\tilde{\phi}_{\lambda}({\mathbf{x}})|\leq\gamma\\}$
for some small $\gamma>0$.

###### Remark 3

The choice of $\epsilon_{1}$ will depend upon the total time used for
simulation in backward time and the stable eigenvalue with the smallest real
part, $\bar{\lambda}$. The $\epsilon_{1}$ neighborhood around the unstable
equilibrium point will roughly grow in size as
$e^{-\bar{\lambda}T}\epsilon_{1}$. The choice of $\epsilon_{1}$ and $T$ are
connected to cover the stable manifold. The choice of $\gamma$ depends on how
coarsely we can resolve zero. With the value of eigenfunction computed at
different points in the state space, the $\gamma$ will be determined based on
the smallest positive and largest negative value of the eigenfunctions at the
sampled data.

We provide results on the convergence analysis for the finite-dimensional
approximation of the stable manifold. We start with the following assumption
on the basis functions.

###### Assumption 3

We assume that the basis function,
$\Psi({\mathbf{x}})=(\psi_{1}({\mathbf{x}}),\ldots,\psi_{N}({\mathbf{x}}))^{\top}$,
satisfies $\mu$ independence property,

$\displaystyle\mu\\{{\mathbf{x}}\in{\mathbb{R}}^{n}:{\mathbf{c}}^{\top}\Psi({\mathbf{x}})=0\\}=0,$
(21)

for all non-zero ${\mathbf{c}}\in{\mathbb{R}}^{N}$. The measure $\mu$ is
assumed to be equivalent to Lebesgue measure, $m$, on ${\mathbf{X}}$ i.e.,
$\mu(A)=0$ if and only if $m(A)=0$ for any measurable set
$A\subset{\mathbf{X}}$. The initial conditions
$\\{{\mathbf{x}}_{\ell}\\}_{\ell=1}^{L}$ are assumed to be independent
identically distributed (i.i.d) sampled from $\mu$.

Under Assumption 3, it can be shown that the matrix ${\mathbf{G}}$ in Eq. (19)
is invertible. However, due to space constraints we have to omit the proof.
Let ${\cal F}_{N}\subset{\cal L}_{2}({\mathbf{X}},\mu)$ a finite-dimensional
subspace spanned by
$\Psi({\mathbf{x}})=[\psi_{1}({\mathbf{x}}),\ldots,\psi_{N}({\mathbf{x}})]^{\top}$.
For $\phi\in{\cal L}_{2}({\mathbf{X}},\mu)$, we can define the projection of
$\phi$ on the closed subspace ${\cal F}_{N}\subset{\cal
L}_{2}({\mathbf{X}},\mu)$ as

$\displaystyle P_{N}^{\mu}\phi$ $\displaystyle=\arg\min_{f\in{\cal
F}_{N}}\|f-\phi\|_{{\cal L}_{2}({\mathbf{X}},\mu)}$
$\displaystyle=\Psi^{\top}({\mathbf{x}})\arg\min_{{\mathbf{c}}\in{\mathbb{R}}^{n}}\int_{\mathbf{X}}|{\mathbf{c}}^{\top}\Psi({\mathbf{x}})-\phi({\mathbf{x}})|^{2}d\mu({\mathbf{x}}).$

Let $\hat{\mu}_{L}$ be the empirical measure i.e.,
$\hat{\mu}_{L}=\frac{1}{L}\sum_{k}^{L}\delta_{{\mathbf{z}}_{k}}$, with
$\delta_{{\mathbf{z}}_{k}}$ the Dirac-delta measure. We then have

$\displaystyle P_{N}^{\hat{\mu}_{L}}\phi$ $\displaystyle=\arg\min_{f\in{\cal
F}_{N}}\|f-\phi\|_{{\cal L}_{2}({\mathbf{X}},\hat{\mu}_{L})}$
$\displaystyle=\Psi^{\top}({\mathbf{x}})\arg\min_{{\mathbf{c}}\in{\mathbb{R}}^{N}}\frac{1}{L}\sum_{i=1}^{M}\left({\mathbf{c}}^{\top}\Psi({\mathbf{x}}_{i})-\phi({\mathbf{x}}_{i})\right)^{2}.$

###### Lemma 1

We have

$\displaystyle\tilde{\phi}_{\lambda}({\mathbf{x}})\to
P_{N}^{\mu}\hat{\phi}_{\lambda}({\mathbf{x}})\;\;{\rm as}\;\;L\to\infty$ (22)

###### Proof:

From (19), (20), and (III-A), it follows that
$\tilde{\phi}_{\lambda}=P_{N}^{\hat{\mu}_{L}}\hat{\phi}_{\lambda}$. Since the
sample $\\{{\mathbf{x}}_{k}\\}$ are assumed to be drawn i.i.d. it follows from
law or large numbers that $P_{N}^{\hat{\mu}_{L}}\hat{\phi}_{\lambda}\to
P_{N}^{\mu}\hat{\phi}_{\lambda}$ as $L\to\infty$. ∎

##  IV Simulation Results

###  IV-A Bistable Toggle Switch

In the first example, we consider the dynamics of a genetic toggle switch. The
basic dynamics comprise two biological states, representing proteins
repressing each other to apply mutual negative feedback [16].

$\displaystyle\dot{x}_{1}=\frac{\alpha_{1}}{1+x_{2}^{\beta_{1}}}-\eta_{1}x_{1},\;\;\dot{x}_{2}=\frac{\alpha_{2}}{1+x_{1}^{\beta_{2}}}-\eta_{2}x_{2}.$
(23)

Here, for the device parameters
$\alpha_{1},\alpha_{2}=1,\beta_{1}=3.55,\beta_{2}=3.53,$ and
$\eta_{1},\eta_{2}=0.5,$ (23) have two stable equilibrium points and one
saddle point equilibrium at $(1,1)$.

Figure 1: Stability boundary obtained from phase portrait (blue) and estimated
stable manifold of unstable equilibrium point (red).

The linearization of the system at the origin has eigenvalues
$(0.3850,-1.3850)$. The left eigenvector corresponding to the unstable
eigenvalue is given by ${\mathbf{w}}_{\lambda}=(0.7061,-0.7081)^{\top}$. It is
clear from the phase portrait, Fig. 1 (blue curve), that the stability
boundary is almost linear for this example. Therefore, we approximate the
eigenfunction with linear only terms, i.e.,
${\mathbf{w}}_{\lambda}^{\top}{\mathbf{x}}$. The corresponding zero-level
curve of linear only eigenfunction is shown in Fig. 1 (red curve). From the
figure, it is evident that the zero-level curve of eigenfunction matches well
with the stability boundary.

###  IV-B 2-D Speed Control System

We consider a nonlinear speed control system as follows,

$\displaystyle\dot{x}_{1}=x_{2},\;\;\dot{x}_{2}=-K_{d}x_{2}-x_{1}-gx_{1}^{2}\left(\frac{x_{2}}{K_{d}}+x_{1}+1\right)$
(24)

Here, $K_{d}=1$ and $g=6$, the system has three equilibrium points,
$(-0.7886,0)$ and $(0,0)$ are stable equilibrium and $(-0.21135,0)$ is saddle
(type-1) equilibrium point. In this example, we are interested in finding the
stability boundary of $(0,0)$. The stable manifold of $(-0.21135,0)$,
characterized by a zero level set of the unstable eigenfunction of (24),
provides the stability boundary for $(0,0)$. The linearization of the system
around $(-0.21135,0)$ gives eigenvalues $(0.4309,-1.6990)$. The eigenfunction
corresponding to unstable eigenvalue $0.4309$ is the eigenfunction of
interest.

Figure 2: Speed control nonlinear system (a) Set ${\cal U}_{0}$ (blue) and,
set $\cal U$ obtained by backward reachable set (green) (b) Estimated stable
manifold (stability boundary) (red).

We consider a small region around $(-0.21135,0)$ representing ${\cal U}_{0}$
with $500$ points with $\epsilon_{1}=0.2$ as shown in Fig. 2(a) (blue) and
obtain the set $\cal U$ by backward propagation of (24) over time $[0,10]$ as
illustrated in Fig. 2(a) (green). We restrict the data points in ${\cal U}$ in
domain $[-1,1]^{2}$ to compute the eigenfunction values using (18) for $15850$
data points. In Fig. 2(b), we show the estimated stability boundary on top of
the phase portrait. The parameter $\gamma=5e^{-5}$. Clearly, the stability
boundary follows the phase portrait representing the efficacy of the proposed
approach.

###  IV-C Two Generator Infinite Bus Power System

The power system is a classic example of interest to engineering, where
determining stability boundary is crucial. For the power system operator, it
is of interest to know the critical clearing time (CCT), which is defined as
the maximum allowable time to clear the fault without causing instability. The
CCT is the time taken by the trajectory of the on-fault unstable power system
dynamics to cross the stability boundary of the post-fault stable power system
dynamics. The transient stability analysis of the power system to determine
the CCT is then effectively reduced to computing the stability boundary of the
post-fault stable power system dynamics. For this example, we consider a
three-generator system [17], with generator three as the reference bus as
follows.

$\displaystyle\dot{\delta}_{1}$
$\displaystyle=\omega_{1},\;\dot{\omega}_{1}=-\alpha_{1}\sin\delta_{1}-\beta_{1}\sin(\delta_{1}-\delta_{2})-D_{1}\omega_{1}$
(25) $\displaystyle\dot{\delta}_{2}$
$\displaystyle=\omega_{2},\;\dot{\omega}_{2}=-\alpha_{2}\sin\delta_{2}-\beta_{2}\sin(\delta_{2}-\delta_{1})-D_{2}\omega_{1}+P_{m}.$

Here, $\delta_{i}$ and $\omega_{i}$ are the generator rotor angle and angular
velocity, respectively, for $i^{th}$ generator, $\alpha_{i},\beta_{i}$ are the
ratio of generator internal voltage and line impedance with $\alpha_{1}=1$,
$\alpha_{2},\beta_{1},\beta_{2}=0.5$, $D_{1}=0.4$, and $D_{2}=0.5$. The system
has $6$ type-1 equilibrium points with $(\delta_{1},\delta_{2})$ values as
$(3.24,0.31)$, $(3.04,3.24)$, $(0.03,3.10)$, $(-3.03,0.31)$, $(-3.24,-3.03)$,
$(0.03,-3.17)$, and the post-fault stable equilibrium point is $(0.02,0.06)$
with $\omega_{1},\omega_{2}=0$. We consider ${\cal U}_{0}$ with $500$ data
points and $\epsilon_{1}=0.1$ and propagate this backward in time for $5$ sec.
to contain $\cal U$. We compute the stable manifold corresponding to each
unstable equilibrium point (UEP) as a zero-level set of the unstable
eigenfunction as shown in Fig. 3 (green) with a combination of linear and
trigonometric basis. The parameter $\gamma=1e^{-3}$.

Figure 3: Estimated stability boundary projected on $\delta_{1},\delta_{2}$
plane.

The joint level set of all the stable manifolds represents the stability
boundary shown in Fig. 3 (dotted red). We consider an on-fault case by
changing the parameters in (25) as $\alpha_{1},\alpha_{2}=0.01$,
$\beta{1}=0.05$, and $\beta_{2}=0.001$ shown in Fig. 3 (black). The CCT for
this case using time domain simulation is $43.8$sec, and the CCT obtained from
our approach, i.e., the time at which on-fault trajectory first crosses the
stable manifold, is $43.7$sec. Therefore, this example shows the efficacy of
our approach for practical applications of transient stability analysis in
power systems.

##  V Conclusions and Discussion

We have proposed a direct approach for the stability analysis of a dynamical
system using the Koopman spectrum. Compared to the existing approach based on
the Lyapunov function, our proposed approach is more direct as it relies on
the explicit computation of stability boundary using the zero-level curve of
Koopman eigenfunction. We also provide a path-integral formula for the
calculation of Koopman eigenfunction. One of the challenges with the proposed
approach is that it requires the knowledge of unstable equilibrium points on
the stability boundary.

## References

  * [1] H.-D. Chiang and L. F. Alberto, _Stability regions of nonlinear dynamical systems: theory, estimation, and applications_. Cambridge University Press, 2015.
  * [2] H. K. Khalil, _Nonlinear Systems_. New Jersey: Prentice Hall, 1996.
  * [3] P. A. Parrilo, “Structured semidefinite programs and semialgebraic geometry methods in robustness and optimization,” Ph.D. dissertation, California Institute of Technology, Pasadena, CA, 2000.
  * [4] P. Bevanda, M. Beier, S. Kerz, A. Lederer, S. Sosnowski, and S. Hirche, “Diffeomorphically learning stable Koopman operators,” _IEEE Control Systems Letters_ , vol. 6, pp. 3427–3432, 2022.
  * [5] H.-D. Chiang, M. Hirsch, and F. Wu, “Stability regions of nonlinear autonomous dynamical systems,” _IEEE Transactions on Automatic Control_ , vol. 33, no. 1, pp. 16–27, 1988.
  * [6] B. Hou, S. Bose, and U. Vaidya, “Sparse learning of kernel transfer operators,” _Asilomar Conference on Signals, Systems & Computers, 2021_.
  * [7] S. Klus, I. Schuster, and K. Muandet, “Eigendecompositions of transfer operators in reproducing kernel hilbert spaces,” _Journal of Nonlinear Science_ , vol. 30, no. 1, pp. 283–315, 2020.
  * [8] X. Ma, B. Huang, and U. Vaidya, “Optimal quadratic regulation of nonlinear system using koopman operator,” in _2019 American Control Conference (ACC)_ , 2019, pp. 4911–4916.
  * [9] M. Korda and I. Mezić, “On convergence of extended dynamic mode decomposition to the Koopman operator,” _Journal of Nonlinear Science_ , vol. 28, no. 2, pp. 687–710, 2018.
  * [10] ——, “Optimal construction of Koopman eigenfunctions for prediction and control,” _IEEE Transactions on Automatic Control_ , vol. 65, no. 12, pp. 5114–5129, 2020.
  * [11] A. Kumar, B. Umathe, U. Vaidya, and A. Kelkar, “Safe Operating Limits of Vehicle Dynamics Under Parameter Uncertainty Using Koopman Spectrum1,” _ASME Letters in Dynamic Systems and Control_ , vol. 3, no. 2, p. 021008, 10 2023. [Online]. Available: https://doi.org/10.1115/1.4063479
  * [12] A. Lasota and J. Yorke, “On the existence of invariant measure for piecewise monotonic transformations,” _Trans. A.M.S._ , vol. 186, pp. 481–488, 1973.
  * [13] I. Mezić, “Spectrum of the Koopman operator, spectral expansions in functional spaces, and state-space geometry,” _Journal of Nonlinear Science_ , vol. 30, no. 5, pp. 2091–2145, 2020.
  * [14] ——, “Koopman operator, geometry, and learning of dynamical systems,” _Not. Am. Math. Soc_ , vol. 68, no. 7, pp. 1087–1105, 2021.
  * [15] S. A. Deka, S. S. Narayanan, and U. Vaidya, “Path-integral formula for computing Koopman eigenfunctions,” _arXiv preprint arXiv:2307.06805, Accpted for publication in IEEE CDC_ , 2023.
  * [16] S. P. Nandanoori, S. Sinha, and E. Yeung, “Data-driven operator theoretic methods for phase space learning and analysis,” _Journal of Nonlinear Science_ , vol. 32, no. 6, p. 95, 2022.
  * [17] H.-D. Chiang, _Direct methods for stability analysis of electric power systems: theoretical foundation, BCU methodologies, and applications_. John Wiley & Sons, 2011.