Challenges in Dynamic Mode Decomposition

# Challenges in Dynamic Mode Decomposition

Ziyou Wu  Steven L. Brunton  & Shai Revzen

###### Contents

  1. 1 Abstract
  2. 2 Introduction
  3. 3 Background
    1. 3.1 Dynamic mode decomposition (DMD)
    2. 3.2 Koopman theory
  4. 4 Methods
    1. 4.1 A “simple” linear system
    2. 4.2 Nonlinear observations of the system
    3. 4.3 Dataset preparation
  5. 5 Results
    1. 5.1 Fixed spectrum with different system matrix normality
    2. 5.2 Different spectrum with orthogonal eigenspace
  6. 6 Discussion and Conclusion
  7. 7 Acknowledgements
  8. 8 References
  9. A Appendix
    1. A.1 Density of linearizable systems
    2. A.2 Diagonal linear systems
    3. A.3 All 3x3 cases are covered
    4. A.4 9D linear system with non-resonant modes
    5. A.5 Other DMD algorithms
    6. A.6 5D system

##  1 Abstract

Dynamic Mode Decomposition (DMD) is a powerful tool for extracting spatial and
temporal patterns from multi-dimensional time series, and it has been used
successfully in a wide range of fields, including fluid mechanics, robotics,
and neuroscience. Two of the main challenges remaining in DMD research are
noise sensitivity and issues related to Krylov space closure when modeling
nonlinear systems. Here, we investigate the combination of noise and
nonlinearity in a controlled setting, by studying a class of systems with
linear latent dynamics which are observed via multinomial observables. Our
numerical models include system and measurement noise. We explore the
influences of dataset metrics, the spectrum of the latent dynamics, the
normality of the system matrix, and the geometry of the dynamics. Our results
show that even for these very mildly nonlinear conditions, DMD methods often
fail to recover the spectrum and can have poor predictive ability. Our work is
motivated by our experience modeling multilegged robot data, where we have
encountered great difficulty in reconstructing time series for oscillatory
systems with slow transients, which decay only slightly faster than a period.

##  2 Introduction

The dynamic mode decomposition (DMD) is a dimensionality reduction and
modeling approach that was developed in the fluids community [1, 2]. DMD is a
modal decomposition technique [3, 4] that decomposes a matrix of high-
dimensional time series data into a set of spatial coherent structures that
exhibit the same linear dynamics in time (e.g., oscillations, exponential
growth/decay). In this way, the DMD may be thought of as an algorithmic
descendent of the principal components analysis (PCA) in space with the
Fourier transform in time [5, 6, 7, 8]. Since its introduction, DMD has been
rigorously connected to nonlinear dynamical systems [9], as an approximation
of the infinite dimensional Koopman operator [10, 11]. The Koopman operator,
introduced in 1931, provides an alternative operator-theoretic perspective for
dynamical systems, and it has experienced a resurgence because of its utility
for data-driven analysis [12, 13, 14, 15, 16, 17].

Since its introduction, DMD has been applied to a wide range of fields beyond
fluid mechanics [2, 18], including robotics [19], disease modeling [20],
neuroscience [21], finance [22], and plasma physics [23, 24], to name only a
few [7]. One of the key reasons for the broad success of DMD is its
formulation in terms of linear algebra, statistical regression, and
optimization, making it relatable to the recent work in machine learning [25]
and highly extensible. Within a short time, DMD has seen dozens of algorithmic
extensions, including for control [26, 27], for large and streaming data sets
[28, 29, 30, 31] incorporating ideas of randomized linear algebra [32, 33],
for model reduction [34, 35], for limited measurements [36, 37, 38, 39] and
time delays [40, 41], for multiresolution analysis [42], for selecting a
sparse subset of dominant modes [43], and for Bayesian formulations [44],
among others [45, 46]. The PyDMD [47] software package has emerged as an open
source tool for several of these extensions.

Two of the main avenues of DMD research have centered around noise robustness
and closure issues related to modeling nonlinear systems. It was recognized
early that DMD is quite sensitive to noise [48], and several approaches have
been proposed to address this sensitivity, including a forward-backward
averaging [49], total least-squares regression [50], variable projection [51],
variational approaches [45], and robust principal component analysis [52].
However, noise sensitivity is still a leading challenge, and is one that we
directly investigate in this paper. Similarly, soon after the initial
connection between DMD and Koopman theory [9], there were several follow-up
studies addressing how and when a linear regression model can capture
essential features of a nonlinear system. The extended DMD algorithm was
developed to augment the DMD state with nonlinear observables to enrich the
regression [53, 54, 55], and theoretical convergence results prove that in the
infinite data limit, these models converge to the projection of the Koopman
operator on this observable subspace [56]. However, often the observable
subspace is not closed under action of the Koopman operator, so that it does
not form an invariant subspace, resulting in poorly predicted dynamics [57].

In our own work we have sometimes encountered great difficulty in
reconstructing time series with DMD. A particularly difficult class seems to
be the reconstruction of intermediate transients – transient phenomena that
decay faster than the slowest mode, but not much faster. Here we investigate
the possible sources of these difficulties by studying special class of
systems for which both the latent linearization and the linearizing
observables are completely known – the class systems with linear latent
dynamics which are observed via multinomial observables:

$\displaystyle\mathbf{x}_{k+1}=\mathbf{A}\mathbf{x}_{k};$
$\displaystyle\mathbf{y}=P(\mathbf{x})$ (1)

where $P(\mathbf{x})$ is a (typically low order) multinomial. However, since
we are considering models of physical systems, we will assume both system
noise, and measurement noise:

$\displaystyle\mathbf{x}_{k+1}=\mathbf{A}\mathbf{x}_{k}+(\text{system
noise});$ $\displaystyle\mathbf{y}=P(\mathbf{x})+(\text{measurement noise}).$
(2)

In this paper, we explore the combination of noise and nonlinearity in this
controlled setting, where it is possible to isolate and analyze each effect
with incrementally increasing severity. Specifically, we explore a
particularly benign class of such systems, having 3-dimensional latent
dynamics and observables that are monotone polynomials in individual
coordinates. We demonstrate that even for these very mildly nonlinear systems,
observed over a subspace closed under the Koopman operator, DMD methods often
fail to recover the spectrum and can have poor predictive value. We also
explore the influences of dataset metrics, the spectrum of the latent
dynamics, the normality of the system matrix, and the geometry of the
dynamics. Based on results from this self-contained framework, we give
generalizable recommendations regarding dataset properties for DMD analysis.
All code used to generate these results are available as open source at: .

##  3 Background

Here we will introduce the dynamic mode decomposition (DMD), which is a recent
technique for linear system identification. DMD is closely related to Koopman
operator theory, which provides a representation of a dynamical system in
terms of the evolution its “observables”. We will also introduce Koopman
theory, as it relates to performing DMD on data that is augmented with
nonlinear observables.

###  3.1 Dynamic mode decomposition (DMD)

Here we present the _exact DMD_ formulation of Tu et al. [6], along with
several leading variants to handle noisy data [50, 49, 51]. Exact DMD is
modified from the original algorithm of Schmid [1], and both algorithms seek a
best-fit linear operator $\mathbf{A}\in\mathbb{R}^{n\times n}$ that evolves a
high-dimensional spatial measurement vector $\mathbf{x}\in\mathbb{R}^{n}$
forward in time:

$\displaystyle\mathbf{x}_{k+1}=\mathbf{A}\mathbf{x}_{k}.$ (3)

More specifically, DMD seeks the leading eigenvalues and eigenvectors of the
matrix $\mathbf{A}$. The DMD eigenvectors $\boldsymbol{\phi}$, or _modes_ ,
have the dimension of the original data $\mathbf{x}$, and correspond to
spatial structures that behave coherently in time according to the
corresponding eigenvalue $\lambda$.

The DMD algorithm begins with time-series data
$\\{\mathbf{x}_{1},\mathbf{x}_{2},\cdots,\mathbf{x}_{m}\\}$, which is arranged
into two data matrices:

$\displaystyle\mathbf{X}=\begin{bmatrix}\mathbf{x}_{1}&\mathbf{x}_{2}&\cdots&\mathbf{x}_{m-1}\end{bmatrix}$
$\displaystyle\mathbf{X}^{\prime}=\begin{bmatrix}\mathbf{x}_{2}&\mathbf{x}_{3}&\cdots&\mathbf{x}_{m}\end{bmatrix}.$
(4)

The best-fit linear operator $\mathbf{A}$ may be solved for in the data
equation

$\displaystyle\mathbf{X}^{\prime}\approx\mathbf{A}\mathbf{X}$ (5)

by a least-squares minimization:

$\displaystyle\mathbf{A}=\text{argmin}_{\mathbf{A}}\|\mathbf{X}^{\prime}-\mathbf{A}\mathbf{X}\|_{F}=\mathbf{X}^{\prime}\mathbf{X}^{\dagger}$
(6)

where $\|\cdot\|_{F}$ is the Frobenius norm and $\mathbf{X}^{\dagger}$ is the
pseudoinverse of $\mathbf{X}$, computed via the singular value decomposition
(SVD):

$\displaystyle\mathbf{X}=\mathbf{U}\boldsymbol{\Sigma}\mathbf{V}^{\mathsf{T}}\quad\Longrightarrow\quad\mathbf{X}^{\dagger}=\mathbf{V}\boldsymbol{\Sigma}^{-1}\mathbf{U}^{\mathsf{T}}.$
(7)

In practice, the SVD of $\mathbf{X}$ is typically truncated so that only
$r\ll\min({m,n})$ modes are retained:
$\mathbf{X}\approx\mathbf{U}_{r}\boldsymbol{\Sigma}_{r}\mathbf{V}^{\mathsf{T}}_{r}$.

DMD was introduced for high-dimensional systems, such as fluid dynamics, where
the state dimension $n$ may be in the millions. In this case, the matrix
$\mathbf{A}$ may be intractably large, and instead of constructing
$\mathbf{A}$ directly, it’s low-dimensional representation
$\tilde{\mathbf{A}}$ is analyzed in the $r$-dimensional subspace given by the
columns of $\mathbf{U}_{r}$:

$\displaystyle\tilde{\mathbf{A}}=\mathbf{U}_{r}^{\mathsf{T}}\mathbf{A}\mathbf{U}_{r}=\mathbf{U}_{r}^{\mathsf{T}}\mathbf{X}^{\prime}\mathbf{V}_{r}\boldsymbol{\Sigma}_{r}^{-1}.$
(8)

The eigenvalues of $\tilde{\mathbf{A}}$ are the same as the eigenvalues of
$\mathbf{A}$, and are given by the diagonal elements of the matrix
$\boldsymbol{\Lambda}$

$\displaystyle\tilde{\mathbf{A}}\mathbf{W}=\mathbf{W}\boldsymbol{\Lambda}.$
(9)

The corresponding DMD modes, or eigenvectors of the high-dimensional system,
are computed from the projected eigenvectors $\mathbf{W}$ as

$\displaystyle\boldsymbol{\Phi}=\mathbf{X}^{\prime}\mathbf{V}_{r}^{\mathsf{T}}\boldsymbol{\Sigma}_{r}^{-1}\mathbf{W}.$
(10)

Tu et al. [6] showed that the columns of $\boldsymbol{\Phi}$ are exact
eigenvectors of the high-dimensional matrix $\mathbf{A}$.

Once the DMD eigenvalues and eigenvectors are computed, it is possible to
approximate the time series data through the DMD expansion, which is closely
related to the Koopman mode decomposition:

$\displaystyle\mathbf{x}_{k}\approx\sum_{j=1}^{r}\boldsymbol{\phi}_{j}\lambda^{k-1}b_{j}=\boldsymbol{\Phi\Lambda}^{k-1}\mathbf{b}.$
(11)

The vector $\mathbf{b}$ contains the amplitudes of each DMD mode, and is often
computed as

$\displaystyle\mathbf{b}=\boldsymbol{\Phi}^{\dagger}\mathbf{x}_{1}.$ (12)

Written in matrix form, the DMD expansion becomes

$\displaystyle\mathbf{X}\approx\underbrace{\begin{bmatrix}\vline&&\vline\\\
\boldsymbol{\phi}_{1}&\cdots&\boldsymbol{\phi}_{r}\\\
\vline&&\vline\end{bmatrix}}_{\boldsymbol{\Phi}}\underbrace{\begin{bmatrix}b_{1}&&\\\
&\ddots&\\\
&&b_{r}\end{bmatrix}}_{\text{diag}(\mathbf{b})}\underbrace{\begin{bmatrix}1&\cdots&\lambda_{1}^{m-2}\\\
\vdots&\ddots&\vdots\\\
1&\cdots&\lambda_{r}^{m-2}\end{bmatrix}}_{\mathbf{T}(\boldsymbol{\Lambda})}.$
(13)

The estimation of $\mathbf{b}$ above is based on the first snapshot of data
only. It is possible to improve this estimate with the following optimization
over all snapshots:

$\displaystyle\text{argmin}_{\mathbf{b}}\|\mathbf{X}-\boldsymbol{\Phi}\text{diag}(\mathbf{b})\mathbf{T}(\boldsymbol{\Lambda})\|_{F}.$
(14)

Jovanovic et al. [43] further add a sparsity promoting penalty to make the
vector $\mathbf{b}$ as sparse as possible.

Askham and Kutz [51] introduced the _optimized_ DMD algorithm that
simultaneously optimizes over the modes and eigenvalues using variable
projection. In this approach, they combine the modes and amplitudes
$\boldsymbol{\Phi}_{\mathbf{b}}=\boldsymbol{\Phi}\text{diag}(\mathbf{b})$ and
solve the following optimization

$\displaystyle\min_{\boldsymbol{\Lambda},\boldsymbol{\Phi}_{\mathbf{b}}}\|\mathbf{X}-\boldsymbol{\Phi}_{\mathbf{b}}\mathbf{T}(\boldsymbol{\Lambda})\|_{F}.$
(15)

The optimized DMD approach has been proven to be quite robust to noisy data
and data that is not perfectly approximated by a linear system.

The sensitivity of DMD to noise was recognized early by Bagheri [48], and
several DMD variants have been proposed to improve the robustness of the
algorithm. Two leading approaches include the forward-backward (FB) DMD [49]
and the total least-squares (TLS) DMD [50]. Both algorithms correct for sensor
noise and process noise. Forward-backward DMD averages the standard DMD
operator $\mathbf{A}$ and the inverse of a reverse-time DMD operator
$\mathbf{A}^{b}$ computed by switching the order of $\mathbf{X}$ and
$\mathbf{X}^{\prime}$ in the algorithm. The averaged operator

$\displaystyle\mathbf{A}_{\text{FB}}=\left(\mathbf{A}\left(\mathbf{A}^{b}\right)^{-1}\right)^{1/2}$
(16)

cancels out the bias introduced from noise. Similarly, total least-squares DMD
makes note of the fact that if DMD is viewed as a regression from
$\mathbf{x}_{k}$ to $\mathbf{x}_{k+1}$, then noise will affect both the
dependent and independent variables. Thus, replacing the standard SVD-based
least-squares regression with a total least-squares algorithm will improve the
robustness to noise and remove bias. In the TLS algorithm, the data
$\mathbf{X}$ and $\mathbf{X}^{\prime}$ are stacked into a larger matrix that
is used to compute the low-dimensional subspace

$\displaystyle\mathbf{Z}=\begin{bmatrix}\mathbf{X}\\\
\mathbf{X}^{\prime}\end{bmatrix}=\mathbf{U_{Z}}\boldsymbol{\Sigma}_{\mathbf{Z}}\mathbf{V_{Z}}^{\mathsf{T}}.$
(17)

The $\mathbf{U_{Z}}$ matrix is partitioned into four sub-matrices

$\displaystyle\mathbf{U_{Z}}=\begin{bmatrix}\mathbf{U}_{\mathbf{Z},a}&\mathbf{U}_{\mathbf{Z},b}\\\
\mathbf{U}_{\mathbf{Z},c}&\mathbf{U}_{\mathbf{Z},d}\end{bmatrix}$ (18)

and the debiased estimate of the DMD operator $\mathbf{A}$ is given by

$\displaystyle\mathbf{A}_{\text{TLS}}=\mathbf{U}_{\mathbf{Z},c}\mathbf{U}_{\mathbf{Z},a}^{\dagger}.$
(19)

In this work, we will compare the exact DMD algorithm, along with the forward-
backward, total-least squares, and optimized DMD variants.

###  3.2 Koopman theory

DMD was quickly connected to nonlinear dynamical systems through Koopman
operator theory [9]. Here we provide a brief introduction of some important
concepts that will be relevant when performing DMD on a state augmented with
nonlinear observables.

Let $\mathbf{F}:\mathbb{X}\times\mathbb{R}\to\mathbb{X}$ be the flow of a
continuous- or discrete- time dynamical system, denoted
$\mathbf{x}(t)=\mathbf{F}^{t}(\mathbf{x}(0))$ where
$\mathbf{x}(0),\mathbf{x}(t)\in\mathbb{X}$. By “flow” we mean that
$\mathbf{F}$ is a representation of the group $(\mathbb{R},+)$, i.e. $\forall
t,s\in\mathbb{R}:\,\mathbf{F}^{t}\circ\mathbf{F}^{s}=\mathbf{F}^{t+s}$ and
$\mathbf{F}^{0}$ is the identity map on $\mathbb{X}$.

An “observable” is a map $\mathbf{g}:\mathbb{X}\to\mathbb{V}$ into some space
of observation values $\mathbb{V}$. We will usually assume $\mathbb{V}$ is of
the form $\mathbb{C}^{d}$: some vector space over the field of complex
numbers.

The “Koopman Operators” (a.k.a. “composition operators”) $\mathcal{K}^{t}$ act
on observables (e.g. $\mathbf{g}$) by:

$\forall
t\in\mathbb{R},\,\mathbf{x}\in\mathbb{X}:\,(\mathcal{K}^{t}\mathbf{g})(\mathbf{x})=\mathbf{g}(\mathbf{F}^{t}(\mathbf{x})).$
(20)

Our interest lies in finding observables $\mathbf{g}$ such that
$\mathbf{g}\circ\mathbf{F}^{t}$ is a linear system in $\mathbb{V}$, i.e. there
exists a matrix $\mathbf{M}\in\text{Aut}(\mathbb{V})$ such that
$\exp(t\mathbf{M})\cdot\mathbf{g}=\mathbf{g}\circ\mathbf{F}^{t}=\mathcal{K}^{t}\mathbf{g}$.
Such $\mathbf{g}$ spans an “eigenspace” of $\mathcal{K}$ with a finite-
dimensional evolution operator $\mathbf{M}$. In the case where $\mathbf{M}$ is
a scalar, this coincides with the familiar definition of eigenvectors of a
linear operator.

Note that for any two eigenfunctions $g_{{}_{\mathbf{A}}},g_{{}_{\mathbf{B}}}$
with commuting (matrix-valued) eigenvalues $\mathbf{A}$ and $\mathbf{B}$, the
function $g(\mathbf{x}):=g_{{}_{\mathbf{A}}}(\mathbf{x})\cdot
g_{{}_{\mathbf{B}}}(\mathbf{x})$ satisfies

$g\circ\mathbf{F}^{t}=(g_{{}_{\mathbf{A}}}\circ\mathbf{F}^{t})\cdot(g_{{}_{\mathbf{B}}}\circ\mathbf{F}^{t})=\exp(t\mathbf{A})\cdot
g_{{}_{\mathbf{A}}}\cdot\exp(t\mathbf{B})\cdot
g_{{}_{\mathbf{B}}}=\exp(t(\mathbf{A}+\mathbf{B}))\cdot(g_{{}_{\mathbf{A}}}\cdot
g_{{}_{\mathbf{B}}})=\exp(t(\mathbf{A}+\mathbf{B}))\cdot g,$ (21)

and thus the pointwise product of eigenfunctions whose eigenvalues commute is
itself an eigenfunction whose eigenvalues are the sum of the constituent
eigenvalues. This also implies (trivially) that non-negative integer powers of
eigenfunctions are eigenfunctions.

Together these observations imply that any monomial of eigenfunctions
$g_{1}^{\alpha_{1}}\cdots g_{m}^{\alpha_{m}}$ whose respective matrix
eigenvalues were $\mathbf{M}_{k}$ for $g_{k}$, and for whom all
$\mathbf{M}_{k}$ commute, is itself an eigenfunction which has the matrix
eigenvalue $\alpha_{1}\mathbf{M}_{1}+\ldots+\alpha_{m}\mathbf{M}_{m}$.

It is also possible to compute DMD on an augmented state that includes
nonlinear functions of the measured state $\mathbf{x}$, in a procedure called
_extended DMD_ [53, 54, 55]. By including nonlinear functions, it is possible
to approximate the projection of the Koopman operator onto a larger space of
functions, where it may be possible to better approximate the relevant
eigenfunctions. It was recently shown that extended DMD is equivalent to the
earlier _variational approach of conformation dynamics_ (VAC) [58, 59, 60],
introduced by Noé and Nüske to simulate molecular dynamics with a broad
separation of timescales. Further connections between eDMD and VAC and between
DMD and the time lagged independent component analysis (TICA) are explored in
a recent review [61]. This work will explore extended DMD where the state is
augmented with monomial functions.

##  4 Methods

Generic stable linear systems can exhibit a great many behaviors. To explore
several of these, including some of the interactions between different
eigenvalues and eigenspace geometries, we constructed a 3-dimensional linear
system.

###  4.1 A “simple” linear system

We created a stable real linear system, with states
$\mathbf{x}\in\mathbb{R}^{3}$, and considered specifically the case of a
spectrum that consists of a (stable) complex-conjugate pair and a single
(stable) real eigenvalue. Such a system has two invariant subspaces – a two-
dimensional subspace in which the complex-conjugate pair generate a spiral
sink, and a one-dimensional subspace in which the real eigenvalue induces
exponential convergence. We chose, without loss of generality, to make the two
dimensional invariant space the XY plane with the major and minor axes of the
ellipsoidal spiral to be aligned with the X and Y axes. Thus, the most general
3D system matrix $\mathbf{A}$, with a spectrum consisting of one complex-
conjugate pair $\alpha\pm i\beta$ and one real eigenvalue $\lambda$, can be
constructed as follows. We took $\mathbf{A}$ to be constructed from an
eigenvector matrix $\mathbf{Q}$, a diagonal matrix $\mathbf{S}$ to rescale the
major vs the minor axis in the XY plane, and a block diagonal eigenvalue
matrix $\mathbf{\Lambda}$,

$\mathbf{A}:=\mathbf{QS\Lambda}\mathbf{S}^{-1}\mathbf{Q}^{-1}.$ (22)

We parameterized $\mathbf{Q}$, $\mathbf{S}$, $\mathbf{\Lambda}$ by

$\mathbf{Q}:=\left[\begin{array}[]{ccc}1&0&\sin(\theta)\cos(\phi)\\\
0&1&\sin(\theta)\sin(\phi)\\\
0&0&\cos(\theta)\end{array}\right],\mathbf{S}:=\left[\begin{array}[]{ccc}s&0&0\\\
0&1&0\\\
0&0&1\end{array}\right],\mathbf{\Lambda}:=\left[\begin{array}[]{ccc}\alpha&\beta&0\\\
-\beta&\alpha&0\\\ 0&0&\lambda\end{array}\right].$ (23)

The parameters $\phi$, $\theta$, and $s$ exhaust all possible 3D linear
dynamics with the spectrum of $\mathbf{\Lambda}$ up to an orthogonal
similarity transform, as proved in A.3.

###  4.2 Nonlinear observations of the system

We took as our nonlinear observations of the system the output

$\displaystyle\mathbf{y}:=[x_{1}+0.1(x_{1}^{2}+x_{2}x_{3}),~{}x_{2}+0.1(x_{2}^{2}+x_{1}x_{3}),~{}x_{3}+0.1(x_{3}^{2}+x_{1}x_{2})]^{\mathsf{T}}.$
(24)

For the range of $\mathbf{x}$ values we considered (initial conditions on the
unit sphere), $\mathbf{y}$ coordinates were strictly monotone in the
corresponding $\mathbf{x}$ coordinates, and observations corresponded to
unique states. As figure 1 illustrates, the trajectories of $\mathbf{y}$ and
of $\mathbf{x}$ are so similar that for a casual observer they would seem
nearly indistinguishable. However, as we will show in 5, these nonlinearities
can have profound impact on our ability to model the system using DMD.

Figure 1: Example trajectories of our 3D system eqn. 2 for various parameters
of eqn. 23 (colors). Eigenvalues were $0.9(\sqrt{3}/2\pm 0.5i)$ and $0.6$
throughout, and we used the same initial condition. We plotted trajectories
without noise added for both the latent $\mathbf{x}$ (dots), and the mildly
nonlinear observable $\mathbf{y}$ (eqn. 24; stars), in 3D perspective (left)
and in an XY plane projection (right). In the 3D plot we indicated the
direction of the $0.6$ eigenspace controlled by $\theta$ and $\phi$ and given
by the last column of $Q$ in eqn. 23 (dashed line). We showed the effects of
$s$ in the $XY$ plane using ellipsoids (dashed line) whose major and minor
axes have ratio $s^{2}$. The trajectories demonstrate that our nonlinear
observations are quite similar to the latent linear dynamics, and that
trajectories vary their geometry with the parameters of eqn. 23 even while the
spectrum remains the same.

Because $\mathbf{A}$ had distinct eigenvalues, it was diagonalizable over
$\mathbb{C}$. Consider the diagonal matrix $\mathbf{\tilde{\Lambda}}$ with 9
diagonal elements comprising the spectrum of $\mathbf{A}$ and all 6 pairwise
sums of eigenvalues of $\mathbf{A}$. Every monomial of the form
$x_{1}^{k_{1}}x_{2}^{k_{2}}x_{3}^{k_{3}}$ with $1\leq k_{1}+k_{2}+k_{3}\leq 2$
evolves exponentially over time as one of the diagonal elements of
$\mathbf{\tilde{\Lambda}}$. We conclude that the dynamics of $\mathbf{y}$
could just as well have been written as
$\mathbf{z}_{n+1}:=\mathbf{P}\mathbf{\tilde{\Lambda}}\mathbf{P}^{-1}\mathbf{z}_{n}$,
$\mathbf{y}:=\mathbf{W}\mathbf{z}$ for matrices $\mathbf{P}$ and $\mathbf{W}$.

This is of utmost importance, since it implies that _this nonlinearly observed
3D linear system is exactly reproduced in a 9D linear system observed
linearly_. More generally, all monomial observables of order $m$ or less can
be represented using eigenspaces whose eigenvalues can be obtained as a sum of
no more than $m$ of the original eigenvalues; see A.2 for details. We may
therefore rest assured that a DMD method that allowed for 9 modes could in
fact predict the values of $\mathbf{y}$ to without any truncation errors due
to the nonlinearity of the observations.

###  4.3 Dataset preparation

Our numerical experiments used the 3D systems, with Gaussian system noise and
Gaussian measurement noise both set to zero mean and a standard deviation of
$0.05$. The initial conditions were chosen uniformly distributed on the unit
sphere.

We explored all combinations of $s\in\\{0.1,0.5,1.0\\}$,
$\phi\in\\{0,0.79,\pi/2\\}$,
$\theta\in\\{0,1.1,1.2,1.3,1.4,1.5,1.52,1.53,1.560\\}$ with a fixed eigenvalue
spectrum $0.5(\sqrt{3}/2\pm 0.5i)$ and $0.8$. We also explored different
spectra with complex conjugate eigenvalue pair $0.5(\sqrt{3}/2\pm 0.5i)$ and a
real eigenvalue varying from 0.01 to 1, with orthogonal eigenvector space
$s=1,\theta=0$. Furthermore, we considered datasets consisting of $100$ points
each with: (i) $N=50$ initial conditions for trajectories of length $L=2$;
(ii) $N=10$ initial conditions for trajectories of length $L=10$; (iii) $N=2$
initial conditions for trajectories of length $L=50$. Because of the presence
of system noise, all of these datasets meet the condition of persistence of
excitation [62].

In the following, we present contour plots showing the distribution of
eigenvalue estimates. Wherever we show such plots, they smoothed with a
Gaussian kernel of $\sigma=0.01$. We first sorted all eigenvalues by their
imaginary parts, and then registered each to its nearest true eigenvalue bin.
To ensure the spectrum symmetry, we discarded estimates with more real
eigenvalues than the true spectrum. If there were any extra complex
eigenvalues pairs, we registered two copies of the average value. Furthermore,
we added data to these contour plots in batches of $300$ numerical
experiments. This number was sufficient so that the KL-divergence in $P_{1}$
and $P_{2}$, the probability distributions before and after a new estimated
eigenvalue $\hat{\lambda}$ was added, was smaller than 0.001 for all
experiments, i.e.
$\sum_{\hat{\lambda}}P_{1}(\hat{\lambda})\log\left(\frac{P_{1}(\hat{\lambda})}{P_{2}(\hat{\lambda})}\right)<0.001$.

##  5 Results

In this section, we present the eigenvalue estimation by exact DMD on a 3D
linear system with mild nonlinear measurements and observation noise. We
examine the effects of system normality, spectrum, mild nonlinearity in the
observables, and varying number of initial conditions through their effect on
the recovered spectrum and the prediction of the dynamics. We found that
mildly nonlinear observations and slight non-normality caused large eigenvalue
estimation errors. For a system with a well-posed system matrix, having a
dataset with more initial conditions and shorter trajectories can
significantly improve the prediction. With a slightly ill-conditioned system
matrix, a moderate trajectory length improves the spectrum recovery.

###  5.1 Fixed spectrum with different system matrix normality

We fixed the spectrum of the system matrix to be $0.5(\sqrt{3}/2\pm 0.5i)$ and
$0.8$, and varied its $\phi$, $s$, and $\theta$ as described in eqn. 22. The
linear states evolved with Gaussian system noise $\sigma=0.05$, and they were
observed via monomials with Gaussian observation noise $\sigma=0.05$. In
figure 2 and figure 4, we show the DMD eigenvalue density plot under five
different normality settings, using first order and second order monomial
observables. In the first column, the system matrix $\mathbf{A}$ is well-
conditioned, with orthogonal eigenvectors and a centrosymmetric spiral sink.
In the second column, the eigenvectors remain orthogonal, and we stretch the
spiral sink by varying the parameter $s$, so that the resulting trajectory
converges following an ellipsoidal spiral with a minor-to-major axes ratio of
$s^{2}$. The eigenvalue density contour is stretched vertically when $s$
decreases, which yields a larger variance in the complex frequency estimation.
In the third and fifth column, we varied the angle between the 1-dimensional
stable subspace corresponding to the real eigenvalue and the 2-dimensional
subspace corresponding to the spiral sink, with the spiral sink remaining
centrosymmetric. The fourth column exhibits both non-orthogonal eigenvectors
and a non-centrosymmetric spiral sink. When eigenvectors are non-orthogonal, a
complicated multi-modal error structure shows up in the
$L=10,50,\theta=1.40,1.56$ cases.

To thoroughly explore the effect of system matrix normality on spectrum
recovery and dynamics prediction, in figure 3 and figure 5 we plotted the
standard deviation of each estimated DMD eigenvalue versus the condition
number. The complex eigenvalues always come out as conjugate pairs from exact
DMD, so we give one plot per pair. The discarded trials plot (right lower
corner in figure 3 and figure 5) gives the percentage of DMD estimates
resulting in a structurally incorrect estimation, i.e. where the number of
real eigenvalues is larger than the ground truth system, because it breaks the
symmetry pattern of the complex conjugate pairs. When the condition number is
below 100, the eigenvalue estimation from a dataset with more initial
conditions and short trajectories (squares) has lower std, and less discarded
trials. When the condition number is larger than 100, the cases with moderate
trajectory length (circles) performs better.

The estimation error caused by system matrix non-normality can be amplified
through higher order observables. In the first two columns of figure 2, we
find that the eigenvalue contours are close to the true eigenvalues. The
second order observables in figure 4 amplified the effect of slight changes in
the condition number of the system matrix. Only the slowest modes (red and
purple contours) were identifiable. Furthermore, in figure 5, the standard
deviation of the eigenvalue estimates increases to 0.25 for most eigenvalues
with a modest condition number of 4.

Figure 2: Eigenvalue density plot with first order linear observables.
Different rows have different trajectory lengths and number of trajectories,
while holding the same total data points. The system $\mathbf{A}$ matrix is
constructed by $\phi=0$ and $\theta,s$ specified on top of each column.
$\kappa(\mathbf{A})$ is the condition number of ground truth system matrix,
and $\kappa(\tilde{\mathbf{A}})$ is the mean of estimated condition number
found by exact DMD. Red crosses are ground truth eigenvalues. Figure 3: The
standard deviation of each estimated eigenvalue versus system $\mathbf{A}$
matrix condition number plot. The colors of different eigenvalues match the
color of contours in figure 2. Figure 4: Eigenvalue density plot with monomial
observables up to 2nd order. Different rows have different trajectory lengths
and number of trajectories, while holding the same total data points. The
system $\mathbf{A}$ matrix is constructed by $\phi=0$ and $\theta,s$ specified
on top of each column. $\kappa(\mathbf{A})$ is the condition number of ground
truth system matrix, and $\kappa(\tilde{\mathbf{A}})$ is the mean of estimated
condition number found by exact DMD. Red crosses are ground truth eigenvalues.
Figure 5: The standard deviation of each estimated eigenvalue versus system
$\mathbf{A}$ matrix condition number plot. The colors of different eigenvalues
match the color of contours in figure 4.

###  5.2 Different spectrum with orthogonal eigenspace

In this section, we kept an orthogonal eigenvector space with $s=1,\theta=0$
and a complex conjugate eigenvalue pair $0.5(\sqrt{3}/2\pm 0.5i)$. We varied
the real eigenvalue from 0 to 1. The linear states evolved with Gaussian
system noise $\sigma=0.05$, and were observed via monomials having Gaussian
observation noise $\sigma=0.05$. In figure 6, we show the DMD eigenvalue
density plot under three different spectra using data prepared as described in
4.3. In figure 7 and figure 8, we show the standard deviation of the estimated
DMD eigenvalues versus varying the eigenvalue $\lambda$.

More initial conditions with shorter trajectories uniformly outperforms less
initial conditions with longer trajectories, resulting in more accurate and
precise estimations. In the case of second order monomial observables with a
trajectory length of L=2, the standard deviation of the estimated eigenvalues
increases when $\lambda$ moves closer to the conjugate pair. This suggests a
dataset with short trajectories might be more sensitive to a small spectral
distance. Nonetheless, this case still has lower std comparing to less initial
conditions with longer trajectories. When $\lambda$ moves from 0.2 to 1.0, the
percentage of discarded trials for L=2 decreases and stays at 0, for L=10 this
percentage remains about the same, and for L=50 it increases to 50%.

Figure 6: Eigenvalue density plot with linear observables (left 9 plots) and
observables up to second order monomials (right 9 plots). Different rows have
different trajectory length, and different columns have different real
eigenvalue. Red crosses are ground truth eigenvalues. Figure 7: The standard
deviation versus real eigenvalue. The colors of different eigenvalues match
the color of contours in the left 9 plots of figure 6.  Figure 8: The standard
deviation versus real eigenvalue. The colors of different eigenvalues match
the color of contours in the right 9 plots of figure 6.

##  6 Discussion and Conclusion

We have shown that under conditions of mild system and observation noise, and
mild nonlinearity, DMD can encounter significant difficulties in recovering a
linear model with the correct spectrum and an accurate prediction of the
dynamics, even when an exact linear model is known to exist. To study these
effects, we constructed a self-contained 3D linear system where it is possible
to independently control the eigenvalue spectrum, the normality of the system
matrix, the nonlinearity of the observations, the system and observation
noise, and the number and length of the trajectories used as training data.

In this study, we find that DMD is quite sensitive to non-normality (see figs.
2, 3, 4, 5, and supplemental figs. 12, 13), and to the existence of higher
order nonlinear terms (see figs. 4, 6 and supplemental figs. 11, 13). As to
the question of which data-set metrics govern this sensitivity, in some cases
shorter data with more initial conditions performs better (as in figs. 6, 7,
8), and in some data fewer, longer time series produce better results (see 2
4). The challenge of non-normality has implications for the identification of
fluid systems, as highly sheared flows are generally characterized by non-
normal linearized dynamics, which was an original motivation for DMD [1].

We first fixed the spectrum of our linear system, varying the normality of
system matrix and the geometry of trajectories. The results suggest that an
ill-condition system matrix and mildly-nonlinear observables cause DMD to fail
to recover the full spectrum; this phenomenon is especially evident in figure
4. Furthermore, DMD estimation can be very sensitive to noise when the angle
between eigenspaces is small. For the case of linear observables, when we
decrease the angle between invariant subspaces, a multi-modal error structure
appears in figure 2. The complex eigenvalue pair and the real eigenvalue
appear to have switched real parts. The cause of this error structure is not
understood and remains an interesting avenue of future work.

We then explored systems with different eigenvalue spectra, where we changed
the distance between the eigenvalues. For short trajectories ($L=2$), closer
eigenvalues result in higher standard deviation in the predicted spectrum,
using second order observables. Generally, with enough initial conditions, DMD
yields descent prediction of the spectra for all well-conditioned systems
explored here.

DMD under second order observations generally does not provide good
estimations for transient oscillatory modes. To ensure this is not the effect
of spectrum resonance caused by first and second order monomial observations,
we compute the same contour plot on a 9D non-resonant linear system with
similar spectrum in figure 9. The result is similar to the right 9 plots of
figure 6. Further, we show that using the latest denoising DMD algorithms in
§A.5, the difficulty in spectrum recovery often persists, although each method
performs well for certain cases.

The 3D linear system analyzed here is a special case of a larger class of
systems designed specifically to test the sensitivity of DMD to the challenges
listed above. The class of diagonal linear systems with multinomial
observables provides an algebraically convenient yet sufficiently general
class of systems to use to study these problems. As the theorems in appendix
sections A.1 and A.2 show, diagonal linear systems with multinomial
observables are universal in the sense of being able to approximate any
smooth, stable dynamical system (details in Theorem 1). Furthermore, these
systems have a finite dimensional linear representation which is exact – a
linear observation of a linear dynamical system (see Theorem 2) – and
computable in close form. Thus, the class of diagonal linear systems with
multinomial observables is a useful space to test any technique that
identifies linear models, including the many variants of DMD.

There are several future directions that are motivated by this work. It will
be important to explore how these challenges scale to higher dimensional
systems. There is also the potential to design sampling strategies to improve
the conditioning of the DMD procedure. In addition to understanding and
characterizing these challenges on numerical examples, it will be interesting
to apply these careful studies to low-dimensional mechanical systems.

##  7 Acknowledgements

The authors acknowledge funding support from the Army Research Office
(W911NF-17-1-0306).

##  8 References

## References

  * [1] P. J. Schmid, “Dynamic mode decomposition of numerical and experimental data,” _Journal of Fluid Mechanics_ , vol. 656, pp. 5–28, Aug. 2010.
  * [2] P. Schmid, L. Li, M. Juniper, and O. Pust, “Applications of the dynamic mode decomposition,” _Theoretical and Computational Fluid Dynamics_ , vol. 25, no. 1-4, pp. 249–259, 2011.
  * [3] K. Taira, S. L. Brunton, S. Dawson, C. W. Rowley, T. Colonius, B. J. McKeon, O. T. Schmidt, S. Gordeyev, V. Theofilis, and L. S. Ukeiley, “Modal analysis of fluid flows: An overview,” _AIAA Journal_ , vol. 55, no. 12, pp. 4013–4041, 2017.
  * [4] K. Taira, M. S. Hemati, S. L. Brunton, Y. Sun, K. Duraisamy, S. Bagheri, S. Dawson, and C.-A. Yeh, “Modal analysis of fluid flows: Applications and outlook,” _AIAA Journal_ , vol. 58, no. 3, pp. 998–1022, 2020.
  * [5] K. K. Chen, J. H. Tu, and C. W. Rowley, “Variants of dynamic mode decomposition: Boundary condition, Koopman, and Fourier analyses,” _Journal of Nonlinear Science_ , vol. 22, no. 6, pp. 887–915, 2012.
  * [6] J. H. Tu, C. W. Rowley, D. M. Luchtenburg, S. L. Brunton, and J. N. Kutz, “On dynamic mode decomposition: theory and applications,” _Journal of Computational Dynamics_ , vol. 1, no. 2, pp. 391–421, 2014.
  * [7] J. N. Kutz, S. L. Brunton, B. W. Brunton, and J. L. Proctor, _Dynamic Mode Decomposition: Data-Driven Modeling of Complex Systems_. SIAM, 2016.
  * [8] S. L. Brunton and J. N. Kutz, _Data-Driven Science and Engineering: Machine Learning, Dynamical Systems, and Control_. Cambridge University Press, 2019.
  * [9] C. W. Rowley, I. Mezić, S. Bagheri, P. Schlatter, and D. Henningson, “Spectral analysis of nonlinear flows,” _J. Fluid Mech._ , vol. 645, pp. 115–127, 2009.
  * [10] B. O. Koopman, “Hamiltonian systems and transformation in Hilbert space,” _Proceedings of the National Academy of Sciences_ , vol. 17, no. 5, pp. 315–318, 1931.
  * [11] B. Koopman and J. v. Neumann, “Dynamical systems of continuous spectra,” _Proceedings of the National Academy of Sciences of the United States of America_ , vol. 18, no. 3, p. 255, 1932.
  * [12] I. Mezić and A. Banaszuk, “Comparison of systems with complex behavior,” _Physica D: Nonlinear Phenomena_ , vol. 197, no. 1, pp. 101–133, 2004.
  * [13] I. Mezić, “Spectral properties of dynamical systems, model reduction and decompositions,” _Nonlinear Dynamics_ , vol. 41, no. 1-3, pp. 309–325, 2005\.
  * [14] M. Budišić, R. Mohr, and I. Mezić, “Applied Koopmanism a),” _Chaos: An Interdisciplinary Journal of Nonlinear Science_ , vol. 22, no. 4, p. 047510, 2012.
  * [15] I. Mezic, “Analysis of fluid flows via spectral properties of the Koopman operator,” _Annual Review of Fluid Mechanics_ , vol. 45, pp. 357–378, 2013\.
  * [16] S. Bagheri, “Koopman-mode decomposition of the cylinder wake,” _Journal of Fluid Mechanics_ , vol. 726, pp. 596–623, 2013.
  * [17] S. L. Brunton, M. Budišić, E. Kaiser, and J. N. Kutz, “Modern Koopman theory for dynamical systems,” _arXiv preprint arXiv:2102.12086_ , 2021.
  * [18] A. Towne, O. T. Schmidt, and T. Colonius, “Spectral proper orthogonal decomposition and its relationship to dynamic mode decomposition and resolvent analysis,” _Journal of Fluid Mechanics_ , vol. 847, pp. 821–867, 2018.
  * [19] E. Berger, M. Sastuba, D. Vogt, B. Jung, and H. B. Amor, “Estimation of perturbations in robotic behavior using dynamic mode decomposition,” _Journal of Advanced Robotics_ , vol. 29, no. 5, pp. 331–343, 2015.
  * [20] J. L. Proctor and P. A. Eckhoff, “Discovering dynamic patterns from infectious disease data using dynamic mode decomposition,” _International health_ , vol. 7, no. 2, pp. 139–145, 2015.
  * [21] B. W. Brunton, L. A. Johnson, J. G. Ojemann, and J. N. Kutz, “Extracting spatial-temporal coherent patterns in large-scale neural recordings using dynamic mode decomposition,” _Journal of Neuroscience Methods_ , vol. 258, pp. 1–15, 2016.
  * [22] J. Mann and J. N. Kutz, “Dynamic mode decomposition for financial trading strategies,” _Quantitative Finance_ , pp. 1–13, 2016.
  * [23] R. Taylor, J. N. Kutz, K. Morgan, and B. Nelson, “Dynamic mode decomposition for plasma diagnostics and validation,” _Review of Scientific Instruments_ , vol. 89, no. 053501, 2018.
  * [24] A. A. Kaptanoglu, K. D. Morgan, C. J. Hansen, and S. L. Brunton, “Characterizing magnetized plasmas with dynamic mode decomposition,” _Physics of Plasmas_ , vol. 27, no. 3, p. 032108, 2020.
  * [25] S. L. Brunton, B. R. Noack, and P. Koumoutsakos, “Machine learning for fluid mechanics,” _Annual Review of Fluid Mechanics_ , vol. 52, pp. 477–508, 2020\.
  * [26] J. L. Proctor, S. L. Brunton, and J. N. Kutz, “Dynamic mode decomposition with control,” _SIAM Journal on Applied Dynamical Systems_ , vol. 15, no. 1, pp. 142–161, 2016.
  * [27] E. Kaiser, J. N. Kutz, and S. L. Brunton, “Sparse identification of nonlinear dynamics for model predictive control in the low-data limit,” _Proceedings of the Royal Society of London A_ , vol. 474, no. 2219, 2018\.
  * [28] M. S. Hemati, M. O. Williams, and C. W. Rowley, “Dynamic mode decomposition for large and streaming datasets,” _Physics of Fluids (1994-present)_ , vol. 26, no. 11, p. 111701, 2014.
  * [29] T. Sayadi and P. J. Schmid, “Parallel data-driven decomposition algorithm for large-scale datasets: with application to transitional boundary layers,” _Theoretical and Computational Fluid Dynamics_ , pp. 1–14, 2016.
  * [30] S. D. Pendergrass, J. N. Kutz, and S. L. Brunton, “Streaming GPU singular value and dynamic mode decompositions,” _arXiv preprint arXiv:1612.07875_ , 2016.
  * [31] H. Zhang, C. W. Rowley, E. A. Deem, and L. N. Cattafesta, “Online dynamic mode decomposition for time-varying systems,” _arXiv preprint arXiv:1707.02876_ , 2017.
  * [32] D. Bistrian and I. Navon, “Randomized dynamic mode decomposition for non-intrusive reduced order modelling,” _International Journal for Numerical Methods in Engineering_ , 2016.
  * [33] N. B. Erichson, L. Mathelin, J. N. Kutz, and S. L. Brunton, “Randomized dynamic mode decomposition,” _SIAM Journal on Applied Dynamical Systems_ , vol. 18, no. 4, pp. 1867–1891, 2019.
  * [34] G. Tissot, L. Cordier, N. Benard, and B. R. Noack, “Model reduction using dynamic mode decomposition,” _Comptes Rendus Mécanique_ , 2014.
  * [35] B. R. Noack, W. Stankiewicz, M. Morzynski, and P. J. Schmid, “Recursive dynamic mode decomposition of a transient cylinder wake,” _Journal of Fluid Mechanics_ , vol. 809, pp. 843–872, 2016.
  * [36] S. L. Brunton, J. L. Proctor, J. H. Tu, and J. N. Kutz, “Compressed sensing and dynamic mode decomposition,” _Journal of Computational Dynamics_ , vol. 2, no. 2, pp. 165–191, 2015.
  * [37] F. Gueniat, L. Mathelin, and L. Pastur, “A dynamic mode decomposition approach for large and arbitrarily sampled systems,” _Physics of Fluids_ , vol. 27, no. 2, p. 025113, 2015.
  * [38] Z. Bai, E. Kaiser, J. L. Proctor, J. N. Kutz, and S. L. Brunton, “Dynamic mode decomposition for compressive system identification,” _AIAA Journal_ , vol. 58, no. 2, pp. 561–574, 2017.
  * [39] B. Kramer, P. Grover, P. Boufounos, M. Benosman, and S. Nabi, “Sparse sensing and DMD based identification of flow regimes and bifurcations in complex flows,” _SIAM Journal on Applied Dynamical Systems_ , vol. 16, no. 2, pp. 1164–1196, 2017.
  * [40] S. L. Brunton, B. W. Brunton, J. L. Proctor, E. Kaiser, and J. N. Kutz, “Chaos as an intermittently forced linear system,” _Nature Communications_ , vol. 8, no. 19, pp. 1–9, 2017.
  * [41] H. Arbabi and I. Mezic, “Ergodic theory, dynamic mode decomposition, and computation of spectral properties of the Koopman operator,” _SIAM Journal on Applied Dynamical Systems_ , vol. 16, no. 4, pp. 2096–2126, 2017.
  * [42] J. N. Kutz, X. Fu, and S. L. Brunton, “Multi-resolution dynamic mode decomposition,” _SIAM Journal on Applied Dynamical Systems_ , vol. 15, no. 2, pp. 713–735, 2016.
  * [43] M. R. Jovanović, P. J. Schmid, and J. W. Nichols, “Sparsity-promoting dynamic mode decomposition,” _Physics of Fluids_ , vol. 26, no. 2, p. 024103, 2014.
  * [44] N. Takeishi, Y. Kawahara, Y. Tabei, and T. Yairi, “Bayesian dynamic mode decomposition,” _Twenty-Sixth International Joint Conference on Artificial Intelligence_ , 2017.
  * [45] O. Azencot, W. Yin, and A. Bertozzi, “Consistent dynamic mode decomposition,” _SIAM Journal on Applied Dynamical Systems_ , vol. 18, no. 3, pp. 1565–1585, 2019.
  * [46] A. G. Nair, B. Strom, B. W. Brunton, and S. L. Brunton, “Phase-consistent dynamic mode decomposition from multiple overlapping spatial domains,” _Phys. Rev. Fluids_ , vol. 5, p. 074702, 2020.
  * [47] N. Demo, M. Tezzele, and G. Rozza, “PyDMD: Python dynamic mode decomposition,” _The Journal of Open Source Software_ , vol. 3, no. 22, p. 530, 2018.
  * [48] S. Bagheri, “Effects of weak noise on oscillating flows: Linking quality factor, Floquet modes, and Koopman spectrum,” _Physics of Fluids_ , vol. 26, no. 9, p. 094104, 2014.
  * [49] S. T. Dawson, M. S. Hemati, M. O. Williams, and C. W. Rowley, “Characterizing and correcting for the effect of sensor noise in the dynamic mode decomposition,” _Experiments in Fluids_ , vol. 57, no. 3, pp. 1–19, 2016\.
  * [50] M. S. Hemati, C. W. Rowley, E. A. Deem, and L. N. Cattafesta, “De-biasing the dynamic mode decomposition for applied Koopman spectral analysis,” _Theoretical and Computational Fluid Dynamics_ , vol. 31, no. 4, pp. 349–368, 2017.
  * [51] T. Askham and J. N. Kutz, “Variable projection methods for an optimized dynamic mode decomposition,” _SIAM Journal on Applied Dynamical Systems_ , vol. 17, no. 1, pp. 380–416, 2018.
  * [52] I. Scherl, B. Strom, J. K. Shang, O. Williams, B. L. Polagye, and S. L. Brunton, “Robust principal component analysis for modal decompositionn of corrupt fluid flows,” _Physical Review Fluids_ , vol. 5, no. 054401, 2020\.
  * [53] M. O. Williams, C. W. Rowley, I. Mezić, and I. G. Kevrekidis, “Data fusion via intrinsic dynamic variables: An application of data-driven Koopman spectral analysis,” _EPL (Europhysics Letters)_ , vol. 109, no. 4, p. 40007, 2015.
  * [54] M. O. Williams, C. W. Rowley, and I. G. Kevrekidis, “A kernel approach to data-driven Koopman spectral analysis,” _Journal of Computational Dynamics_ , vol. 2, no. 2, pp. 247–265, 2015.
  * [55] Q. Li, F. Dietrich, E. M. Bollt, and I. G. Kevrekidis, “Extended dynamic mode decomposition with dictionary learning: A data-driven adaptive spectral decomposition of the koopman operator,” _Chaos: An Interdisciplinary Journal of Nonlinear Science_ , vol. 27, no. 10, p. 103111, 2017.
  * [56] M. Korda and I. Mezić, “On convergence of extended dynamic mode decomposition to the Koopman operator,” _arXiv preprint arXiv:1703.04680_ , 2017.
  * [57] S. L. Brunton, B. W. Brunton, J. L. Proctor, and J. N. Kutz, “Koopman invariant subspaces and finite linear representations of nonlinear dynamical systems for control,” _PLoS ONE_ , vol. 11, no. 2, p. e0150171, 2016.
  * [58] F. Noé and F. Nuske, “A variational approach to modeling slow processes in stochastic dynamical systems,” _Multiscale Modeling & Simulation_, vol. 11, no. 2, pp. 635–655, 2013.
  * [59] F. Nüske, B. G. Keller, G. Pérez-Hernández, A. S. Mey, and F. Noé, “Variational approach to molecular kinetics,” _Journal of chemical theory and computation_ , vol. 10, no. 4, pp. 1739–1752, 2014.
  * [60] F. Nüske, R. Schneider, F. Vitalini, and F. Noé, “Variational tensor approach for approximating the rare-event kinetics of macromolecular systems,” _J. Chem. Phys._ , vol. 144, no. 5, p. 054105, 2016.
  * [61] S. Klus, F. Nüske, P. Koltai, H. Wu, I. Kevrekidis, C. Schütte, and F. Noé, “Data-driven model reduction and transfer operator approximation,” _Journal of Nonlinear Science_ , 2018.
  * [62] M. Green and J. B. Moore, “Persistence of excitation in linear systems,” in _1985 American Control Conference_ , 1985, pp. 412–417.
  * [63] M. D. Kvalheim, D. Hong, and S. Revzen, “Generic properties of koopman eigenfunctions for stable fixed points and periodic orbits,” 2020. [Online]. Available: https://arxiv.org/abs/2010.04008

##  Appendix A Appendix

At first, it might appear that the class of systems described by multinomial
observations of diagonal linear dynamics is rather restrictive. As we show in
Section A.1, the converse is true – these systems are quite general.
Furthermore, we show in Section A.2 that such systems can be viewed instead of
as a multinomial observation of a non-resonant linear system, as a linear
observation of an extended linear systems whose eigenvalues fall on a lattice.
Thus a system that identifies linear observations of linear dynamics is
actually quite general.

###  A.1 Density of linearizable systems

Consider a physical model in the form of a finite dimensional ordinary
differential equation $\dot{\mathbf{x}}=F_{0}(\mathbf{x})$, with
$x\in\mathbb{X}$, and the “flow” $\mathbf{F}_{0}$ associated with
$F_{0}(\cdot)$ as defined in §3.2.

The vast majority of control theory as used in practical engineering
applications, and a large part of the mathematical modeling done in the
physical sciences, deals with the behavior of such systems around equilibria,
i.e. points of the form $F_{0}(\mathbf{x}_{0})=0$. Without loss of generality,
we shall assume that the equilibrium of interest is at $0$, i.e. $F_{0}(0)=0$.
Furthermore such models are usually of stable systems, i.e. have a
neighborhood $U$ in which for all $\mathbf{x}\in U$,
$\lim_{t\to\infty}\mathbf{F}_{0}^{t}(\mathbf{x})=0$ for the flow
$\mathbf{F}_{0}$. We will restrict our discussion to a compact neighborhood
$U$ of $0$ for which the dynamics are stable.

Borrowing from the definitions in [63], we define
$\mathfrak{X}^{\infty}_{\text{fix}}(U)$ as the set of all
$\mathcal{C}^{\infty}$ smooth vector fields with a stable fixed point at $0$
with stability basin $U$. Within this set there exists
$\mathcal{G}_{\text{fix}}(U)\in\mathfrak{X}^{\infty}_{\text{fix}}(U)$ which
are those vector fields whose Jacobian at $0$ has distinct eigenvalues.

######  Theorem 1.

For any choice of $\mathcal{C}^{\infty}$ smooth vector field
$F_{0}\in\mathfrak{X}^{\infty}_{\text{fix}}(U)$ which converges to
$0\in\mathbb{R}^{n}$, any choice of finite observation time $T>0$, and any
desired accuracy $\delta>0$, there exist:

  1. 1.

a linear system with a diagonal system matrix $\Lambda\in\mathbf{M}_{n\times
n}(\mathbb{C})$

  2. 2.

a multinomial observable of that system $P:\mathbb{C}^{n}\to U$

  3. 3.

an invertible mapping of initial conditions $S:U\to\mathbb{C}^{n}$

such that the trajectory of the flow $\mathbf{F}_{0}^{t}$ starting at
$\mathbf{x}$ is $\delta$-close to the $P$ observations of the trajectory of
the linear system $\dot{\mathbf{y}}=\Lambda\mathbf{y}$ starting at
$S(\mathbf{x})$.

###### Proof.

Using [63, Theorem 5], we learn that $\mathcal{G}_{\text{fix}}(U)$ is (among
other things) dense in $\mathfrak{X}^{\infty}_{\text{fix}}(U)$ with respect to
the infinity norm, i.e. for every $\varepsilon>0$ and
$F_{0}\in\mathfrak{X}^{\infty}_{\text{fix}}(U)$ there is a
$F_{\varepsilon}\in\mathcal{G}_{\text{fix}}(U)$ such that
$\|F_{\varepsilon}-F_{0}\|_{\infty}<\varepsilon$. From [63, Theorem 16], we
obtain that the $n$ principal Koopman eigenfunctions of $F_{\varepsilon}$
exist, they are uniquely defined, and together they provide a
$\mathcal{C}^{\infty}$ diffeomorphism from $\Psi:U\to\mathbb{C}^{n}$, that
conjugates the flow $\mathbf{F}_{\varepsilon}$:
$\mathbf{F}^{t}_{\varepsilon}=\Psi^{-1}\circ\exp(t\mathrm{\mathbf{D}}F_{\varepsilon}(0))\cdot\Psi$
for all $t$. From the Stone-Weierstrass theorem it follows that there exists a
multinomial approximation $\Psi^{-1}_{\varepsilon}$ to $\Psi^{-1}$ such that
within $\Psi(U)$,
$\|\Psi^{-1}_{\varepsilon}-\Psi^{-1}\|_{\infty}<\varepsilon$. This gives
$\|\mathbf{F}^{t}_{\varepsilon}-\Psi^{-1}_{\varepsilon}\circ\exp(t\mathrm{\mathbf{D}}F_{\varepsilon}(0))\cdot\Psi\|_{\infty}<\varepsilon$.
Note that because it is $\mathcal{C}^{\infty}$ on a compact set $U$, $F_{0}$
is Lipschitz; assume its Lipschitz constant is $L$. Consider the difference
between observing trajectories of the true flow $\mathbf{F}_{0}^{t}$ and
trajectories of the multinomial observations $\Psi^{-1}_{\varepsilon}$ of the
diagonal linear system $\exp(t\mathrm{\mathbf{D}}F_{\varepsilon}(0))$ over a
finite time horizon $T$. There is a well known result bounding the difference
between flows using the infinity norm of the difference between the vector
fields:

$\displaystyle\|\mathbf{F}^{t}_{0}-\mathbf{F}^{t}_{\varepsilon}\|_{\infty}\leq\frac{\exp(tL)-1}{L}\|F_{0}-F_{\varepsilon}\|_{\infty}<\varepsilon
C$

with $C:=(\exp(TL)-1)/L$. Using the triangle inequality this yields

$\displaystyle\|\mathbf{F}^{t}_{0}-\Psi^{-1}_{\varepsilon}\circ\exp(t\mathrm{\mathbf{D}}F_{\varepsilon}(0))\cdot\Psi\|_{\infty}<(C+1)\varepsilon$
(25)

We complete the proof by choosing $\varepsilon<\delta/(C+1)$; diagonalizing
$\mathrm{\mathbf{D}}F_{\varepsilon}(0)=V^{-1}\Lambda V$ with its eigenvector
matrix $V$ (which is possible because
$F_{\varepsilon}\in\mathcal{G}_{\text{fix}}(U)$); defining
$P:=\Psi^{-1}_{\varepsilon}\circ V^{-1}$; and defining $S:=V\cdot\Psi$. ∎

###  A.2 Diagonal linear systems

We now further strengthen the result of Theorem 1, to show what is in some
ways a proof that finite Koopman linearizations can approximate all
$\mathfrak{X}^{\infty}_{\text{fix}}$ systems.

######  Theorem 2.

Consider the linear dynamical system $\dot{\mathbf{x}}=\Lambda\mathbf{x}$ with
diagonal system matrix $\Lambda$ on $\mathrm{\mathbf{X}}=\mathbb{C}^{n}$. Let
$h:\mathrm{\mathbf{X}}\to\mathbb{C}$ be an observable which is a multinomial
(multivariate polynomial) of order $m$ in the components of
$\mathbf{x}\in\mathrm{\mathbf{X}}$. Then there exists an $N\in\mathbb{N}$,
$N\geq n$ and a linear system with a diagonal system matrix
$\tilde{\Lambda}\in\mathbf{M}_{N\times N}(\mathbb{C})$, a linear observable of
that system $H:\mathbb{C}^{N}\to\mathbb{C}$, and a mapping of initial
conditions $\Psi:U\to\mathbb{C}^{N}$, such that for all $t$ and $\mathbf{x}$:

$h(\exp(t\Lambda)\cdot\mathbf{x})=H\cdot\exp(t\tilde{\Lambda})\cdot\Psi(\mathbf{x}).$

###### Proof.

The proof is primarily a matter of notation.

Let $\text{diag}(\Lambda)=\lambda_{1},\ldots,\lambda_{n}$:

$\exp(t\Lambda)\cdot\mathbf{x}=\exp(t\Lambda)\cdot(x_{1},\ldots,x_{n})=(\mathrm{e}^{t\lambda_{1}}x_{1},\ldots,\mathrm{e}^{t\lambda_{n}}x_{n}).$
(26)

An $h$ can be written as

$h(\mathbf{x})=\sum_{\alpha_{1}+\ldots+\alpha_{n}\leq
m}H_{(\alpha_{1},\ldots,\alpha_{n})}x_{1}^{\alpha_{1}}\cdots
x_{n}^{\alpha_{n}}=\sum_{\|\alpha\|_{1}\leq m}H_{\alpha}\mathbf{x}^{\alpha}$
(27)

where the latter uses multi-index notation for $\alpha$. The number of non-
zero $H_{\alpha}$ coefficients is the $N$ needed.

However, this same class of systems can be viewed in a different light. Define
$\mathrm{\mathbf{X}}^{[m]}$ the space given by monomial features of order $m$
or less from $\mathrm{\mathbf{X}}$, i.e. the image of $\mathrm{\mathbf{X}}$
under the map

$\Psi(\mathbf{x}):=\bigoplus_{\|\alpha\|_{1}\leq m}\mathbf{x}^{\alpha}$ (28)

where $\|\cdot\|_{1}$ is the 1-norm. We denote by $\Psi_{\alpha}(\mathbf{x})$
the $\alpha$ component of this result. This definition allows us to view $\sum
H_{\alpha}\mathbf{x}^{\alpha}$ as an inner product, and $H$ as a linear
functional $L(\mathrm{\mathbf{X}}^{[m]},\mathbb{C})$, giving

$h(\mathbf{x})=\sum_{\|\alpha\|_{1}\leq
m}H_{\alpha}\mathbf{x}^{\alpha}=H\cdot\Psi(\mathbf{x}).$ (29)

We have “lifted” $\mathrm{\mathbf{X}}$ into a feature space
$\mathrm{\mathbf{X}}^{[m]}$ in which $h$ admits a linear representation.

Given that
$\exp(t\Lambda)\cdot\mathbf{x}=(\mathrm{e}^{t\lambda_{1}}x_{1},\ldots,\mathrm{e}^{t\lambda_{n}}x_{n})$,
note that for any multi-index $\alpha$ we can define
$\lambda_{[{\alpha}]}:=\sum_{i=1}^{n}\lambda_{i}\alpha_{i}$ and have

$\Psi_{\alpha}(\exp(t\Lambda)\cdot\mathbf{x})=(\exp(t\Lambda)\cdot\mathbf{x})^{\alpha}=\mathrm{e}^{t\lambda_{1}\alpha_{1}}x_{1}^{\alpha_{1}}\cdots\mathrm{e}^{t\lambda_{n}\alpha_{n}}x_{n}^{\alpha_{n}}=\mathrm{e}^{t\lambda_{[{\alpha}]}}x^{\alpha}=\mathrm{e}^{t\lambda_{[{\alpha}]}}\Psi_{\alpha}(\mathbf{x})$
(30)

In other words, if we define a diagonal matrix $\tilde{\Lambda}$ whose
diagonal elements are $\lambda_{[{\alpha}]}$ for the same collection of multi-
indices $\alpha$ used to create $\Psi$,
$\Psi(\exp(t\Lambda)\cdot\mathbf{x})=\exp(t\tilde{\Lambda})\Psi(\mathbf{x})$.
∎

Using Theorem 2, we have shown that the approximation obtained from Theorem 1
in terms of non-linear observations of a linear system can in fact be viewed
as a linear observation of an extended linear system. This extended linear
system has a very special structure – its eigenvalues fall on the lattice
$\lambda_{[{\alpha}]}$ generated from the original system’s eigenvalues
$\lambda$. It is also initialized only on the (dynamically invariant)
algebraic variety which is the image of $\Psi$,
$\mathrm{\mathbf{X}}^{[m]}\subset\mathbb{C}^{N}$. Nonetheless, the result
tells us that all systems in $\mathfrak{X}^{\infty}_{\text{fix}}(U)$ can be
approximated to any chosen accuracy $\delta>0$ and for any time horizon $T>0$
by a linear observation of a linear system.

###  A.3 All 3x3 cases are covered

######  Lemma A.1.

Any $3\times 3$ real matrix $\mathbf{M}$ with one complex-conjugate pair
$\alpha\pm i\beta$ and one real eigenvalue $\lambda$ can be parametrized by
eqn. 22 and eqn. 23 up to an orthogonal similarity transform $\mathbf{U}$.
i.e.

$\mathbf{M}=\mathbf{UAU}^{\mathsf{T}}=\mathbf{UQS\Lambda
S}^{-1}\mathbf{Q}^{-1}\mathbf{U}^{\mathsf{T}}.$ (31)

###### Proof.

Any $3\times 3$ real matrix $\mathbf{M}$ with this spectrum has three distinct
eigenvalues, and therefore, is block diagonalizable. Let $\mathbf{v}_{1}$ be a
(complex) eigenvector associated with $\alpha+i\beta$ and $\mathbf{v}_{2}$ be
a (real) eigenvector associated with $\lambda$. The eigendecomposition of $M$
gives $\mathbf{M}=\mathbf{P\Lambda P^{-1}}$, where
$\mathbf{P}=[\text{Re}(\mathbf{v}_{1}),\text{Im}(\mathbf{v}_{1}),\mathbf{v}_{2}]$.
W.l.o.g., we can transform $\mathbf{M}$ by a rotation matrix $\mathbf{U}_{1}$,
so that the transformed matrix
$\mathbf{M^{\prime}}=\mathbf{U_{1}MU_{1}^{T}}=\mathbf{P^{\prime}\Lambda
P^{\prime-1}}$ has complex eigenvector pair in $xy$-plane, i.e.

$\mathbf{P}^{\prime}=\left[\begin{array}[]{ccc}\text{Re}(\mathrm{v}_{1,1}^{\prime})&\text{Im}(\mathrm{v}_{1,1}^{\prime})&\vdots\\\
\text{Re}(\mathrm{v}_{1,2}^{\prime})&\text{Im}(\mathrm{v}_{1,2}^{\prime})&\mathbf{v}_{2}^{\prime}\\\
0&0&\vdots\end{array}\right].$

The proof is the same as showing there always exists an orthogonal similarity
transform $\mathbf{U}_{2}$,
s.t.
$\mathbf{M^{\prime}=U}_{2}\mathbf{QS\Lambda}\mathbf{S}^{-1}\mathbf{Q}^{-1}\mathbf{U}^{-1}_{2}$.
Let $\mathbf{U}_{2}$ be a rotation matrix in $xy$-plane:

$\mathbf{U}_{2}=\left[\begin{array}[]{ccc}\cos(\omega)&-\sin(\omega)&0\\\
\sin(\omega)&\cos(\omega)&0\\\
0&0&1\end{array}\right],\mathbf{U}_{2}\mathbf{QS}=\left[\begin{array}[]{ccc}s\cos(\omega)&-\sin(\omega)&\sin(\theta)\cos(\omega+\phi)\\\
s\sin(\omega)&\cos(\omega)&\sin(\theta)\sin(\omega+\phi)\\\
0&0&\cos(\theta)\end{array}\right].$

The third column of $\mathbf{U}_{2}\mathbf{QS}$ can match any eigenvector
$\mathbf{v}_{2}^{\prime}$ by picking angle $\theta$ and $\phi$. The following
Mathematica code demonstrates that we can always pick $\omega$ and $s$ to
match the upper-left $2\times 2$ block. In this code, $k_{1}=\frac{ly+xz}{lx-
yz},k_{2}=\frac{l^{2}+z^{2}}{lx-yz}$ and $\left[\begin{array}[]{cc}x&y\\\
z&l\end{array}\right]=\left[\begin{array}[]{cc}\text{Re}(\mathrm{v}_{1,1}^{\prime})&\text{Im}(\mathrm{v}_{1,1}^{\prime})\\\
\text{Re}(\mathrm{v}_{1,2}^{\prime})&\text{Im}(\mathrm{v}_{1,2}^{\prime})\end{array}\right]$.
Thus a solution satisfying $0<s\leq 1$ and $0\leq\omega<2\pi$ always exists.

⬇ Lam = {{a,b},{-b,a}};  U ={{Cos[w],-Sin[w]},{Sin[w],Cos[w]}};  T = {{1,1},
{-I,I}};  S ={{s,0},{0,1}};  A = Simplify[U.S.Lam.Inverse[S].Inverse[U]];  Vi
= Transpose[Simplify[Eigenvectors[A]]];  Vip = Simplify[Vi.Inverse[T]];
Simplify[Solve[Vip[[1,1]]==k1 && Vip[[1,2]]==k2, {w,s}]]
---

∎

###  A.4 9D linear system with non-resonant modes

We present a 9D linear system with linear observables that having spectrum
similar to the right 9 plots of figure 6. Each eigenvalue differs less than
0.01. The eigenvalue density plots in figure 9 and figure 6 are similar as
well. This shows the eigenvalue estimation error is not cause by eigenvalue
resonance with second order observables, but the number and spectrum distance.

Figure 9: Eigenvalue density plot using slightly perturbed non-resonant
eigenvalues. Different rows have different trajectory length, and different
columns have similar eigenvalue distribution as right 9 plots of figure 6. Red
crosses are ground truth eigenvalues.

###  A.5 Other DMD algorithms

We present the results with the same setting as in section 5, using forward-
backward (FB) DMD [49], total-least square (TLS) DMD [50] and optimized DMD
[51]. The difficulties in spectrum recovery, and sensitivity to non-normality
and higher order observables persist. The optimized DMD does not support
multiple trajectories as input data, so it only uses one trajectory with
specified length throughout the plots.

In figure 10-figure 13, we show the contour plots with the same system as the
first and fourth column of figure 2 and figure 4. FB DMD and TLS DMD can
correct the underestimation bias caused by observation noise in exact DMD [7],
but they overestimate the modes under both system and observation noise in
figure 10. We also did not see bias in exact DMD on a dataset with more
initial conditions. In figure 12, FB DMD removes the multi-modal error, and
gives better estimation on the real eigenvalue with $L=10,50$. In figure 11
and figure 13, all the algorithms encountered difficulties to recovery the
full spectrum with second order observables. Optimized DMD with $L=50$ in
figure 11 gives correct frequency estimation for the slowest decaying
oscillatory mode (purple contour), and in figure 13 the purple contour is
surrounding all of the oscillatory modes. This might indicate the slowest
decaying mode from Optimized DMD can always give some reliable information
about the dynamics. Furthermore, Optimized DMD only used one trajectory with
length 50, i.e. half of the data comparing with the other algorithms. A
promising future work is to make Optimized DMD support multiple trajectories.

Figure 10: Eigenvalue density plot with first order observables. System $A$
matrix is constructed by $\theta=0,\phi=0,s=1$. Different rows have different
trajectory length, and different columns use different algorithms. Red crosses
are ground truth eigenvalues. Optimized DMD uses only 1 trajectory. Figure 11:
Eigenvalue density plot with monomial observables up to 2nd order. System $A$
matrix is constructed by $\theta=0,\phi=0,s=1$. Different rows have different
trajectory length, and different columns use different algorithms. Red crosses
are ground truth eigenvalues. Optimized DMD uses only 1 trajectory. Figure 12:
Eigenvalue density plot with first order observables. System $A$ matrix is
constructed by $\theta=1.4,\phi=0,s=0.1$. Different rows have different
trajectory length, and different columns use different algorithms. Red crosses
are ground truth eigenvalues. Optimized DMD uses only 1 trajectory. Figure 13:
Eigenvalue density plot with monomial observables up to 2nd order. System $A$
matrix is constructed by $\theta=1.4,\phi=0,s=0.1$. Different rows have
different trajectory length, and different columns use different algorithms.
Red crosses are ground truth eigenvalues. Optimized DMD uses only 1
trajectory.

###  A.6 5D system

We present a 5D system with spectrum more reminiscent of multi-legged systems,
shown in figure 14.

Figure 14: Eigenvalue density plot with two slow oscillatory modes added (5d
system with linear observation). Different rows have different trajectory
length, and different columns use different algorithms. Red crosses are ground
truth eigenvalues.


# Enhanced Flight Envelope Protection: A Novel Reinforcement Learning Approach

Akin Catak  Ege C. Altunkaya  Mustafa Demir  Emre Koyuncu  Ibrahim Ozkol
Aviation Institute, Aerospace Research Center, Istanbul Technical University,
Istanbul, Türkiye (e-mail<EMAIL_ADDRESS><EMAIL_ADDRESS>Turkish Aerospace, Air Vehicle Technologies Research Center, Istanbul,
Türkiye, (e-mail<EMAIL_ADDRESS> Department of Aeronautics
Engineering, Istanbul Technical University, Istanbul, Türkiye, (e-mail:
<EMAIL_ADDRESS> Department of Aeronautics Engineering, Istanbul
Technical University, Istanbul, Türkiye, (e-mail<EMAIL_ADDRESS>
###### Abstract

This paper introduces a flight envelope protection algorithm on a longitudinal
axis that leverages reinforcement learning (RL). By considering limits on
variables such as angle of attack, load factor, and pitch rate, the algorithm
counteracts excessive pilot or control commands with restoring actions. Unlike
traditional methods requiring manual tuning, RL facilitates the approximation
of complex functions within the trained model, streamlining the design
process. This study demonstrates the promising results of RL in enhancing
flight envelope protection, offering a novel and easy-to-scale method for
safety-ensured flight.

###### keywords:

Flight envelope protection, reinforcement learning, nonlinear flight control.

††thanks: This work has been submitted to IFAC for possible publication.

##  1 Introduction

Ensuring the safety of an aircraft’s operational capabilities has always been
an irreplaceable priority. Regardless of the class of the aircraft, stability
and safety must be proven at every single operation point within the flight
envelope; for this purpose, an intensive effort is made throughout the design.
Nevertheless, the control system design undertakes a significant portion of
this workload, particularly in the design of fighter jets, which have been
designed as aerodynamically unstable for over four decades. Within this
context, the design of flight envelope protection algorithms is one of the
spearhead requirements in terms of ensuring both safety and stability.

The principal purpose of flight envelope protection algorithms is averting the
aircraft from the deviation of the predetermined safe flight envelope, thereby
mitigating the risk of entering an upset condition. Generally, the objective
is achieved by limiting the angle of attack, sideslip angle, load factor,
angular rates, bank angle, and true airspeed. The distinctions between the
algorithms arise at the point of how to limit the aircraft within the safe
boundaries and how to recover the aircraft if a violation occurs.

###  1.1 State-of-Art

Various approaches to flight envelope protection are available, including
adaptive flight envelope protection algorithms Lombaerts et al. (2017), haptic
feedback design Van Baelen et al. (2020), reachability analysis Yu et al.
(2020), command and reference governor techniques Simon et al. (2017); Ye et
al. (2015); Wang et al. (2022), reference command regeneration Yu et al.
(2022), control barrier functions Altunkaya et al. (2024), and traditional
flight envelope protection algorithms Moreira et al. (2022); Nguyen (1980).
Specifically, the study in Yu et al. (2020) introduces an online method for
determining the flight envelope using reachability analysis and safeguarding
it with optimal control inputs under icing conditions. Given that icing
severely degrades the aerodynamic performance of the aircraft, the reliability
of a conventional flight envelope diminishes; thus, it is imperative to assess
external factors to avoid departing from the safe flight envelope. The
proposed method employs a protection strategy that requires meticulous design
and validation to ensure effectiveness under all possible scenarios.
Additionally, a command governor approach is suggested in Simon et al. (2017)
to limit maneuvers based on load factor and angle of attack. This protection
design is tested in a hardware-in-the-loop environment to demonstrate its
efficacy; however, a significant limitation of this method is the need for
scheduling the weights in the constructed objective function. This scheduling
issue is similarly noted in Wang et al. (2022), where there is a necessity for
the scheduling of the command governor matrices.

###  1.2 Problem Statement

Even though there is a diverse variety of flight envelope protection
algorithms in the literature, their industrial design relies on classical
methods. However, the design of classical flight envelope protection
algorithms requires careful, experienced, and wise evaluation to certify
safety within the entire flight envelope of the aircraft. Otherwise, improper
design may result in high pilot workload, unacceptable flying quality, and
even instability. As a result, the conventional design approach can be
regarded as a tedious and serious workload.

In this regard, the study aims to address both safety and pilot workload
issues while proposing a fast and reliable solution approach. Therefore, a
novel flight envelope protection algorithm using reinforcement learning is
presented to address the aforementioned issues effectively, and as the primary
step of the research, a flight envelope protection algorithm just for the
angle of attack and load factor is developed. As a consequence of the
research, a care-free, i.e., less pilot workload, and safety-ensured flight
envelope protection system is expected.

###  1.3 Contributions

The contributions of this study are itemized as follows;

  * •

A novel flight envelope protection algorithm leveraging reinforcement learning
is introduced.

  * •

The issues encountered throughout the design of conventional flight envelope
protection algorithms are mitigated with the proposed method.

  * •

The performance deterioration of traditional flight envelope protection
algorithms during extreme maneuvers, where coupling effects and nonlinearities
arise remarkably, is mitigated with the proposed method.

Consequently, the efficacy of the proposed method is assessed under diverse
circumstances.

##  2 Preliminaries

The fundamental requirements prior to constructing the principal methodology
are briefly presented in this section.

###  2.1 Nonlinear Flight Dynamics

The minimum needs for developing a nonlinear flight dynamics model for F-16 is
discussed under two basic headlines: the equations of motion as well as
aerodynamics and actuators.

####  2.1.1 Equations of Motion

The nonlinear rigid body dynamics are given in a compact form in Eq. (1) and
Eq. (2).

$\bm{\dot{V}}=m^{-1}[\bm{F}-\bm{\omega}\times m\bm{V}]$ (1)

$\bm{\dot{\omega}}=J^{-1}[\bm{M}-\bm{\omega}\times J\bm{\omega}]$ (2)

where $\bm{V}\in\mathbb{R}^{3\times 1}$ is the body velocity vector,
$\bm{\omega}\in\mathbb{R}^{3\times 1}$ is the angular rate vector,
$\bm{F}\in\mathbb{R}^{3\times 1}$ is the total body force vector,
$\bm{M}\in\mathbb{R}^{3\times 1}$ is the total moment vector,
$J\in\mathbb{R}^{3\times 3}$ is the inertia tensor, and $m$ is the mass of the
aircraft. Additionally, the rotational kinematics are given by Eq. (3).

$\bm{\dot{\Omega}}=\begin{bmatrix}1&\sin\phi\tan\theta&\cos\phi\tan\theta\\\
0&\cos\phi&-\sin\phi\\\
0&\sin\phi\sec\theta&\cos\phi\sec\theta\end{bmatrix}\bm{\omega}$ (3)

where $\bm{\Omega}$ denotes the Euler angles.

####  2.1.2 Aerodynamics and Actuators

The aerodynamic modeling consists of the scheduled wind-tunnel data, provided
in Nguyen (1979), and their corresponding formulation to obtain the force
($C_{x},C_{y},C_{z}$) and moment coefficients ($C_{l},C_{m},C_{n}$). After
calculating the coefficients, the necessary dimensionalization is processed
using the current states and geometric properties of the aircraft, i.e.
dynamic pressure, wing area, span, and mean aerodynamic chord.

In addition, each control surface has a unique actuator modeling, i.e. a
first-order model with a time constant $0.0495s$. Furthermore, the rate
saturation values are $60^{\circ}/s$, $80^{\circ}/s$, and $120^{\circ}/s$,
whereas the position saturation values are $\pm 25^{\circ}$, $\pm
21.5^{\circ}$, and $\pm 30^{\circ}$ for the horizontal tail, aileron, and
rudder respectively.

###  2.2 Flight Control Laws

The control augmentation system consists of a single-loop angular rate control
law using the incremental nonlinear dynamic inversion (INDI). Prior to
deriving the control law, the control-affine form of the Euler’s equations of
motion is given by Eq. (4).

$\bm{\dot{\omega}}=-J^{-1}(\bm{\omega}\times
J\bm{\omega})+\underbrace{J^{-1}\bar{q}_{\infty}S\begin{bmatrix}b&&\\\
&\bar{c}&\\\ &&b\\\
\end{bmatrix}\Phi}_{\begin{subarray}{c}\bm{g}(\bm{x})\end{subarray}}\underbrace{\bm{\delta}}_{\begin{subarray}{c}\bm{u}\end{subarray}}$
(4)

where $\bar{q}_{\infty}$, $S$, $b$, and $\bar{c}$ are the dynamic pressure,
wing area, wing span, and mean aerodynamic chord, respectively. Additionally,
$\Phi\in\mathbb{R}^{3\times n}$ denotes the control effectivity matrix, with
$n$ is the number of control surfaces. Finally, the control law for the
angular rates is presented in Eq. (5).

$\bm{u}=\bm{g}(\bm{x}_{0})^{-1}[\bm{\dot{\omega}}_{c}-\bm{\dot{\omega}}_{0}]+\bm{u}_{0}$
(5)

where the subscript $"0"$ denotes the current state and
$\bm{\dot{\omega}}_{c}\in\mathbb{R}^{3}$ is the virtual input to be designed.
The exact form of the control law is given by Eq. (6).

$\bm{\delta}=\Bigg{\\{}J^{-1}\bar{q}_{\infty}S\begin{bmatrix}b&&\\\
&\bar{c}&\\\ &&b\\\
\end{bmatrix}\Phi\Bigg{\\}}^{-1}[\bm{\dot{\omega}}_{c}-\bm{\dot{\omega}}_{0}]+\bm{\delta}_{0}$
(6)

where $\bm{\delta}\in\mathbb{R}^{3}$ signifies the control surface
deflections, i.e. aileron, horizontal tail, and rudder, respectively.
Furthermore, the virtual input $\bm{\dot{\omega}}_{c}$ is provided by Eq. (7).

$\bm{\dot{\omega}}_{c}=\begin{bmatrix}K_{p}&&\\\ &K_{q}&\\\ &&K_{r}\\\
\end{bmatrix}\begin{bmatrix}p_{c}-p\\\ q_{c}-q\\\ r_{c}-r\\\ \end{bmatrix}$
(7)

where $K_{p}$, $K_{q}$, and $K_{r}$ represent the gains for the roll, pitch,
and yaw channels, respectively.

Consequently, the aforementioned relations complete the flight control law for
the angular rates.

###  2.3 Classical Flight Envelope Protection

The classical protection algorithm for the angle of attack and load factor is
developed based on a simple rationale: applying a counter-action as a
restorative measure to prevent FEP violations.

####  2.3.1 Angle of Attack and Load Factor Protection

The excess angle of attack and load factor boundaries can be prevented through
the implementation of a ’restorative’ pitch rate command. Before reaching the
predefined limits of the angle of attack and load factor, the pilot commands
are gradually diminished to avert the aircraft’s entry into an upset
condition. Additionally, both angle of attack and load factor protection must
be viewed comprehensively to produce a singular and cohesive restorative pitch
response; hence, the load factor boundary is also translated into the angle of
attack as represented by Eq. (8).

$\alpha_{max}^{n_{z}}=\dfrac{Wn_{z_{max}}}{\bar{q}_{\infty}SC_{z_{\alpha}}}$
(8)

where $\alpha_{max}^{n_{z}}$ denotes the angle of attack equivalent to the
maximum allowable load factor $n_{z_{max}}$, $W$ is the weight of the
aircraft, $C_{z_{\alpha}}$ is the force coefficient in z-direction derivative
with the angle of attack. The general framework is depicted in Fig. 1.

Figure 1: The classical angle of attack and load factor protection scheme.

Consequently, the classical flight envelope protection for the angle of attack
and load factor is completed.

###  2.4 Reinforcement Learning

Reinforcement Learning (RL) is a subfield of Machine Learning in which the
problem is represented as a Markov decision Process (MDP) (Van Otterlo and
Wiering (2012)), and the actions taken by the agent are evaluated utilizing
Belman optimality conditions. The objective of the RL is to maximize the
cumulative reward by finding the optimal policy of the agent while dynamically
interacting with the environment. An agent can be described as anything that
learns and makes decisions while interacting with the environment. The
environment is what the agent interacts with, basically everything except the
agent. In RL, learning is achieved in an iterative way by evaluating the
sequence of actions and eventually maximizing the cumulative rewards by
updating the policy such that in each state, the policy selects the near-
optimal actions as training progresses.

In MDP every step can be represented with the tuple $<s,a,p,r>$ consisting of
environments state ($s\in\mathcal{S}$), action ($a\in\mathcal{A}$), transition
function $p$ ($0\leq p\leq 1$ ) between state $s$ and state $s^{\prime}$, and
reward $r$ for taking the action $a$ that takes from state $s$ to state
$s^{\prime}$. Another important concept is the policy $\pi(s|a)$, which can be
thought of as a strategy. The policy determines the probability of each action
at each stage. The objective of Reinforcement learning is to find an optimal
policy $\pi^{*}(s|a)$ in which the cumulative reward is maximized. The Value
function is the combination of immediate Reward and future rewards. Thus, in
other words, the objective can be thought of as finding an optimal policy that
maximizes the expected value function.

###  2.5 Actor Critic - Deep Deterministic Policy Gradient

Actor-critic methods are the natural extension of the idea of reinforcement-
comparison methods to temporal difference learning and to the full
reinforcement-learning problem (Konda and Tsitsiklis (1999)). It combines the
benefits of both value-based and policy-based RL methods. Typically, the
critic is a state-value function. After each action selection, the critic
evaluates the new state to determine whether things have gone better or worse
than expected. Deep Deterministic Policy Gradient is one of the actor-critic
algorithms that is used for continuous actions (Lillicrap et al. (2015); Duan
et al. (2016)). It approximates actor and critic using deep neural networks
and employs deterministic policy gradient for learning.

The action determined by the actor with the policy parameters $(\theta^{\pi})$
is given in Eq. (9) with the added stochastic noise $\mathcal{N}_{t}$

$a_{t}=\pi(s_{t}|\theta^{\pi})+\mathcal{N}_{t}$ (9)

The critic network is responsible for evaluating the Q-value for the state-
action pairs. Critic calculates the temporal difference(TD) error,
$\delta_{t}$, using the immediate reward ($r_{t}$), current Q-Value, and
target models Q-value $Q^{\prime}$ as given in Eq. (10).

$\delta_{t}=r_{t}+\gamma
Q^{\prime}(s_{t+1},\pi^{\prime}(s_{t+1}|\theta^{\pi^{\prime}})|\theta^{Q^{\prime}})-Q(s_{t},a_{t}|\theta^{Q})$
(10)

where the $Q^{\prime}$ and $\pi^{\prime}$ are target models with their
parameters $\theta^{Q^{\prime}}$ and $\theta^{\pi^{\prime}}$, and $\gamma$
being discount factor, basically determining the importance of the future
states.

Critic parameters $\theta^{\pi}$ are updated by minimizing the loss across all
samples as in Eq. (11).

$L=\frac{1}{2M}\sum_{i=1}^{M}\left(\delta_{t}\right)^{2}$ (11)

where $M$ is the mini-batch size.

The actor parameters $\theta^{\pi}$ are updated using the gradient given in
Eq. (12). To maximize expected discounted reward

$\nabla_{\theta^{\pi}}J\approx\frac{1}{2M}\sum_{i=1}^{M}\left(\nabla_{a}Q(s,a|\theta^{Q})|_{a=\pi(s|\theta^{\pi})}\nabla_{\theta^{\pi}}\pi(s|\theta^{\pi})\right)$
(12)

where $\nabla_{a}Q(s,a|\theta^{Q})$ is the gradient of the critic output with
respect to action ${a=\pi(s|\theta^{\pi})}$ generated by actor network and
$\nabla_{\theta^{\pi}}\pi(s|\theta^{\pi})$ is the gradient of the actor
network’s output with respect to the actor parameters $\theta^{\pi}$.

The target networks are updated with smoothing factor $\tau$ as given in Eq.
(13).

$\displaystyle\theta^{Q^{\prime}}\leftarrow\tau\theta^{Q}+(1-\tau)\theta^{Q^{\prime}}$
(13)
$\displaystyle\theta^{\pi^{\prime}}\leftarrow\tau\theta^{\pi}+(1-\tau)\theta^{\pi^{\prime}}$
(14)

where small $\tau=0.001$ enables smooth updates for target networks.

##  3 Methodology and Reward Structure

The proposed method focuses on developing a flight envelope protection scheme
given in Fig. 2 using RL. In this method, the environment consists of the
nonlinear F-16 model and the flight control law. This environment enables
capturing the complex flight dynamics and control challenges. Since we have
continuous actions, the Deep Deterministic Policy Gradient (DDPG) algorithm,
an off-policy actor-critic method that is also sample efficient, is employed.
The action of the proposed envelope protection scheme in the longitudinal axis
is the restorative action $q_{rest}$ within the range $-20^{\circ}/s\leq
q_{rest}\leq 30^{\circ}/s$.

Figure 2: Actor-critic scheme.

For the proposed method, observation space should contain the relevant states
of the environment, in this case, the nonlinear F-16 model.

$S=[n_{z},\alpha,\phi,p,q,e_{q},\bar{q}_{\infty}]$ (15)

where states $n_{z},\alpha$ and $q$ that need to be protected are included as
well as $\phi$ and $p$ due to coupling affects. The tracking error,
$e_{q}=q_{cmd}-q$, is included to give information about tracking error and
dynamic pressure $\bar{q}_{\infty}$ is also included since its substantial
effects over the dynamics. In training, the aircraft starts from level flight
trim condition with 0.6 Mach speed and 500-meter altitude. Each episode is 10
seconds long and has constant pilot input on the pitch channel determined
randomly at the beginning of each episode. The training is stopped when the
average reward of 150 episodes reaches the desired stopping value.

The reward used in the proposed method is in the form of feedback taken from
the environment, which evaluates how good or bad the current policy is given
the desired objectives. For our problem, the Reward function is in the form of
Eq. (16).

$R=r_{s}+W_{t}r_{t}+W_{a}r_{a}+W_{n}r_{n}+W_{q}r_{q}+r_{p}$ (16)

where $r$s represent individual rewards and $w$s are corresponding weights.
$r_{s}$ is the positive constant reward obtained at each step supporting the
agent to survive. $r_{t}$ is the tracking cost added for increasing the track
performance. $r_{a}$ is the cost of the angle of attack excess, which will be
defined later. Similar to $r_{t}$, $r_{n}$ and $r_{q}$ are the costs for
exceeding load factor and pitch rate limits respectively. The last component
of $R$ is the $r_{p}$, which is the so-called isdone reward penalizing the
worst policies.

Reward related to survival is $r_{s}=0.1$. Since it is a constant it does not
have a weight associated with it. The cost related to tracking performance is
given in Eq. (17).

$r_{t}=\left(\frac{\left(|q|-\left|q_{cmd}\right|\right)}{\left|q_{cmd}\right|+\epsilon}\right)^{2}$
(17)

where $q$ is the pitch rate, $q_{cmd}$ is commanded pitch rate and $\epsilon$
is a small number for avoiding zero division.

The cost related to exceeding the angle of attack limit is given in Eq. (18).

$r_{a}=\left\\{\begin{array}[]{cl}-\left(\frac{\left(|\alpha|-0.9\left|\alpha_{\max}\right|\right)}{0.9\left|\alpha_{\max}\right|}\right)^{2}&\text{if
}\left|\alpha\right|\geq 0.9\left|\alpha_{\max}\right|\\\ 0&\text{ otherwise
}\end{array}\right.$ (18)

where $\alpha$ is the angle of attack and its limit is $\alpha_{max}$.
Multiplication with $0.9$ is added to the reward to prevent the agent from
getting close to the limit. The cost related to exceeding the load factor
limit is given in Eq. (19).

$r_{n}=\left\\{\begin{array}[]{cl}-\left(\frac{\left(|n_{z}|-\left|n_{zmax}\right|\right)}{\left|n_{zmax}\right|}\right)^{2}&\text{if
}\left|n_{z}\right|\geq\left|n_{zmax}\right|\\\ 0&\text{ otherwise
}\end{array}\right.$ (19)

where $n_{z}$ is the load factor component and $n_{zmax}$ is the maximum
allowable load factor. Similarly, the pitch rate cost can be defined as well
in Eq. (20).

$r_{q}=\left\\{\begin{array}[]{cl}-\left(\frac{\left(|q|-\left|q_{max}\right|\right)}{\left|q_{max}\right|}\right)^{2}&\text{if
}\left|q\right|\geq\left|q_{max}\right|\\\ 0&\text{ otherwise
}\end{array}\right.$ (20)

The last component of the $R$ is the penalty cost, which is designed to
penalize the excess of the envelope for a certain time as given in Eq. (21).

$r_{p}=\left\\{\begin{array}[]{ll}-400&\text{if
}cond_{1}\operatorname{holds}\\\ -600&\text{if
}cond_{2}\operatorname{holds}\end{array}\right\\}$ (21)

where the ${cond}_{1}$ is the first condition indicating at least 2 of the
envelope parameter($\alpha,n_{z},q$) exceed the limit for 2 seconds and
${cond}_{2}$ is the second condition indicating at least one of the envelope
parameter($\alpha,n_{z},q$) exceed their corresponding limits by $50\%$. If
this cost takes place then the episode is terminated.

##  4 Results

In this work, training was conducted using MATLAB®’s Reinforcement Learning
Toolbox. The computer on which the training was conducted has an Intel i7
$\nth{5}$ series processor with a clock speed of 2.8 GHz and 16 GB of memory.
The noise variance is taken as 0.001, the variance decay rate is 1e-9, the
actor learn rate is 1e-5, the critic learn rate is 1e-4, and the discount
factor is trained as 0.99. The critic network has observation and action paths
with having 2 and 1 hidden layer, respectively. The hidden layer on the action
path has 40 neurons, while hidden layers on the observation path have 80 and
40 neurons, respectively. Those layers are added together to obtain the Q
value. The actor network has one hidden layer with 40 neurons.

The performance of the proposed envelope protection method is examined under
three different pilot commands as tabulated in Table 1.

Table 1: Simulation scenarios | Command
---|---
| $p,^{\circ}/s$ | $q,^{\circ}/s$ | $r,^{\circ}/s$
Scenario-$\\#1$ | 0 | 25 | 0
Scenario-$\\#2$ | 0 | -10 | 0
Scenario-$\\#3$ | 60 | $q(t)$ | 0

All scenarios include the aforementioned classical envelope protection
algorithm as a benchmark to compare. The first scenario includes
$q_{cmd}=25^{\circ}/s$, which is the maximum pilot command used in the
training, and it is evaluated to show the performance of the proposed method
in the maximum pilot command, as depicted in Fig. 3.

Figure 3: The comparison of envelope protection algorithms: scenario-$\\#1$.

The red lines indicate the limits for each parameter to be protected. The
abbreviation RL in the legend denotes the proposed method, whereas CM stands
for the classical approach for the sake of comparison in terms of flight
envelope protection performances. This scenario is the most challenging
because the angle of attack protection and command tracking performance
metrics are conflicting. While the tracking performance is better in the
initial phase of the command for the CM, the tracking performance is better
for the rest of the simulation, and the angle of attack limit is used more
aggressively in the RL method. Both methods can protect the flight envelope in
this scenario, but when performance requirements arise, the results are
debatable. The results of the second scenario are illustrated in Fig. 4.

Figure 4: The comparison of envelope protection algorithms: scenario-$\\#2$.

The second scenario is conducted to evaluate the maximum negative pilot
command used in the training of the RL, which is actually the negative pitch
rate limit of the aircraft. Both methods can protect the envelope but the
tracking performance of the RL is slightly better in addition to the better
usage of the limits. The third and last scenario is presented in Fig. 5.

Figure 5: The comparison of envelope protection algorithms: scenario-$\\#3$.

The third scenario serves as a litmus test, where the proposed RL-based
approach is noticeably better in terms of performance and envelope protection.
In the initial part of the simulation, involving a positive pitch rate coupled
with a roll rate, both methods perform satisfactorily, though CM is slightly
better in tracking performance. The crucial part of the simulation is the
transition from a positive pitch rate to a negative pitch rate command. This
transition is successfully handled by the proposed RL method, while the
classical approach fails to protect the envelope. Besides the poor performance
in envelope protection, the tracking performance of the classical method is
noticeably worse than that of the proposed approach.

Furthermore, Fig. 6 shows the results for different constant pitch rate
commands, demonstrating the proposed method’s efficacy in envelope protection
under varying circumstances. Each simulation is conducted for different pitch
rate commands over an interval of $q_{cmd}=[-10^{\circ}/s\hskip
4.26773pt25^{\circ}/s]$ with an increment of $\Delta q_{cmd}=0.5^{\circ}/s$.

Figure 6: Monte Carlo simulation results: pitch-only assessments.

The lowest pitch rate command exceeds the negative limit for a short time;
however, since the metric for failure is defined as exceeding the limit for at
least 2 seconds, this excess is acceptable.

Finally, the same analysis is conducted to see the protection performance
under pitch-roll coupled maneuver circumstances. In this regard, Fig. 7 is
depicted to show the efficacy of the proposed method.

Figure 7: Monte Carlo simulation results: pitch-roll coupled assessments.

The addition of roll rate command induces oscillation in $\alpha$, $n_{z}$,
and $q$ values. Although the $p$ and $\phi$ are inside the observation space,
the weights of the neurons corresponding to these states are relatively small
due to zero $p_{cmd}$ throughout the simulated maneuver. It is safe to
conclude that the performance of the proposed method is satisfactory in the
presence of coupled commands since all the commands are handled except the two
red lines.

##  5 Conclusion

In this work, a novel flight envelope protection algorithm for the
longitudinal motion of the aircraft is introduced based on reinforcement
learning. The algorithm allows the aircraft to conduct considerably aggressive
maneuvers without violating the predetermined safe flight envelope.
Furthermore, the proposed method surpasses the designed classical flight
envelope protection algorithm, especially in pitch-roll coupled maneuvers. As
a final assessment, two distinct Monte Carlo simulations are carried out in
order to prove the efficacy of the proposed protection algorithm. Impressions
from the simulations are not only satisfactory but also show a significant
potential for more advanced envelope protection algorithms building on the
proposed approach, contributing remarkably to carefree flight. Future plans
encompass the studies based on a complete and coherent envelope protection
algorithms for both longitudinal and lateral motion to ensure safety, even in
the most complex, challenging, and aggressive maneuvers.

## References

  * Altunkaya et al. (2024) Altunkaya, E.C., Çatak, A., Demir, M., Koyuncu, E., and Özkol, İ. (2024).  Stability and safety assurance of an aircraft: A practical application of control lyapunov and barrier functions.  _Available at SSRN 4823223_.
  * Duan et al. (2016) Duan, Y., Chen, X., Houthooft, R., Schulman, J., and Abbeel, P. (2016).  Benchmarking deep reinforcement learning for continuous control.  In _International conference on machine learning_ , 1329–1338. PMLR.
  * Konda and Tsitsiklis (1999) Konda, V. and Tsitsiklis, J. (1999).  Actor-critic algorithms.  _Advances in neural information processing systems_ , 12.
  * Lillicrap et al. (2015) Lillicrap, T.P., Hunt, J.J., Pritzel, A., Heess, N., Erez, T., Tassa, Y., Silver, D., and Wierstra, D. (2015).  Continuous control with deep reinforcement learning.  _arXiv preprint arXiv:1509.02971_.
  * Lombaerts et al. (2017) Lombaerts, T., Looye, G., Ellerbroek, J., and Martin, M.R.y. (2017).  Design and piloted simulator evaluation of adaptive safe flight envelope protection algorithm.  _Journal of Guidance, Control, and Dynamics_ , 40(8), 1902–1924.
  * Moreira et al. (2022) Moreira, M.A., Gripp, J.A., Yoneyama, T., and Marinho, C.M. (2022).  Longitudinal flight control law design with integrated envelope protection.  _Journal of Guidance, Control, and Dynamics_ , 45(9), 1739–1749.
  * Nguyen (1979) Nguyen, L.T. (1979).  _Simulator study of stall/post-stall characteristics of a fighter airplane with relaxed longitudinal static stability_ , volume 12854.  National Aeronautics and Space Administration.
  * Nguyen (1980) Nguyen, L.T. (1980).  _Control-system techniques for improved departure/spin resistance for fighter aircraft_ , volume 791083.  National Aeronautics and Space Administration, Scientific and Technical Information Branch.
  * Simon et al. (2017) Simon, D., Härkegård, O., and Löfberg, J. (2017).  Command governor approach to maneuver limiting in fighter aircraft.  _Journal of Guidance, Control, and Dynamics_ , 40(6), 1514–1527.
  * Van Baelen et al. (2020) Van Baelen, D., Ellerbroek, J., Van Paassen, M., and Mulder, M. (2020).  Design of a haptic feedback system for flight envelope protection.  _Journal of Guidance, Control, and Dynamics_ , 43(4), 700–714.
  * Van Otterlo and Wiering (2012) Van Otterlo, M. and Wiering, M. (2012).  Reinforcement learning and markov decision processes.  In _Reinforcement learning: State-of-the-art_ , 3–42. Springer.
  * Wang et al. (2022) Wang, L., Yang, K., Liu, H., and Yue, T. (2022).  Command governor for constrained attitude angle protection of the wing-in-ground effect craft near sea surface.  _Journal of Aerospace Engineering_ , 35(6), 04022080.
  * Ye et al. (2015) Ye, H., Chen, M., Wu, Q., et al. (2015).  Flight envelope protection control based on reference governor method in high angle of attack maneuver.  _Mathematical Problems in Engineering_ , 2015.
  * Yu et al. (2022) Yu, X., Zhou, X., Zhang, Y., Guo, L., Ye, S., and Peng, X. (2022).  Safety control design with flight envelope protection and reference command generation.  _IEEE Transactions on Aerospace and Electronic Systems_ , 58(6), 5835–5848.
  * Yu et al. (2020) Yu, Z., Li, Y., Zhang, Z., Xu, W., and Dong, Z. (2020).  Online safe flight envelope protection for icing aircraft based on reachability analysis.  _International Journal of Aeronautical and Space Sciences_ , 1–11.