# Bounded nonlinear forecasts of partially observed geophysical systems with
physics-constrained deep learning

Said Ouala1, Steven L. Brunton2, Bertrand Chapron3
Ananda Pascual4, Fabrice Collard5, Lucile Gaultier5, Ronan Fablet1
$(1)$ IMT Atlantique; Lab-STICC, 29200 Brest, France
$(2)$ University of Washington, USA
$(3)$ Ifremer, LOPS, 29200 Brest, France
$(4)$ IMEDEA, UIB-CSIC, 07190 Esporles, Spain
$(5)$ ODL, 29200 Brest, France
{said.ouala<EMAIL_ADDRESS>
<EMAIL_ADDRESS><EMAIL_ADDRESS>
<EMAIL_ADDRESS>
{dr.fab<EMAIL_ADDRESS>
###### Abstract

The complexity of real-world geophysical systems is often compounded by the
fact that the observed measurements depend on hidden variables. These latent
variables include unresolved small scales and/or rapidly evolving processes,
partially observed couplings, or forcings in coupled systems. This is the case
in ocean-atmosphere dynamics, for which unknown interior dynamics can affect
surface observations. The identification of computationally-relevant
representations of such partially-observed and highly nonlinear systems is
thus challenging and often limited to short-term forecast applications. Here,
we investigate the physics-constrained learning of implicit dynamical
embeddings, leveraging neural ordinary differential equation (NODE)
representations. A key objective is to constrain their boundedness, which
promotes the generalization of the learned dynamics to arbitrary initial
condition. The proposed architecture is implemented within a deep learning
framework, and its relevance is demonstrated with respect to state-of-the-art
schemes for different case-studies representative of geophysical dynamics.

_Keywords_ Partially-observed systems $\cdot$ Embedding $\cdot$ Boundedness
$\cdot$ Deep learning $\cdot$ Neural ODE $\cdot$ Forecasting

##  1 Introduction

Modeling the dynamics of complex systems is one of the key drivers of
technological development. Dynamical systems provide interpretable models of
the evolution of real-world phenomena, enabling understanding and prediction.
There are a number of practical and theoretical limitations that arise when
characterizing and modeling physical phenomenon, such as the lack of known
governing equations [1], the presence of hidden variables that are not
measured directly [2], and difficulties in assessing and enforcing known
symmetries and stability [3, 4, 5]. These limitations motivate the development
of new techniques to advance the current state-of-the-art in model discovery.
The current literature is largely based on optimization techniques [6] and
artificial intelligence (AI) [7] algorithms with the aim of discovering new
models, governing equations, conservation laws and symmetries for complex
systems [8, 9, 10, 11, 1, 12, 13] for a variety of applications, ranging from
system identification [14, 1, 15, 16, 17, 18], forecasting [19, 20] and
reconstruction [21, 22, 23] to control [24, 25].

One of the new challenges of data-driven methods and learning based approaches
concerns the exploitation and development of state-of-the-art techniques to
tackle the modeling of complex systems such as the atmosphere, the ocean, and
the climate. These fields encountered a sensing revolution in the last years,
which sets up the community with databases that enable the use of AI. In this
context, significant advances in model reduction techniques [26, 27],
interpolation and reconstruction of gap-free observations [21, 28, 29] and
short term forecast [30] have been achieved by AI models. These models, when
compared to classical physical models, issued for instance from the Navier-
Stokes equations, are easier to handle from a computational point of view,
they also often offer a more flexible application range with the drawback of
been bad learners. AI models often fail generalizing to situations (or initial
conditions) that were never seen in the training data. These model also often
fail at matching physical model in their ability to simulate the dynamics of
the observations. From a pure learning perspective, generalization, and long
term simulation performance of data-driven schemes are linked to physical
properties of the underlying dynamics (such as conserved quantities) and most
of the time, these properties are not revealed by the complexity of the
training data.

In this context, a central question concerns the ability of data-driven models
to capture the behavior of a system from partial and imperfect observations.
Takens theorem [31] states the conditions under which a delay embedding
representation of the observations is diffeomorphic to the underlying hidden
states. Such a representation, when determined, captures essential information
about the hidden states, which makes it possible to learn a model [2].
Unfortunately, the derivation of a dynamical system from such a representation
remains a challenge, since no explicit relationships between the reconstructed
phase space variables and the original hidden variables have been clearly
identified. To address this limitation, the Neural embedding of Dynamical
Systems (NbedDyn), proposed in [32], allows for both reconstructing and
forecasting the phase space of the unseen dynamical system. However, this
model suffers from generalization and stability issues, especially for initial
conditions that are not near the data used to train the model.

It has been shown that embedding prior physical knowledge, such as
invariances, symmetries, and conservation laws, into the learning process of
dynamical systems may improve the generalizability and stability of the
resulting models [3, 33, 34, 35, 36, 4]. Within this new field of research,
there is an increasing focus on incorporating stability and boundendess
constraints into data-driven models. This interest is driven by a combination
of works, on different classes of models ranging from sparse regression
techniques [4] to deep learning architectures [37, 38], showing that such
constraints can be necessary in order to promote stability and generalization
of the models.

Following these developments, we propose a robust framework for the modeling
of geophysical dynamical systems from data. This framework relies on the
learning of an embedding of the observations governed by a neural ODE (NODE)
model [13]. This model is constrained to be bounded by restricting the form of
the NODE to be linear-quadratic which enables the use of a generalization of
the direct Lyapunov method given by the Schlegel boundedness theorem [39],
which has recently been used to assess the stability of data driven models [5]
and has been incorporated into sparse nonlinear models [4]. On different case-
studies, we show that once these constraints are satisfied, the trained
representation can reproduce realistic trajectories with respect to the
training data using a closed-loop prediction setting. The boundedness
constraint also improves the generalization capabilities of the trained model
for arbitrary initial conditions, even those lying outside the attractor
spanned by the training data. This latter property is extremely important,
since the generalization performance is a critical feature for data-driven
schemes. Overall, key aspects of the proposed framework are three-fold:

  * •

we define a physics-constrained learning-based workflow for the modeling and
forecasting of partially observed dynamical systems, with an emphasis on
geophysical dynamics;

  * •

the proposed framework relies on learning an implicit higher-dimensional
embedding. The dynamics of this embedding are constrained to be bounded, which
promotes the generalization of models for arbitrary initial conditions, even
those far from the training data;

  * •

the relevance of the proposed scheme is tested and demonstrated with respect
to state-of-the-art methods for short-term forecasting and long-term
simulation on several case studies, namely the Lorenz 63 and Lorenz 96 systems
and the shallow water equation dynamics. Specifically, our model is compared
to the sparse regression (SR) technique [1], to a recurrent neural network
(RNN), to the Latent-ODE model as presented in [13] and to an extended dynamic
mode decomposition (EDMD) technique based on delay embedding observables [40].
Through all experiments, the proposed framework is the only model able to
achieve simultaneously both an accurate short-term forecast and long-term
simulation of the dynamics.

Regarding the data-driven identification of real-world geophysical dynamics,
we believe that this proposed framework provides a promising approach for
learning consistent models in terms of short-term and long-term forecasts
through the implementation of physical constraints from prior knowledge of the
conservation laws governing the dynamics.

The paper is organized as follows. We start in section 2 by introducing state-
of-the-art literature on data-driven dynamical systems and physics based
learning. We then briefly discuss, in section 3, stability criteria, both from
a dynamical systems theory and an identification perspective. We also
introduce a generalization of the direct Lyapunov theorem as proposed in [39].
Section 4 presents the proposed framework, followed by the experiments and
results in Section 5. We close this paper with conclusions and perspectives
for future works in Section 6.

##  2 Related work

Data-driven identification of dynamical representations. Stated above, recent
advances in learning-based approaches motivate numerous investigations to
apply these tools for the identification of dynamical representations from
data. A full review of state-of-the-art techniques is out of the scope of this
paper, and we concentrate on providing a high-level overview, organized by the
nature of the training data. Specifically, when provided with observations
that form an embedding of a deterministic time-evolving system, several works
including symbolic or dictionary-based approaches [41, 42, 43, 1, 44], neural-
network and deep learning approaches [45, 46, 13, 47] and model-free
representations [48, 12] demonstrated the ability of data-driven techniques to
reverse engineer the governing equations from a sequence of observations.

When provided with noisy and sparse observations, studies have investigated
either deep learning generative architectures, such as variational
autoencoders [49, 50, 51], or classical Bayesian filtering schemes such as the
ensemble Kalman filter [52], in deriving reconstructions of the irregular
observations that can be useful from an identification perspective [53, 54,
17].

When provided with partial observations, i.e. with no one-to-one mapping
between the observation space and the underlying dynamical system, most state-
of-the-art techniques do not apply since the variability of the observations
is influenced by some missing components, making the temporal evolution of the
observations stochastic in the observation space. From this point of view, the
deep learning literature in time series forecasting mainly relies on some sort
of recurrent neural network (RNN) architectures [55, 56, 57] which can be
motivated from a dynamical system perspective as a temporal propagator of a
delay embedding of the observations. Recent work has connected SINDy with
time-delay autoencoders to reconstruct hidden variables, with promising
results [2].

Methods to learn linear propagators of non-linear, and potentially high-
dimensional dynamical systems have also emerged in the context of the modern
Koopman theory [58]. Finding a suitably chosen transformation/embedding of the
observations that can be approximately propagated linearly in time has
motivated several works [59, 60, 61, 62, 63, 64, 40, 65, 66, 67] on a variety
of applications, including control [68] and prediction [65].

Physics-informed data-driven dynamical representations. In the context of
dynamical model identification, several works have shown that including
physical knowledge has a positive impact on the data-driven models [69, 3, 46,
70, 35, 71, 34, 36, 72, 73, 74, 75, 76, 77]. For instance, [1] proposed a
sparse regression framework for the derivation of interpretable dynamical
representations from data and [3] extended this framework to satisfy physical
constraints, such as known energy preserving symmetries in the dynamics. In
neural-networks and deep learning, significant effort has been focused on
understanding residual networks [78, 79] as numerical integration schemes of
differential equations [80, 15, 47, 81, 82]. In this line of study, neural
ordinary differential equations [13, 83, 84, 85, 86, 87] have shown great
success in merging neural networks and ODE identification techniques.

There is also a growing interest in developing hybrid representations that
merge both a physical model and a machine learning component. Typical examples
can be found in the context of closure modeling [88, 89, 90, 91, 92, 27] in
which, broadly speaking, a machine learning model is optimized on data to
correct a physical model.

Stability and boundedness of data-driven dynamical systems. Exploiting
stability constraints when learning dynamical systems from data usually
improves the long-term stability and generalization performance of the trained
models [93, 94, 95, 96, 4]. Recent works further investigated Lyapunov
functions of dynamical systems for control applications [97] or to force the
stability around a global equilibrium point [37]. In this context, [38]
proposed a new framework that guarantees exponential stability of an NODE
model. Promoting Lyapunov stability in autoencoders was also investigated in
[98] with applications to reduced-order models of two-dimensional turbulent
fluid flows. Beyond the direct Lyapunov method, the Schlegel boundedness
theorem [39] was also considered in [4] to promote the global stability of
linear quadratic data-driven dynamical systems.

##  3 From Stability to Boundedness, in the context of identification

This section introduces some common stability criteria for the analysis of
non-linear models. We briefly discuss the implementation of these criteria and
focus on boundedness constraints as we further explore their relevance to
regularize the learning of chaotic dynamics.

###  3.1 Stability of limit-sets

Let us assume a continuous $s$-dimensional dynamical system
$\boldsymbol{\mathrm{z}}_{t}$ governed by an autonomous ODE
$\dot{\boldsymbol{\mathrm{z}}}_{t}=f(\boldsymbol{\mathrm{z}}_{t})$ with
$\Phi_{t}$ the corresponding flow
$\Phi_{t}(\boldsymbol{\mathrm{z}}_{t_{0}})=\boldsymbol{\mathrm{z}}_{t_{0}}+\int_{t_{0}}^{t}f(\boldsymbol{\mathrm{z}}_{w})dw$.
The trajectories of this dynamical system are assumed to be asymptotic to a
limit-set $\mathcal{S}$ of dimension ${p}$ contained in $\mathbb{R}^{s}$.

Stability theory addresses the characterization of the asymptotic behavior of
a set of solutions of a differential equation with respect to a given limit-
set $\mathcal{S}$. Formally, and assuming for the sake of simplicity that
$\mathcal{S}$ is an equilibrium point $\boldsymbol{\mathrm{z}}_{eq}$, we may
distinguish the following stability definitions:

  * •

The equilibrium point $\boldsymbol{\mathrm{z}}_{eq}$ is globally (respectively
locally) stable if $\forall\boldsymbol{\mathrm{z}}_{0}\in\mathbb{R}^{s}$
(respectively
$\boldsymbol{\mathrm{z}}_{0}\in\mathcal{U}\subset\mathbb{R}^{s}$),
$\Phi_{t}(\boldsymbol{\mathrm{z}}_{0})\longrightarrow\boldsymbol{\mathrm{z}}_{eq}+\epsilon$
with $|\epsilon|>0$ and finite;

  * •

The equilibrium point $\boldsymbol{\mathrm{z}}_{eq}$ is globally (respectively
locally) asymptotically stable if
$\forall\boldsymbol{\mathrm{z}}_{0}\in\mathbb{R}^{s}$ (respectively
$\boldsymbol{\mathrm{z}}_{0}\in\mathcal{U}\subset\mathbb{R}^{s}$),
$\Phi_{t}(\boldsymbol{\mathrm{z}}_{0})\longrightarrow\boldsymbol{\mathrm{z}}_{eq}$
as $t\longrightarrow\infty$;

  * •

The equilibrium point $\boldsymbol{\mathrm{z}}_{eq}$ is globally (respectively
locally) exponentially stable if
$\forall\boldsymbol{\mathrm{z}}_{0}\in\mathbb{R}^{s}$ (respectively
$\boldsymbol{\mathrm{z}}_{0}\in\mathcal{U}\subset\mathbb{R}^{s}$),
$|\Phi_{t}(\boldsymbol{\mathrm{z}}_{0})-\boldsymbol{\mathrm{z}}_{eq}|\leq
C|\boldsymbol{\mathrm{z}}_{0}-\boldsymbol{\mathrm{z}}_{eq}|e^{-\alpha t}$ with
$\alpha$ a convergence rate, $C$ a positive constant and $|\cdot|$ a given
norm in $\mathbb{R}^{s}$.

Broadly speaking, we may state that a limit-set is i) stable if two nearby
trajectories stay close by the action of the vector field, ii) asymptotically
stable if a sufficiently close trajectory is attracted to the limit-set, and
finally iii) exponentially stable if a trajectory converges to the limit-set
with an exponential decay rate.

Depending on the limit-set of a given ODE, there are several approaches to
analyze its stability. A great introduction to classical stability criteria is
given in [99], starting from the classical eigenvalues of a linear (or
linearized) system around an equilibrium point and finishing with the Lyapunov
exponents of chaotic attractors. From an identification perspective, and given
some observation data, stability criteria can be used as constraints, either
to reproduce an observed asymptotic behavior or to avoid blowups, as long as
the attractor revealed by the observations is not strange i.e. chaotic.
Chaotic solutions of differential equation are only revealed through criteria
that exploit long-term simulations of the system such as Lyapunov exponents.
Therefore, they cannot be characterized based on the dynamical equation (such
as the eigenvalues of an equilibrium point).

Overall, one cannot explicitly constrain a system to be chaotic, simply by
looking or manipulating its differential equation. This point appears
critical, since numerous geophysical systems can be chaotic. Alternatively, we
propose to relax the problem by moving from constraining chaos to constraining
the boundedness of a system. Such a boundedness constraint applies to every
single limit-set mentioned above and would avoid blowups of data-driven
models. As highlighted by [99], this constraint is also a natural feature of
real-world geophysical systems as, from an experimental perspective, blowups
cannot be observed in nature. For this purpose, we start by introducing the
direct Lyapunov stability criterion, since it provides global stability
properties on non-linear dynamics. A generalization of this approach to the
boundedness of linear quadratic ODEs as proposed in [39] is then presented.

###  3.2 Lyapunov stability of dynamical systems

The direct Lyapunov stability method [100] was introduced to study the
stability of any dynamical system that admits an equilibrium point at the
origin. It uses a scalar function of the state
$V:\mathbb{R}^{s}\longrightarrow\mathbb{R}$ as follows:

$\displaystyle V(\boldsymbol{\mathrm{z}})=0\mbox{ if and only if
}\boldsymbol{\mathrm{z}}=0$ (1) $\displaystyle
V(\boldsymbol{\mathrm{z}})>0\mbox{ if and only if }\boldsymbol{\mathrm{z}}\neq
0$ $\displaystyle\dot{V}(\boldsymbol{\mathrm{z}})\mbox{ }\leq 0\mbox{
}\forall\boldsymbol{\mathrm{z}}\neq 0$

If $V(\boldsymbol{\mathrm{z}})$ satisfies the above conditions
$\forall\boldsymbol{\mathrm{z}}\in\mathbb{R}^{s}$ (respectively
$\forall\boldsymbol{\mathrm{z}}\in\mathcal{U}\subset\mathbb{R}^{s}$) the
system is globally (respectively locally) stable. Furthermore, if
$\dot{V}(\boldsymbol{\mathrm{z}})<0\forall\boldsymbol{\mathrm{z}}\neq 0$ the
asymptotic stability of the system is also guaranteed.

We can apply this criterion to any non-linear model without resorting to any
linearization, which makes it particularly appealing. Yet, finding an
appropriate function is far from being straightforward, and several works
proposed candidate Lyapunov functions for various types of problems [101].

###  3.3 Generalization to boundedness of LQM

The direct Lyapunov stability method is restricted to dynamical systems with
an equilibrium point at the origin. This property is restrictive since it does
not apply to other dynamical regimes such as periodic and chaotic solutions.
Furthermore, the choice of the Lyapunov function being non-systematic, the
exploitation of a direct Lyapunov constraint for the data-driven
identification of dynamical systems is restricted to a small class of
parametric models. Fortunately, Schlegel and Noack [39] proposed a
generalization of the direct Lyapunov method on a class of parametric
differential equations. This work also introduces the associated Lyapunov
function, along with a condition for the existence of a monotonically
attracting trapping region in the phase space. These trapping regions were
later incorporated into nonlinear model identification to improve stability
[4].

A trapping region is a domain in the phase space that locks trajectories of a
dynamical system, i.e., once a trajectory enters a trapping region, it will
stay in this domain as the system evolves [102]. When this region is globally
monotonically attractive, all trajectories in the phase space will converge to
the trapping region. We may point out that a trapping region can contain a
single or multiple limit-sets. The class of models for which the proposed
criterion in [39] is valid are linear quadratic models (LQMs). They can be
encountered, for instance, by a spectral discretization of the Navier-Stokes
equation, which make them appealing in the context of reduced order models
(ROMs). Formally, as proposed in [39], let us rewrite the dynamical system
governing $\boldsymbol{\mathrm{z}}$ as a linear quadratic model

$\dot{\boldsymbol{\mathrm{z}}}_{t}=\boldsymbol{\mathrm{c}}+\boldsymbol{\mathrm{L}}\boldsymbol{\mathrm{z}}_{t}+[\boldsymbol{\mathrm{z}}^{T}_{t}\boldsymbol{\mathrm{Q}}^{(1)}\boldsymbol{\mathrm{z}}_{t},...,\boldsymbol{\mathrm{z}}^{T}_{t}\boldsymbol{\mathrm{Q}}^{(s)}\boldsymbol{\mathrm{z}}_{t}]^{T}$
(2)

where $\boldsymbol{\mathrm{c}}\in\mathbb{R}^{s}$,
$\boldsymbol{\mathrm{L}}\in\mathbb{R}^{s\times s}$ and
$\boldsymbol{\mathrm{Q}}^{(i)}=[q_{i,j,k}]^{s}_{j,k=1},i=1,...,s$. The $s$
symmetric quadratic matrices $\boldsymbol{\mathrm{Q}}^{(i)}$, $i=1,...,s$ are
supposed to be energy preserving, i.e.

$q^{i}_{j,k}+q^{j}_{i,k}+q^{k}_{i,j}=0,i,j,k=1,...,s$ (3)

Let us also consider a shifted variable
$\bar{\boldsymbol{\mathrm{z}}}=\boldsymbol{\mathrm{z}}-\boldsymbol{\mathrm{m}}$
with $\boldsymbol{\mathrm{m}}\in\mathbb{R}^{s}$ an arbitrary finite state. The
dynamical equation of the shifted state can be written as

$\dot{\bar{\boldsymbol{\mathrm{z}}}}_{t}=\boldsymbol{\mathrm{d}}+\boldsymbol{\mathrm{A}}\bar{\boldsymbol{\mathrm{z}}}_{t}+[\bar{\boldsymbol{\mathrm{z}}}^{T}_{t}\boldsymbol{\mathrm{Q}}^{(1)}\bar{\boldsymbol{\mathrm{z}}}_{t},...,\bar{\boldsymbol{\mathrm{z}}}^{T}_{t}\boldsymbol{\mathrm{Q}}^{(s)}\bar{\boldsymbol{\mathrm{z}}}_{t}]^{T}$
(4)

with

$\boldsymbol{\mathrm{d}}=\boldsymbol{\mathrm{c}}+\boldsymbol{\mathrm{L}}\boldsymbol{\mathrm{m}}+[{\boldsymbol{\mathrm{m}}}^{T}\boldsymbol{\mathrm{Q}}^{(1)}{\boldsymbol{\mathrm{m}}},...,{\boldsymbol{\mathrm{m}}}^{T}\boldsymbol{\mathrm{Q}}^{(s)}{\boldsymbol{\mathrm{m}}}]^{T}$
(5)

and

$\boldsymbol{\mathrm{A}}=\bigg{(}a_{ij}\bigg{)}=\bigg{(}l_{ij}+\sum^{s}_{k=1}(q_{i,j,k}+q_{i,k,j})m_{k}\bigg{)}$
(6)

The evolution of the fluctuation energy
$K=\frac{1}{2}\sum_{i=1}^{s}\bar{\boldsymbol{\mathrm{z}}}_{i}^{2}$ of the
shifted system is considered as a Lyapunov function [39]. The time derivative
of this quantity can be written as:

$\dot{K}=[\nabla_{\bar{\boldsymbol{\mathrm{z}}}}\boldsymbol{\mathrm{K}}]^{T}\dot{\bar{\boldsymbol{\mathrm{z}}}}=\bar{\boldsymbol{\mathrm{z}}}^{T}\boldsymbol{\mathrm{A}}_{s}\bar{\boldsymbol{\mathrm{z}}}+\boldsymbol{\mathrm{d}}^{T}\bar{\boldsymbol{\mathrm{z}}}$
(7)

where
$\boldsymbol{\mathrm{A}}_{s}=\frac{1}{2}(\boldsymbol{\mathrm{A}}+\boldsymbol{\mathrm{A}}^{T})$.
The contribution of the quadratic terms to the fluctuation energy $K$ is zero
due to the energy-preserving condition (3). A sufficient condition for the
existence of a monotonically attracting trapping region is the existence of a
finite $\boldsymbol{\mathrm{m}}$ such that $\boldsymbol{\mathrm{A}}_{s}$ has
only negative eigenvalues. In the next section, we enforce this condition in
the proposed data-driven scheme for partially-observed systems in order to
train long-term bounded dynamical models.

##  4 Bounded Neural Embedding

This section first briefly outlines the neural embedding of dynamical systems
—NbedDyn—[32], and introduces the bounded neural embedding framework.

###  4.1 Neural embedding of dynamical systems

Let us consider a dynamical system governed by an autonomous ODE:

$\dot{\boldsymbol{\mathrm{z}}}_{t}=f(\boldsymbol{\mathrm{z}}_{t})$ (8)

For most applications, the true state
$\boldsymbol{\mathrm{z}}_{t}\in\mathbb{R}^{s}$ of the system is unknown and
limited to a series of observations $\\{\boldsymbol{\mathrm{x}}_{t}\\}$:

$\boldsymbol{\mathrm{x}}_{t}=\mathcal{H}(\boldsymbol{\mathrm{z}}_{t})$ (9)

where $\mathcal{H}:\mathbb{R}^{s}\rightarrow\mathbb{R}^{n}$ is an observation
operator that does not satisfy the conditions [103] under which the
predictable deterministic dynamics expressed in the space of
$\boldsymbol{\mathrm{z}}$ are still deterministic in the observation space.

The NbedDyn technique tackles this problem by searching an augmented space,
where the states are governed by diffeomorphic flows and can be mapped to the
observations $\boldsymbol{\mathrm{x}}_{t}$. For any given operator
$\mathcal{H}$ of a deterministic dynamical system, Takens theorem [31]
guarantees that such an augmented space exists. However, instead of using a
delay embedding, NbedDyn defines a $d_{E}$-dimensional augmented space with
states $\boldsymbol{\mathrm{u}}_{t}\in\mathbb{R}^{d_{E}}$ as follows:

$\boldsymbol{\mathrm{u}}_{t}^{T}=[\mathcal{M}(\boldsymbol{\mathrm{x}}_{t}^{T}),\boldsymbol{\mathrm{y}}_{t}^{T}]$
(10)

where $\mathcal{M}$ is an order reduction operator such as
$\mathcal{M}(\boldsymbol{\mathrm{x}}_{t})\in\mathbb{R}^{r}$ with $r\leq n$ and
$\boldsymbol{\mathrm{y}}_{t}\in\mathbb{R}^{d_{E}-r}$ are stated as the latent
states. They account for the information of the unobserved components of the
true state $\boldsymbol{\mathrm{z}}_{t}$.

The corresponding dynamics and observation operator are defined as:

$\displaystyle\dot{\boldsymbol{\mathrm{u}}}_{t}$
$\displaystyle=f_{\theta}(\boldsymbol{\mathrm{u_{t}}})$ (11)
$\displaystyle\boldsymbol{\mathrm{x}}_{t}$
$\displaystyle=\mathcal{M}^{-1}(\boldsymbol{\mathrm{G}}\boldsymbol{\mathrm{u_{t}}})$
(12)

where the dynamical operator $f_{\theta}$ belongs to a given family of neural
network operators parameterized by a parameter vector ${\theta}$.
$\boldsymbol{\mathrm{G}}$ is a projection matrix that satisfies
$\mathcal{M}(\boldsymbol{\mathrm{x}}_{t})=\boldsymbol{\mathrm{G}}\boldsymbol{\mathrm{u_{t}}}$.
Using an integration scheme (in all the experiments, the flow of the proposed
framework is built on a classical Runge-Kutta 4 integration scheme), we can
associate $f_{\theta}$ with a diffeomorphic mapping, to derive the flow
$\Phi_{\theta}$ of the NODE model (11) as follows:

$\Phi_{\theta,t}(\boldsymbol{\mathrm{u}}_{t-1})=\boldsymbol{\mathrm{u}}_{t-1}+\int_{t-1}^{t}f_{\theta}(\boldsymbol{\mathrm{u}}_{t-1})$
(13)

From Eqs. (11), (12) and (13), we define a state space model:

$\
\left\\{\begin{array}[]{ccl}\boldsymbol{\mathrm{u}}_{t}&=&\Phi_{\theta,t}(\boldsymbol{\mathrm{u}}_{t-1})\\\
\boldsymbol{\mathrm{x}}_{t}&=&\mathcal{M}^{-1}(\boldsymbol{\mathrm{G}}\boldsymbol{\mathrm{u_{t}}})\end{array}\right.$
(14)

Given an observation time series of size $N+1$
$\\{\boldsymbol{\mathrm{x}}_{0},\ldots,\boldsymbol{\mathrm{x}}_{t_{N}}\\}$,
the Neural Embedding scheme aims to minimize the forecasting error of the
observations with respect to the model parameters and the latent states as
follows:

$\displaystyle\hat{\theta},\boldsymbol{\mathrm{y}}_{t_{0}:t_{N-1}}=\displaystyle\arg\displaystyle\min_{\theta,\\{\boldsymbol{\mathrm{y}}_{t}\\}_{t}}\sum_{i=1}^{N}\|\boldsymbol{\mathrm{x}}_{t_{i}}-\mathcal{M}^{-1}(\boldsymbol{\mathrm{G}}\Phi_{\theta,t_{i}}\left(\boldsymbol{\mathrm{u}}_{t_{i-1}})\right)\|^{2}+\lambda_{1}\|\boldsymbol{\mathrm{u}}_{t_{i}}-\Phi_{\theta,t_{i}}(\boldsymbol{\mathrm{u}}_{t_{i-1}})\|^{2}$
(15)

with $\lambda_{1}$ a weighting parameter and $\|\cdot\|$ the $L^{2}$ norm.

###  4.2 Constrained dynamical embedding

Figure 1:  Sketch of the proposed architecture. Given a sequence of
observations issued from an unknown higher dimensional dynamical system, the
observations are projected into a higher dimensional space parameterized by a
linear quadratic neural ODE model. The parameters of the models, highlighted
in red, as well as the latent states are optimized to minimize the forecasting
of the observations while satisfying the boundedness constraints as proposed
in [39].

When considering the parametrization of the model $f_{\theta}$, a linear-
quadratic formulation is suitable for the identification of reduced order
models of incompressible flows [104]. Furthermore, this architecture makes the
application of the attracting trapping region condition, introduced in the
previous section, tractable using the Schlegel boundedness theorem [39].
Regarding the approximation properties of this class of models, the
combination of Volterra theory, on the solutions of non-linear differential
equations [105], with lifting transformation techniques [106, 107] shows that
linear quadratic models can faithfully represent (in a higher dimension) the
dynamics of non-linear differential equations. In this context, the operator
$f_{\theta}$ in (11) is restricted to be linear quadratic, and can be written,
similarly to (2) as follows:

$\dot{\boldsymbol{\mathrm{u}}}_{t}=\boldsymbol{\mathrm{c}}+\boldsymbol{\mathrm{L}}\boldsymbol{\mathrm{u}}_{t}+[\boldsymbol{\mathrm{u}}^{T}_{t}\boldsymbol{\mathrm{Q}}^{(1)}\boldsymbol{\mathrm{u}}_{t},...,\boldsymbol{\mathrm{u}}^{T}_{t}\boldsymbol{\mathrm{Q}}^{({d_{E}})}\boldsymbol{\mathrm{u}}_{t}]^{T}$
(16)

where $\boldsymbol{\mathrm{c}}\in\mathbb{R}^{d_{E}}$,
$\boldsymbol{\mathrm{L}}\in\mathbb{R}^{{d_{E}}\times{d_{E}}}$ and
$\boldsymbol{\mathrm{Q}}^{(i)}=[q_{i,j,k}]^{d_{E}}_{j,k=1},i=1,...,{d_{E}}$.
The above approximate model is shifted according to
$\bar{\boldsymbol{\mathrm{u}}}=\boldsymbol{\mathrm{u}}-\boldsymbol{\mathrm{m}}$
with $\boldsymbol{\mathrm{m}}\in\mathbb{R}^{d_{E}}$. The approximate dynamical
equation of the shifted state can be written as

$\dot{\bar{\boldsymbol{\mathrm{u}}}}_{t}=\boldsymbol{\mathrm{d}}+\boldsymbol{\mathrm{A}}\bar{\boldsymbol{\mathrm{u}}}_{t}+[\bar{\boldsymbol{\mathrm{u}}}^{T}_{t}\boldsymbol{\mathrm{Q}}^{(1)}\boldsymbol{\mathrm{u}}_{t},...,\bar{\boldsymbol{\mathrm{u}}}^{T}_{t}\boldsymbol{\mathrm{Q}}^{({d_{E}})}\bar{\boldsymbol{\mathrm{u}}}_{t}]^{T}$
(17)

with $\boldsymbol{\mathrm{d}}$ and $\boldsymbol{\mathrm{A}}$ computed
according to (5) and (6) respectively.

The training setting comes to jointly learning the model parameters
$\theta=\\{\boldsymbol{\mathrm{c}},\boldsymbol{\mathrm{L}},\boldsymbol{\mathrm{Q^{1}}},\boldsymbol{\mathrm{Q^{2}}},\cdots,\boldsymbol{\mathrm{Q^{d_{E}}}},\boldsymbol{\mathrm{m}}\\}$
and the latent states $\boldsymbol{\mathrm{y}}$ according to the following
constrained optimization problem

$\displaystyle\begin{aligned}
\hat{\theta},\boldsymbol{\mathrm{y}}_{t_{0}:t_{N-1}}=\displaystyle\arg\min_{\theta,\\{\boldsymbol{\mathrm{y}}_{t}\\}_{t}}&\displaystyle\sum_{i=1}^{N}\|\boldsymbol{\mathrm{x}}_{t_{i}}-\mathcal{M}^{-1}(\boldsymbol{\mathrm{G}}(\Phi_{\theta,t_{i}}\left(\boldsymbol{\mathrm{u}}_{t_{i-1}})+\boldsymbol{\mathrm{m}})\right)\|^{2}\\\
&+\lambda_{1}\|\boldsymbol{\mathrm{u}}_{t_{i}}-\Phi_{\theta,t_{i}}(\boldsymbol{\mathrm{u}}_{t_{i-1}})\|^{2}\\\
~{}\\\ &+\lambda_{2}\mathcal{C}_{1}\\\ ~{}\\\
&+\lambda_{3}\mathcal{C}_{2}\end{aligned}$ (18)

with
$\mathcal{C}_{1}=\sum_{i,j,k=1}^{s}\|q_{i,j,k}+q_{i,k,j}+b_{j,i,k}+b_{j,k,i}+b_{k,i,j}+b_{k,j,i}\|^{2}$
and
$\mathcal{C}_{2}=\sum_{i=1}^{s}\mathrm{Max}(\alpha_{i},0)/\mathrm{Max}(\alpha_{i}+1,0)$
where $\alpha_{i},i=1,...,{d_{E}}$ the eigenvalues of the matrix
$\boldsymbol{\mathrm{A}}_{s}=\frac{1}{2}(\boldsymbol{\mathrm{A}}+\boldsymbol{\mathrm{A}}^{T})$.
The variables $\lambda_{1,2,3}$ are constant weighting parameters. This loss
function corresponds to the initial NbedDyn loss given by (15) with two
additional constraints $\mathcal{C}_{1}$ and $\mathcal{C}_{2}$. The first
constraint $\mathcal{C}_{1}$ steams from the energy-preserving condition given
in (3). It forces the contribution of the quadratic terms of $f_{\theta}$ to
the fluctuation energy to sum up to zero. The second constraint,
$\mathcal{C}_{2}$, ensures that the eigenvalues of
$\boldsymbol{\mathrm{A}}_{s}$ are negative. Satisfying these constraints
guarantees that the model $f_{\theta}$ is bounded through the existence of a
monotonically attracting trapping region that includes the limit-set revealed
by the minimization of the forecasting loss. Overall, Fig. 1 highlights the
main mechanisms of the proposed framework.

##  5 Numerical experiments

Numerical experiments are performed to evaluate and demonstrate the relevance
of the proposed framework for the forecasting and simulation of partially-
observed systems. We consider three case-studies: the Lorenz-63, Lorenz-96,
and Shallow Water Equation (SWE) dynamics. The code associated to this work
can be found at https://github.com/CIA-Oceanix/Bounded-NbedDyn.

###  5.1 Lorenz 63

The Lorenz-63 dynamical system is a 3-dimensional model governed by the
following ODE:

$\left\\{\begin{array}[]{ccl}\frac{dz_{t,1}}{dt}&=&\sigma\left(z_{t,2}-z_{t,2}\right)\\\
\frac{dz_{t,2}}{dt}&=&\rho z_{t,1}-z_{t,2}-z_{t,1}z_{t,3}\\\
\frac{dz_{t,3}}{dt}&=&z_{t,1}z_{t,2}-\beta z_{t,3}\end{array}\right.$ (19)

Under parametrization $\sigma=10$, $\rho=28$ and $\beta=8/3$, this system
exhibits chaotic dynamics with a strange attractor [108].

We simulate a Lorenz-63 state sequence of $5000$ time steps, sampled regularly
at $dt=0.01$. This simulation is computed using the LOSDA ODE solver [109].
The training sequence consists of the initial $4000$ measurements of the first
Lorenz-63 variable, i.e., $\mathcal{H}=[1,0,0]$ and
$\boldsymbol{\mathrm{x}}_{t}=\mathcal{H}(z_{t})=\boldsymbol{\mathrm{z}}_{t,1}$.
The remaining time steps are used as a test set. In this experiment, the
projection operator in (10) is the identity matrix, i.e.
$\mathcal{M(\cdot)}=I_{1}$.

Parameterization of the data-driven models: For benchmarking purposes, we
perform a quantitative comparison with other leading approaches using delay
embedding representations [31]. The parameters of the delay embedding
representation, namely the lag $\tau$ and the dimension $d_{E}$ of the
augmented space, were chosen using standard techniques. Specifically, the lag
parameter was computed using both the mutual information and correlation
techniques [110], respectively denoted as $\tau_{MI}$ and $\tau_{Corr}$.
Regarding the dimension of the embedding representation, we used the Whitney’s
embedding condition $d_{E}=2\mathrm{p}+1$ with $\mathrm{p}$ the dimension of
the hidden limit-set. The delay embedding dimension was also computed using
the false nearest neighbors (FNN) method [111]. We also tested arbitrary
parameters for the delay embedding dimension. Given the delay embedding
representation, we tested the Sparse Regression (SR) method [1] on a second
order polynomial representation of the delay embedding states. The time
integration of the SR model is carried using the LOSDA solver [109]. Regarding
deep learning models, we compare our method to a stacked Bidirectional LSTM
(RNN) and to the latent-ODE model [13]. The RNN model includes 4 LSTM layers
with 10 hidden units. The input sequence of the model is an observation
sequence of size 30. The Latent-ODE model involves a recognition network, a
decoder and a dynamical model. The recognition network is an RNN with 40
hidden units. We use a 3-dimensional latent space that evolves in time
according to a linear quadratic neural ODE model. The decoder network is a
fully-connected neural network with one hidden layer and 15 hidden units.

The proposed framework is tested with a dimension of the augmented state space
$d_{E}=3$, i.e., two latent states are concatenated to the observed variable.
Overall the augmented state $\boldsymbol{\mathrm{u}}$ is constructed as
$\boldsymbol{\mathrm{u}}_{t}=[\boldsymbol{\mathrm{z}}_{1,t},\boldsymbol{\mathrm{y}}_{1,t},\boldsymbol{\mathrm{y}}_{2,t}]^{T}$.
The parameters of the data-driven model governing $\boldsymbol{\mathrm{u}}$,
as well as the latent states $\boldsymbol{\mathrm{y}}_{1,t}$ and
$\boldsymbol{\mathrm{y}}_{2,t}$ are optimized based on both i) the initial
NbedDyn formulation presented as proposed in [32] (i.e., using the
optimization criterion given by (15)) and ii) the constrained version
introduced in this work, illustrated by criterion (18).

Model $t_{0}+dt$ $t_{0}+4dt$ $\lambda_{1}$ $\lambda_{1}$ SR $\tau_{MI}=16$,
$d_{E}(FNN)=3$ $2.8E-2$ $1.4E-1$ $NaN$ $NaN$ $\tau_{MI}=16$,
$d_{E}$(Whitney)$=6$ $1.3E-1$ $3.0E-1$ $NaN$ $NaN$ $\tau_{Corr}=27$,
$d_{E}$(FNN)$=3$ $6.6E-1$ $2.9$ $NaN$ $NaN$ $\tau_{Corr}=27$,
$d_{E}$(Whitney)$=6$ $1.9E-1$ $9.0E-1$ $NaN$ $NaN$ $\tau=6$, $d_{E}=3$
$4.0E-3$ $0.02.7E-2$ $NaN$ $NaN$ $\tau=10$, $d_{E}=3$ $0.01.3E-2$ $0.07.3E-2$
$NaN$ $NaN$ Latent-ODE $6.5E-2\pm 4.4E-2$ $1.5E-1\pm 1.1E-1$ $NaN$ $NaN$ RNN
$2.7E-1\pm 5.7E-2$ $3.1E-1\pm 1.4E-1$ $-5.4\pm 0.0$ NbedDyn $d_{E}=3$ 8.8E-5
$\pm$ 2.9E-5 5.8E-4 $\pm$ 1.2E-4 0.88 $\pm$ 0.02 $NaN$ Constrained NbedDyn
$d_{E}=3$ 1.2E-4 $\pm$ 4.1 E-5 3.5E-4 $\pm$ 7.4 E-5 0.91 $\pm$ 0.02 0.90 $\pm$
0.02

Table 1: Forecasting performance on the test set of the benchmarked data-
driven models for Lorenz-63 case-study: first two columns : root mean square
error (RMSE) for different forecasting time-steps, third column : largest
Lyapunov exponent of a predicted series of length of 10000 time-steps (The
true largest Lyapunov exponent of the Lorenz 63 model is 0.91 [112]).

Forecasting performance of the proposed data-driven models: We evaluate in
Table 1 the short-term forecast as well as the topological structure of the
simulated limit-sets (illustrated for instance through the largest Lyapunov
exponent). Regarding the short-term forecast, both versions of the NbedDyn
model outperform benchmarked state-of-the-art techniques in terms of root mean
square error (RMSE). The constrained version results in a slightly larger RMSE
score due to the additional boundedness terms in the optimization criterion.
In the other hand, when considering the long term simulation performance of
the tested data-driven models, when the initial condition is inside the
spanned attractor of the augmented states, the dynamical model optimized using
the criterion (15) gives trajectories that are bounded and with similar
topological characteristics to the true Lorenz 63 model. However, when the
initial condition is far from the spanned attractor, the model optimized
through (15) diverges to infinity. From a machine learning perspective, this
is the direct consequence of a poor generalization performance to states that
are far from the attractor spanned by the training data. In this situation,
our model involves several attracting regions of chaotic and unstable
solutions. When the initial condition is far from the spanned attractor, the
state evolution is dominated by a positive energy growth, which makes our
model diverge to infinity. The constrained model in the other hand, satisfies
elementary conservation constraints that are present in the actual Lorenz 63
system and leads to a bounded behavior with a large attracting region of the
chaotic limit-set. We further illustrate these conclusions with forecasting
examples depicted in Figure 2 with initial conditions lying both inside and
outside the attractor.

Figure 2: Forecasting performance of the NbedDyn models for Lorenz 63 case-
study. First row, NbedDyn model as presented in [32]; second row, proposed
constrained NbedDyn model; third row, true Lorenz 63 model. The left (resp.
right) panel illustrates the forecasting performance of the NbedDyn for an
initial condition inside (resp. outside) the attractor. It is worth noting
that the scales of the forecasted latent states $\boldsymbol{\mathrm{y}}_{1}$
and $\boldsymbol{\mathrm{y}}_{2}$ are different from the ones of the unseen
states $\boldsymbol{\mathrm{z}}_{1}$ and $\boldsymbol{\mathrm{z}}_{2}$. This
stems from the fact that given partial observations of the system, several
models can reflect the variability of the observations while being
diffeomorphic to the actual governing dynamics.

###  5.2 Lorenz 96

The Lorenz-96 system [113] with periodic boundary conditions is governed by
the following ODE:

$\frac{d\boldsymbol{\mathrm{z}}_{t,i}}{dt}=(\boldsymbol{\mathrm{z}}_{t,i+1}-\boldsymbol{\mathrm{z}}_{t,i-2})\boldsymbol{\mathrm{z}}_{t,i-1}-\boldsymbol{\mathrm{z}}_{t,i}+F$
(20)

with $\boldsymbol{\mathrm{z}}_{t,-1}=\boldsymbol{\mathrm{z}}_{t,s}$,
$\boldsymbol{\mathrm{z}}_{t,s+1}=\boldsymbol{\mathrm{z}}_{t,1}$.

The Lorenz-96 state sequences is also simulated using the LOSDA ODE solver
[109] with $F=8$ and a dimension $s=40$. Similar to the previous experiment,
the numerical simulation consists of $5000$ time steps sampled at $dt=0.01$.
The initial $4000$ measurements of the first 20 states, i.e.,
$\boldsymbol{\mathrm{x}}_{t}=\mathcal{H}(z_{t})=[\boldsymbol{\mathrm{z}}_{1,t},\cdots,\boldsymbol{\mathrm{z}}_{20,t}$]T
are used as training data. The remaining time steps are used as a test set. In
this experiment, the projection operator in (10) is the identity matrix, i.e.
$\mathcal{M(\cdot)}=I_{20}$.

Parametrization of the data-driven models : The proposed framework, is tested
with a dimension of the augmented state space $d_{E}=40$ i.e., 20 latent
states are concatenated to the observed variable. In this experiment, the
quadratic part of the NODE model $f_{\theta}$ is convolutional and corresponds
to the true non-linear interactions in (20) with trainable weights. Such
representation significantly accelerates the computation of the constraint
$\mathcal{C}_{1}$ in (18) 111Other experiments, not shown here, demonstrate
that a full linear quadratic NODE model also converges (but at a slower rate)
to solutions with the same properties as the convolutional model evaluated in
this experiment.. Similar to the Lorenz 63 experiment, both the constrained
and unconstrained versions of the model are compared to the sparse regression
technique, to a stacked bidirectional LSTM (RNN) and to the Latent-ODE model
[13]). The SR model is based on a second-order polynomial representation. The
state-space of the model is built on a stacked delay embedding of the 20
observed states. For each observed state, the dimension of the embedding is
computed using the FNN method and the lag is chosen according to both the
mutual information and the correlation techniques. The time integration of the
SR model is carried using the LOSDA solver [109]. The RNN model includes 10
LSTM layers with 100 hidden units. The input sequence of the model is an
observation sequence of size 20. Regarding the Latent-ODE model, the
recognition network is an RNN with 100 hidden units. The dimension of the
latent space is set to 40 and the dynamical model is a linear quadratic neural
ODE. The decoder network is a fully-connected neural network with one hidden
layer and 100 hidden units.

Model $t_{0}+dt$ $t_{0}+4dt$ $\lambda_{1}$ $\lambda_{1}$ SR( $\tau_{MI}$ )
$>10$ $>10$ $NaN$ $NaN$ SR( $\tau_{Corr}$ ) $>10$ $>10$ $NaN$ $NaN$ Latent-ODE
0.262 $\pm$ 2.08E-1 0.560 $\pm$ 3.7E-1 NaN NaN RNN 0.186 $\pm$ 3.6E-2 0.231
$\pm$ 5.7E-2 5.091 $\pm$ 1.856 2.451 $\pm$ 2.092 NbedDyn 0.009 $\pm$ 2.4E-3
0.293 $\pm$ 6.4E-2 NaN NaN Constrained NbedDyn 0.012 $\pm$ 2.3E-3 0.036 $\pm$
6.9E-3 1.211 $\pm$ 1.32E-1 1.399 $\pm$ 6.6E-2

Table 2: Forecasting performance on the test set of the data-driven models for
Lorenz-96 dynamics where only the first 20 state variables are observed: first
two columns : mean RMSE for different forecasting time-steps, third column :
largest Lyapunov exponent of a predicted series of length of 10000 time-steps
(The true largest Lyapunov exponent of the Lorenz 96 model is 1.67 [54]).

Forecasting performance of the proposed data-driven models: Similar to the
Lorenz 63 experiment, both versions of the NbedDyn model outperform
benchmarked state-of-the-art techniques in terms of short-term forecast RMSE.
When considering the attractor reconstruction based on long-term simulation of
the data-driven models, the constrained version of the NbedDyn model is the
only model able to unfold the Lorenz 96 attractor. The other data-driven
models, including the unconstrained NbedDyn either diverge to infinity or
generate trajectories that do not match the Lorenz 96 hidden attractor
(illustrated through the largest Lyapunov exponent).

Figure 3:  Mean Fast Fourier Transform (FFT) modulus of the benchmarked data-
driven models for Lorenz-96 case-study. We report the FFT modulus of simulated
time series over the test period for the benchmarked data-driven models. The
solid lines represent the FFT averaged over five training runs for each
scheme. The light-color interval is the standard deviation of the different
models. The unconstrained NbedDyn trajectories diverge after a short
forecasting time, thus, its FFT is omitted.  Figure 4:  Forecasting
performance on the benchmarked data-driven models for Lorenz-96 case-study.
First row: ground-truth Lorenz 96 time series; Second row : data-driven model
simulations over a test period; Third row: associated RMSE map. For each
subplot, the vertical axis refers to the time axis and the horizontal axis to
the 20-dimensional observation state. As unconstrained NbedDyn trajectories
blowup after a short time horizon, no observation state can be depicted.

We further illustrate the quantitative analysis conclusions through the visual
comparison of the FFT modulus as well as a forecasting example in figures 3
and 4 respectively. The proposed architecture shows a better match to the true
trajectory both in the temporal and spectral domains. Specifically, the
constrained NbedDyn trajectory, although diverges from the ground truth (due
to the chaotic nature of the attractor), keeps a similar spatio-temporal
behavior illustrated by the prediction example in Fig. 4 as well as the
spectrum in Fig. 3. Interestingly, this performance is truly due to the
boundedness constraints since the unconstrained version, and even using the
same ODE parameterization, diverges to infinity after a short simulation time.
The RNN trajectories decently capture the time frequency patterns of the true
states but fail to reproduce the Lorenz 96 dynamics in the phase space.

###  5.3 Patch shallow Water equation (PSWE)

As a third case-study, we consider SWE dynamics governed by the following set
of equations:

$\left\\{\begin{aligned} &\frac{\partial v_{x}}{\partial
t}-F_{v_{x}}=-g\frac{\partial\eta}{\partial x}\\\ &\frac{\partial
v_{y}}{\partial t}-F_{v_{y}}=-g\frac{\partial\eta}{\partial y}\\\
&\frac{\partial\eta}{\partial t}+\frac{\partial(\eta+H)v_{x}}{\partial
x}+\frac{\partial(\eta+H)v_{y}}{\partial y}&=0\end{aligned}\right.$ (21)

where $x$ and $y$ represent the 2D directions of the fluid. $\eta$ is the
fluid surface elevation, and $(v_{x},v_{y})$ is the fluid’s horizontal flow
velocity. g is the acceleration due to gravity, which is taken to be equal to
$9.81m/s^{2}$. $H$ corresponds to the total depth (here $H=100m$) of the fluid
and $f=f_{0}+\lambda y$ is the Coriolis forcing where $f_{0}=1E-4s^{-1}$ and
$\lambda=2E^{-11}(m\times s)^{-1}$.

The direct numerical simulation of (21) is carried using a finite difference
method. The length of the domain is set to $1000km\times 1000km$ with a
corresponding regular discretization of $80\times 80$. The temporal step size
$dt$ is set to satisfy the Courant–Friedrichs–Lewy condition ($dt=40.41$
seconds). As observed process, we consider a low-resolution version of a patch
of size $250km\times 250km$ of the sea surface elevation $\eta$. The patch is
located in the center of the 2D domain. The low-resolution projection is
computed from an 8-dimensional EOF decomposition, which amounts to capture
$80\%$ of the total variance. Regarding the training configuration, $100000$
simulation time steps are generated. The transient (first $2500$) time steps
are omitted, and we use the post-transient first 49701 time-steps as training
data. The remaining sequence is used as a test set.

This case study combines numerous layers of complexity when compared to the
previous experiments. Specifically, the dynamics of the considered
$250km\times 250km$ domain depend on the neighboring regions. Its variability
is rigorously governed by the Shallow Water Equation in the full domain. As
such, one cannot derive an autonomous differential system for the subdomain.
Besides, through the EOF-based filtering, the observed process only refers to
coarse-scale dynamics. Reduced-order modeling schemes [114] and learning-based
closure models [115] would be typical solutions to address this issue, which
would rely on the approximation of the observed dynamics by an autonomous
system in the observation space or in a lower-dimensional projection of this
space. We show in this experiment that the proposed framework can account for
unknown boundary conditions and fine-scale dynamics by retrieving the most
appropriate missing information to forecast future coarse-scale observations
within the considered patch. We may remind that the learning of the latent
states only relies on the observed variables.

Parametrization of the data-driven models : All the tested models are built on
the 8-dimensional EOF basis used to filter the observations. In the proposed
framework, this setting corresponds to a specific parameterisation of the
projection operator $\mathcal{M(\cdot)}$ in (10) where
$\mathcal{M}(\boldsymbol{\mathrm{x}}_{t})\in\mathbb{R}^{8}$ is an EOF
projection. The NbedDyn model (both the constrained and unconstrained
settings) is constructed with a dimension of the augmented state $d_{E}=18$.
Our framework is compared to a stacked Bidirectional LSTM (RNN) and to a
linear regression model. The RNN includes 10 LSTM layers with 100 hidden
units. The input sequence of the RNN is an observation sequence of size 40.
The linear model, that we call here Delay EDMD, is built on a Singular Value
Decomposition (SVD) of delay embedding coordinates of the observations [40].
The delay embedding is computed, for every EOF component, using a lag
embedding equal to one time step and a dimension equal to 200. The dimension
of the SVD is set to 150, which accounts for over 99.99% of the total variance
of the delay embedding representation.

Model | $t_{0}+dt$ | $t_{0}+4dt$ |  |
---|---|---|---|---
RNN |  |  | 7.35E-4 $\pm$ 2.55E-4 |  1.28E-3 $\pm$ 3.71E-4 |  |
NbedDyn |  |  | 1.33E-3 $\pm$ 7.88E-5 | 3.66E-3 $\pm$ 1.92E-4 |  |
Delay EDMD |  |  | 1.43E-3 | 1.42E-3 |  |
Constrained NbedDyn |  |  |  6.28E-4 $\pm$ 1.69E-4 | 1.76E-3 $\pm$ 4.7E-04 |  |
Table 3: Forecasting performance of the benchmarked data-driven models for
PSWE case-study: we report the mean RMSE mertrics for different forecasting
time-steps for the tesst set.

Forecasting performance of the proposed data-driven models: regarding the
short-term forecasting performance, reported for instance in table 3, the
benchmarked models achieve similar errors. These results are also highlighted
in the forecasting examples given in figure 5 where the short term forecast,
highlighted in the blue panel, are similar to the true state, up to 27 minutes
where only the Delay EDMD technique is still qualitatively close to the ground
truth. The analysis of the long-term simulated states, bordered by a red panel
in Figure 5, gives in the other hand a different conclusion. While the
proposed constrained NbedDyn model keeps simulating realistic states, both the
unconstrained NbedDyn and RNN models get stuck either at equilibrium points,
or at very slow manifolds. The linear Delay EDMD model converges to the origin
after a finite amount of time, which is expected since finite-dimensional
linear models can not reflect a chaotic behavior. We can draw the same
analysis from the mean FFT modulus of the forecasted EOF modes, as well as the
mean radially averaged PSD of the 2D fields in Fig. 6. When considering the
spatial spectrum, only the proposed framework is able to match perfectly the
spectrum of the ground-truth. This result illustrates the ability of our model
to capture all the spatial scales of the dynamics and to keep simulating these
scales after long simulation times. Regarding the temporal spectrum, the
constrained NbedDyn and the Delay EDMD models match the spectrum of the ground
truth both in the medium and high frequency range. The low-frequency range in
the other hand is not well captured by the model. This suggests that we may
need to consider an optimization objective (18) with longer prediction times
than the one used in this work to embed realistic low-frequency patterns.

Figure 5: Prediction example over the test period of the data-driven models
for the PSWE case-study. The blue panel illustrates a short-term forecast for
the benchmarked models and the red panel highlights the associated long-term
simulations. A spectral analysis of the simulated trajectories, both in space
and time, are given in figure 6.

(a)

(b)

Figure 6:  Spatio-temporal spectrums of the data-driven models with respect to
the ground truth. The FFT modulus of the data-driven models in (a) are
computed from a simulated time series in the test set. The radially-averaged
Power Spectral Density (PSD) in (b) is averaged over a simulation of size
10000 time steps .The lines represent the averaged spectrums, computed over
five optimized models. The light-color intervals refer to the standard
deviation.

##  6 Conclusion

Two central challenges in data-driven modeling involve how to handle partial
observations and how to guarantee stability in the resulting models. In this
paper, we demonstrate that combining an augmented neural ODE formulation and
boundedness constraints, it is possible to considerably improve the data-
driven identification and forecasting of partially-observed dynamics. The
proposed framework is assessed on several case studies with respect to state-
of-the-art forecasting and system identification techniques. Specifically, for
classical chaotic ODEs (Lorenz 63 and 96), the proposed methodology matches
state-of-the-art short-term forecasting performances while ensuring realistic
long-term patterns for partially-observed settings. We also report a more
complex case study, in which we have considered low-resolution observations of
a subdomain of SWE dynamics with unknown boundary conditions. The results from
this study support the potential of the proposed technique to derive faithful
data-based representations of ocean flows, especially ocean surface dynamics.
Such dynamics involve, similar to the PSWE case study, but at a higher level
of complexity, unseen components, missing forcing, unknown boundary conditions
and unresolved spatio-temporal scales. In all numerical experiments, both the
data-driven model formulation in an augmented space and the boundedness
constraints are key features to capture the dynamics underlying the
observations.

An avenue of promising future work will involve extending the proposed
framework to other non-linear models through the generalization of the
boundedness constraints. Foremost, polynomial models could be considered by
imposing, similar to the quadratic terms of the proposed architecture, an
energy-preserving non-linearity. Other non-linear architectures may still need
further investigation, since the fluctuation energy may not be the appropriate
Lyapunov function to consider. The recent works of [37] may be useful for
finding Lyapunov functions that guarantee stability and boundedness.

To recall, explicitly constraining chaos within a data-driven formulation is
rigorously intractable since such a behavior is only revealed, and thus
potentially constrained using long-term simulations. From this point of view,
the proposed framework is not guaranteed to reach an expected chaotic
evolution, since this is not explicitly constrained within the framework.
Specifically, when considering the bounded NbedDyn model, the learning
criterion allows for any bounded limit-set, including stable ones, as long as
the short-term forecast of the observations is minimized. This criterion does
not guarantee a replication of the chaotic observations, and long-term
simulations of the model can lead to undesirable stable limit-sets.
Interestingly, stable limit-sets can be fully characterized by the dynamical
equation, without resorting to brute force long-term simulation of the model.
Thus, these limit sets can be propagated to infinity. However, sending stable
limit-sets to infinity may result in new questions about how observation data
may help to regularize the learning stage.

## Acknowledgements

This work was supported by LEFE program (LEFE MANU project IA-OAC), CNES
(grant OSTST DUACS-HR) and ANR Projects Melody and OceaniX. It benefited from
HPC and GPU resources from Azure (Microsoft Azure grant) and from GENCI-IDRIS
(Grant 2020-101030). SLB acknowledges valuable discussions with Alan
Kaptanoglu.

## References

  * [1] Steven L. Brunton, Joshua L. Proctor, and J. Nathan Kutz.  Discovering governing equations from data by sparse identification of nonlinear dynamical systems.  Proceedings of the National Academy of Sciences, 113(15):3932–3937, April 2016.
  * [2] Joseph Bakarji, Kathleen Champion, J Nathan Kutz, and Steven L Brunton.  Discovering governing equations from partial measurements with deep delay autoencoders.  arXiv preprint arXiv:2201.05136, 2022.
  * [3] Jean-Christophe Loiseau and Steven L Brunton.  Constrained sparse galerkin regression.  Journal of Fluid Mechanics, 838:42–67, 2018.
  * [4] Alan A Kaptanoglu, Jared L Callaham, Aleksandr Aravkin, Christopher J Hansen, and Steven L Brunton.  Promoting global stability in data-driven models of quadratic nonlinear dynamics.  Physical Review Fluids, 6(9):094401, 2021.
  * [5] Boris Kramer.  Stability domains for quadratic-bilinear reduced-order models.  SIAM Journal on Applied Dynamical Systems, 20(2):981–996, 2021\.
  * [6] Ruoyu Sun.  Optimization for deep learning: theory and algorithms.  arXiv preprint arXiv:1912.08957, 2019.
  * [7] Yann LeCun, Yoshua Bengio, and Geoffrey Hinton.  Deep learning.  nature, 521(7553):436–444, 2015.
  * [8] M Grewal and Keith Glover.  Identifiability of linear and nonlinear dynamical systems.  IEEE Transactions on automatic control, 21(6):833–837, 1976.
  * [9] Keith Glover and Jan Willems.  Parametrizations of linear dynamical systems: Canonical forms and identifiability.  IEEE Transactions on Automatic Control, 19(6):640–646, 1974.
  * [10] Bo Wahlberg.  On the identification of continuous time dynamical systems.  IFAC Proceedings Volumes, 21(9):435–440, 1988.
  * [11] S Narendra Kumpati, Parthasarathy Kannan, et al.  Identification and control of dynamical systems using neural networks.  IEEE Transactions on neural networks, 1(1):4–27, 1990.
  * [12] Jaideep Pathak, Brian Hunt, Michelle Girvan, Zhixin Lu, and Edward Ott.  Model-free prediction of large spatiotemporally chaotic systems from data: A reservoir computing approach.  Phys. Rev. Lett., 120:024102, Jan 2018.
  * [13] Tian Qi Chen, Yulia Rubanova, Jesse Bettencourt, and David K Duvenaud.  Neural ordinary differential equations.  In Advances in Neural Information Processing Systems, pages 6571–6583, 2018.
  * [14] Olalekan Ogunmolu, Xuejun Gu, Steve Jiang, and Nicholas Gans.  Nonlinear systems identification using deep dynamic neural networks.  arXiv preprint arXiv:1610.01439, 2016.
  * [15] Ronan Fablet, Said Ouala, and Cedric Herzet.  Bilinear residual Neural Network for the identification and forecasting of dynamical systems.  SciRate, dec 2017.
  * [16] Duong Nguyen, Said Ouala, Lucas Drumetz, and Ronan Fablet.  Assimilation-based learning of chaotic dynamical systems from noisy and partial data.  In ICASSP 2020 - 2020 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), pages 3862–3866, 2020.
  * [17] Duong Nguyen, Said Ouala, Lucas Drumetz, and Ronan Fablet.  Variational deep learning for the identification and reconstruction of chaotic and stochastic dynamical systems from noisy and partial observations.  arXiv preprint arXiv:2009.02296, 2020.
  * [18] Said Ouala, Laurent Debreu, Ananda Pascual, Bertrand Chapron, Fabrice Collard, Lucile Gaultier, and Ronan Fablet.  Learning runge-kutta integration schemes for ode simulation and identification.  arXiv preprint arXiv:2105.04999, 2021.
  * [19] Anne Braakmann-Folgmann, Ribana Roscher, Susanne Wenzel, Bernd Uebbing, and Jürgen Kusche.  Sea Level Anomaly Prediction using Recurrent Neural Networks.  arXiv:1710.07099 [cs], oct 2017.  arXiv: 1710.07099.
  * [20] Bryan Lim and Stefan Zohren.  Time series forecasting with deep learning: A survey, 2020.
  * [21] Said Ouala, Ronan Fablet, Cédric Herzet, Bertrand Chapron, Ananda Pascual, Fabrice Collard, and Lucile Gaultier.  Neural network based kalman filters for the spatio-temporal interpolation of satellite-derived sea surface temperature.  Remote Sens., 10(12):1864, Nov 2018.
  * [22] William Gilpin.  Deep reconstruction of strange attractors from time series, 2020.
  * [23] Ronan Fablet, Bertrand Chapron, Lucas Drumetz, Etienne Mémin, Olivier Pannekoucke, and François Rousseau.  Learning variational data assimilation models and solvers.  arXiv preprint arXiv:2007.12941, 2020.
  * [24] Kangbeom Cheon, Jaehoon Kim, Moussa Hamadache, and Dongik Lee.  On replacing pid controller with deep learning controller for dc motor system.  Journal of Automation and Control Engineering Vol, 3(6), 2015.
  * [25] Eurika Kaiser, J Nathan Kutz, and Steven L Brunton.  Sparse identification of nonlinear dynamics for model predictive control in the low-data limit.  Proceedings of the Royal Society A, 474(2219):20180335, 2018.
  * [26] Omer San and Romit Maulik.  Extreme learning machine for reduced order modeling of turbulent geophysical flows.  Physical Review E, 97(4):042322, 2018.
  * [27] Abhinav Gupta and Pierre FJ Lermusiaux.  Neural closure models for dynamical systems.  Proceedings of the Royal Society A, 477(2252):20201004, 2021.
  * [28] Said Ouala, Cedric Herzet, and Ronan Fablet.  Sea surface temperature prediction and reconstruction using patch-level neural network representations.  arXiv:1806.00144 [cs, stat], may 2018.  arXiv: 1806.00144.
  * [29] Ronan Fablet, Lucas Drumetz, and Francois Rousseau.  Joint learning of variational representations and solvers for inverse problems with partially-observed data, 2020.
  * [30] Suman Ravuri, Karel Lenc, Matthew Willson, Dmitry Kangin, Remi Lam, Piotr Mirowski, Megan Fitzsimons, Maria Athanassiadou, Sheleem Kashem, Sam Madge, et al.  Skilful precipitation nowcasting using deep generative models of radar.  Nature, 597(7878):672–677, 2021.
  * [31] Floris Takens.  Detecting strange attractors in turbulence.  In David Rand and Lai-Sang Young, editors, Dynamical Systems and Turbulence, Warwick 1980, pages 366–381, Berlin, Heidelberg, 1981. Springer Berlin Heidelberg.
  * [32] S Ouala, D Nguyen, L Drumetz, B Chapron, A Pascual, F Collard, L Gaultier, and R Fablet.  Learning latent dynamics for partially observed chaotic systems.  Chaos: An Interdisciplinary Journal of Nonlinear Science, 30(10):103121, 2020.
  * [33] Alan A Kaptanoglu, Kyle D Morgan, Chris J Hansen, and Steven L Brunton.  Physics-constrained, low-dimensional models for mhd: First-principles and data-driven approaches.  arXiv preprint arXiv:2004.10389, 2020.
  * [34] Yifei Guan, Steven L Brunton, and Igor Novosselov.  Sparse nonlinear models of chaotic electroconvection.  arXiv preprint arXiv:2009.11862, 2020.
  * [35] Rui Wang, Robin Walters, and Rose Yu.  Incorporating symmetry into deep dynamics models for improved generalization.  arXiv preprint arXiv:2002.03061, 2020.
  * [36] Rui Wang, Karthik Kashinath, Mustafa Mustafa, Adrian Albert, and Rose Yu.  Towards physics-informed deep learning for turbulent flow prediction.  In Proceedings of the 26th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining, pages 1457–1466, 2020.
  * [37] Gaurav Manek and J. Zico Kolter.  Learning stable deep dynamics models, 2020.
  * [38] Ivan Dario Jimenez Rodriguez, Aaron D. Ames, and Yisong Yue.  Lyanet: A lyapunov framework for training neural odes, 2022.
  * [39] Michael Schlegel and Bernd R. Noack.  On long-term boundedness of galerkin models.  Journal of Fluid Mechanics, 765:325–352, 2015.
  * [40] Mason Kamb, Eurika Kaiser, Steven L Brunton, and J Nathan Kutz.  Time-delay observables for koopman: Theory and applications.  SIAM Journal on Applied Dynamical Systems, 19(2):886–917, 2020\.
  * [41] Michael Schmidt and Hod Lipson.  Distilling free-form natural laws from experimental data.  Science, 324(5923):81–85, 2009.
  * [42] Johan Paduart, Lieve Lauwers, Jan Swevers, Kris Smolders, Johan Schoukens, and Rik Pintelon.  Identification of nonlinear systems using Polynomial Nonlinear State Space models.  Automatica, 46(4):647–656, April 2010.
  * [43] Wen-Xu Wang, Rui Yang, Ying-Cheng Lai, Vassilios Kovanis, and Celso Grebogi.  Predicting catastrophes in nonlinear dynamical systems by compressive sensing.  Phys. Rev. Lett., 106:154101, Apr 2011.
  * [44] Ye Yuan, Xiuchuan Tang, Wei Zhou, Wei Pan, Xiuting Li, Hai-Tao Zhang, Han Ding, and Jorge Goncalves.  Data driven discovery of cyber physical systems.  Nature communications, 10(1):1–9, 2019.
  * [45] Steffen Wiewel, Moritz Becher, and Nils Thuerey.  Latent-space physics: Towards learning the temporal evolution of fluid flow, 2018.
  * [46] Maziar Raissi, Paris Perdikaris, and George E Karniadakis.  Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations.  Journal of Computational Physics, 378:686–707, 2019.
  * [47] S. Ouala, A. Pascual, and R. Fablet.  Residual integration neural network.  In ICASSP 2019 - 2019 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), pages 3622–3626, May 2019\.
  * [48] Guorui Shen, Jürgen Kurths, and Ye Yuan.  Sequence-to-sequence prediction of spatiotemporal systems.  Chaos: An Interdisciplinary Journal of Nonlinear Science, 30(2):023102, 2020.
  * [49] Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun.  Deep Residual Learning for Image Recognition.  arXiv:1512.03385 [cs], december 2015.  arXiv: 1512.03385.
  * [50] Marco Fraccaro, Søren Kaae Sønderby, Ulrich Paquet, and Ole Winther.  Sequential Neural Models with Stochastic Layers.  arXiv:1605.07571 [cs, stat], may 2016.  arXiv: 1605.07571.
  * [51] Rahul G. Krishnan, Uri Shalit, and David Sontag.  Structured Inference Networks for Nonlinear State Space Models.  arXiv:1609.09869 [cs, stat], September 2016.  arXiv: 1609.09869.
  * [52] Geir Evensen.  Data assimilation: the ensemble Kalman filter.  Springer Science & Business Media, 2009.
  * [53] Marc Bocquet, Julien Brajard, Alberto Carrassi, and Laurent Bertino.  Data assimilation as a learning tool to infer ordinary differential equation representations of dynamical models.  Nonlinear Processes in Geophysics, 26(3):143–162, 2019.
  * [54] Julien Brajard, Alberto Carassi, Marc Bocquet, and Laurent Bertino.  Combining data assimilation and machine learning to emulate a dynamical model from sparse and noisy observations: a case study with the lorenz 96 model.  arXiv preprint arXiv:2001.01520, 2020.
  * [55] D Nagesh Kumar, K Srinivasa Raju, and T Sathish.  River flow forecasting using recurrent neural networks.  Water resources management, 18(2):143–161, 2004.
  * [56] Sumit Kumar, Lasani Hussain, Sekhar Banarjee, and Motahar Reza.  Energy load forecasting using deep learning approach-lstm and gru in spark cluster.  In 2018 Fifth International Conference on Emerging Applications of Information Technology (EAIT), pages 1–4. IEEE, 2018.
  * [57] Yi-Ting Tsai, Yu-Ren Zeng, and Yue-Shan Chang.  Air pollution forecasting using rnn with lstm.  In 2018 IEEE 16th Intl Conf on Dependable, Autonomic and Secure Computing, 16th Intl Conf on Pervasive Intelligence and Computing, 4th Intl Conf on Big Data Intelligence and Computing and Cyber Science and Technology Congress (DASC/PiCom/DataCom/CyberSciTech), pages 1074–1079. IEEE, 2018.
  * [58] B. O. Koopman.  Hamiltonian systems and transformations in hilbert space.  Proceedings of the National Academy of Sciences of the United States of America, 17(5):315–318, 1931.
  * [59] Matthew O Williams, Ioannis G Kevrekidis, and Clarence W Rowley.  A data–driven approximation of the koopman operator: Extending dynamic mode decomposition.  Journal of Nonlinear Science, 25(6):1307–1346, 2015.
  * [60] Steven L Brunton, Bingni W Brunton, Joshua L Proctor, and J Nathan Kutz.  Koopman invariant subspaces and finite linear representations of nonlinear dynamical systems for control.  PloS one, 11(2):e0150171, 2016.
  * [61] Steven L Brunton, Bingni W Brunton, Joshua L Proctor, Eurika Kaiser, and J Nathan Kutz.  Chaos as an intermittently forced linear system.  Nature communications, 8(1):1–9, 2017.
  * [62] Naoya Takeishi, Yoshinobu Kawahara, and Takehisa Yairi.  Learning koopman invariant subspaces for dynamic mode decomposition.  arXiv preprint arXiv:1710.04340, 2017.
  * [63] Bethany Lusch, J Nathan Kutz, and Steven L Brunton.  Deep learning for universal linear embeddings of nonlinear dynamics.  Nature communications, 9(1):1–10, 2018.
  * [64] Enoch Yeung, Soumya Kundu, and Nathan Hodas.  Learning deep neural network representations for koopman operators of nonlinear dynamical systems.  In 2019 American Control Conference (ACC), pages 4832–4839. IEEE, 2019.
  * [65] Henning Lange, Steven L Brunton, and Nathan Kutz.  From fourier to koopman: Spectral methods for long-term time series prediction.  arXiv preprint arXiv:2004.00574, 2020.
  * [66] Omri Azencot, N Benjamin Erichson, Vanessa Lin, and Michael Mahoney.  Forecasting sequential data using consistent koopman autoencoders.  In International Conference on Machine Learning, pages 475–485. PMLR, 2020.
  * [67] Julian Rice, Wenwei Xu, and Andrew August.  Analyzing koopman approaches to physics-informed machine learning for long-term sea-surface temperature forecasting.  arXiv preprint arXiv:2010.00399, 2020.
  * [68] Alexandre Mauroy, Igor Mezić, and Yoshihiko Susuki.  The Koopman Operator in Systems and Control: Concepts, Methodologies, and Applications, volume 484.  Springer Nature, 2020.
  * [69] Andrew J Majda and John Harlim.  Physics constrained nonlinear regression models for time series.  Nonlinearity, 26(1):201, 2012.
  * [70] Andreas Mardt, Luca Pasquali, Frank Noé, and Hao Wu.  Deep learning markov and koopman models with physical constraints.  In Mathematical and Scientific Machine Learning, pages 451–475. PMLR, 2020.
  * [71] Hao Ma, Xiangyu Hu, Yuxuan Zhang, Nils Thuerey, and Oskar J. Haidn.  A combined data-driven and physics-driven method for steady heat conduction prediction using deep convolutional neural networks, 2020.
  * [72] Xiaowei Jin, Shengze Cai, Hui Li, and George Em Karniadakis.  Nsfnets (navier-stokes flow nets): Physics-informed neural networks for the incompressible navier-stokes equations.  Journal of Computational Physics, 426:109951, 2021.
  * [73] Shengze Cai, Zhicheng Wang, Sifan Wang, Paris Perdikaris, and George Karniadakis.  Physics-informed neural networks (pinns) for heat transfer problems.  Journal of Heat Transfer, 2021.
  * [74] K Kashinath, M Mustafa, A Albert, JL Wu, C Jiang, S Esmaeilzadeh, K Azizzadenesheli, R Wang, A Chattopadhyay, A Singh, et al.  Physics-informed machine learning: case studies for weather and climate modelling.  Philosophical Transactions of the Royal Society A, 379(2194):20200093, 2021.
  * [75] Mohammadamin Mahmoudabadbozchelou, Marco Caggioni, Setareh Shahsavari, William H Hartt, George Em Karniadakis, and Safa Jamali.  Data-driven physics-informed constitutive metamodeling of complex fluids: A multifidelity neural network (mfnn) framework.  Journal of Rheology, 65(2):179–198, 2021.
  * [76] Ehsan Kharazmi, Zhongqiang Zhang, and George Em Karniadakis.  hp-vpinns: Variational physics-informed neural networks with domain decomposition.  Computer Methods in Applied Mechanics and Engineering, 374:113547, 2021.
  * [77] Liu Yang, Xuhui Meng, and George Em Karniadakis.  B-pinns: Bayesian physics-informed neural networks for forward and inverse pde problems with noisy data.  Journal of Computational Physics, 425:109913, 2021.
  * [78] Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun.  Deep residual learning for image recognition.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 770–778, 2016.
  * [79] Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun.  Identity mappings in deep residual networks.  In European conference on computer vision, pages 630–645. Springer, 2016.
  * [80] E Weinan.  A proposal on machine learning via dynamical systems.  Communications in Mathematics and Statistics, 5(1):1–11, 2017.
  * [81] Lars Ruthotto and Eldad Haber.  Deep neural networks motivated by partial differential equations.  Journal of Mathematical Imaging and Vision, 62(3):352–364, 2020\.
  * [82] François Rousseau, Lucas Drumetz, and Ronan Fablet.  Residual networks as flows of diffeomorphisms.  Journal of Mathematical Imaging and Vision, 62(3):365–375, 2020\.
  * [83] Emilien Dupont, Arnaud Doucet, and Yee Whye Teh.  Augmented neural odes.  arXiv preprint arXiv:1904.01681, 2019.
  * [84] Hanshu Yan, Jiawei Du, Vincent YF Tan, and Jiashi Feng.  On robustness of neural ordinary differential equations.  arXiv preprint arXiv:1910.05513, 2019.
  * [85] Amir Gholami, Kurt Keutzer, and George Biros.  Anode: Unconditionally accurate memory-efficient gradients for neural odes, 2019.
  * [86] Tianjun Zhang, Zhewei Yao, Amir Gholami, Kurt Keutzer, Joseph Gonzalez, George Biros, and Michael Mahoney.  Anodev2: A coupled neural ode evolution framework, 2019.
  * [87] Aaron Tuor, Jan Drgona, and Draguna Vrabie.  Constrained neural ordinary differential equations with stability guarantees.  arXiv preprint arXiv:2004.10883, 2020.
  * [88] Julia Ling, Andrew Kurzawski, and Jeremy Templeton.  Reynolds averaged turbulence modelling using deep neural networks with embedded invariance.  Journal of Fluid Mechanics, 807:155–166, 2016.
  * [89] Andrea Beck, David Flad, and Claus-Dieter Munz.  Deep neural networks for data-driven les closure models.  Journal of Computational Physics, 398:108910, 2019.
  * [90] Laure Zanna and Thomas Bolton.  Data-driven equation discovery of ocean mesoscale closures.  Geophysical Research Letters, 47(17):e2020GL088376, 2020.
  * [91] Hugo Frezat, Guillaume Balarac, Julien Le Sommer, Ronan Fablet, and Redouane Lguensat.  Physical invariance in neural networks for subgrid-scale scalar flux modeling.  arXiv preprint arXiv:2010.04663, 2020.
  * [92] Alexis-Tzianni G Charalampopoulos and Themistoklis P Sapsis.  Machine-learning energy-preserving nonlocal closures for turbulent fluid flows and inertial tracers.  arXiv preprint arXiv:2102.07639, 2021.
  * [93] Klaus Neumann, Andre Lemme, and Jochen J Steil.  Neural learning of stable dynamical systems based on data-driven lyapunov candidates.  In 2013 IEEE/RSJ International Conference on Intelligent Robots and Systems, pages 1216–1222. IEEE, 2013.
  * [94] Mouhacine Benosman, Boris Kramer, Petros T Boufounos, and Piyush Grover.  Learning-based reduced order model stabilization for partial differential equations: Application to the coupled burgers’ equation.  In 2016 American Control Conference (ACC), pages 1673–1678. IEEE, 2016.
  * [95] Giorgos Mamakoukas, Orest Xherija, and TD Murphey.  Memory-efficient learning of stable linear dynamical systems for prediction and control.  arXiv preprint arXiv:2006.03937, 2020.
  * [96] Yu Wang, Qitong Gao, and Miroslav Pajic.  Deep learning for stable monotone dynamical systems.  arXiv preprint arXiv:2006.06417, 2020.
  * [97] Andrew J. Taylor, Victor D. Dorobantu, Hoang M. Le, Yisong Yue, and Aaron D. Ames.  Episodic learning with control lyapunov functions for uncertain robotic systems, 2019.
  * [98] N Benjamin Erichson, Michael Muehlebach, and Michael W Mahoney.  Physics-informed autoencoders for lyapunov-stable fluid flow prediction.  arXiv preprint arXiv:1905.10866, 2019.
  * [99] Thomas S. Parker and Leon O. Chua.  Stability of Limit Sets, pages 57–82.  Springer New York, New York, NY, 1989.
  * [100] Aleksandr Mikhailovich Lyapunov.  The general problem of the stability of motion.  International journal of control, 55(3):531–534, 1992.
  * [101] Michael Malisoff and Frédéric Mazenc.  Constructions of strict Lyapunov functions.  Springer Science & Business Media, 2009.
  * [102] James D Meiss.  Differential dynamical systems.  SIAM, 2007.
  * [103] Tim Sauer, James A. Yorke, and Martin Casdagli.  Embedology.  Journal of Statistical Physics, 65(3):579–616, Nov 1991.
  * [104] Philip Holmes, John L Lumley, Gahl Berkooz, and Clarence W Rowley.  Turbulence, coherent structures, dynamical systems and symmetry.  Cambridge university press, 2012.
  * [105] RE Langer.  Vito volterra, theory of functionals and of integral and integro-differential equations.  Bulletin of the American Mathematical Society, 38(9. P1):623–623, 1932.
  * [106] Chenjie Gu.  Qlmor: A projection-based nonlinear model order reduction approach using quadratic-linear representation of nonlinear systems.  IEEE Transactions on Computer-Aided Design of Integrated Circuits and Systems, 30(9):1307–1320, 2011.
  * [107] Elizabeth Qian, Boris Kramer, Benjamin Peherstorfer, and Karen Willcox.  Lift & learn: Physics-informed machine learning for large-scale nonlinear dynamical systems.  Physica D: Nonlinear Phenomena, 406:132401, 2020.
  * [108] Edward N. Lorenz.  Deterministic Nonperiodic Flow.  Journal of the Atmospheric Sciences, 20(2):130–141, March 1963\.
  * [109] A. C. Hindmarsh.  ODEPACK, a systematized collection of ODE solvers.  IMACS Transactions on Scientific Computation, 1:55–64, 1983.
  * [110] Henry D. I. Abarbanel.  Choosing Time Delays, pages 25–37.  Springer New York, New York, NY, 1996.
  * [111] Henry D. I. Abarbanel.  Choosing the Dimension of Reconstructed Phase Space, pages 39–67.  Springer New York, New York, NY, 1996.
  * [112] Julien Clinton Sprott.  Chaos and Time-Series Analysis.  Oxford University Press, Inc., New York, NY, USA, 2003.
  * [113] Edward N Lorenz.  Predictability: A problem partly solved.  In Proc. Seminar on predictability, volume 1, 1996.
  * [114] Alfio Quarteroni, Andrea Manzoni, and Federico Negri.  Reduced basis methods for partial differential equations: an introduction, volume 92.  Springer, 2015.
  * [115] Hugo Frezat, Guillaume Balarac, Julien Le Sommer, Ronan Fablet, and Redouane Lguensat.  Physical invariance in neural networks for subgrid-scale scalar flux modeling.  Physical Review Fluids, 6(2):024607, 2021.