<figure><img src="image_1.png" /></figure>

HAL
open science

# CD-ROM: Complementary Deep-Reduced Order Model

Emmanuel Menier, Michele Alessandro Bucci, Mouadh Yagoubi, Lionel

Mathelin, Marc Schoenauer

► To cite this version:

Emmanuel Menier, Michele Alessandro Bucci, Mouadh Yagoubi, Lionel Mathelin, Marc Schoenauer.

CD-ROM: Complementary Deep-Reduced Order Model. Computer Methods in Applied Mechanics and Engineering, 2023, 410, pp.115985. 10.1016/j.cma.2023.115985. hal-03846122

HAL Id: hal-03846122

https://hal.science/hal-03846122v1

Submitted on 31 Mar 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# CD-ROM: Complemented Deep - Reduced Order Model

Emmanuel Menier<sup>a,b,*</sup>, Michele Alessandro Bucci<sup>a</sup>, Mouadh Yagoubi<sup>b</sup>,
Lionel Mathelin<sup>c</sup>, Marc Schoenauer<sup>a</sup>

## Abstract

Model order reduction through the POD-Galerkin method can lead to dramatic gains in terms of computational efficiency in solving physical problems. However, the applicability of the method to non linear high-dimensional dynamical systems such as the Navier-Stokes equations has been shown to be limited, producing inaccurate and sometimes unstable models. This paper proposes a deep learning based closure modeling approach for classical POD-Galerkin reduced order models (ROM). The proposed approach is theoretically grounded, using neural networks to approximate well studied operators. In contrast with most previous works, the present CD-ROM approach is based on an interpretable continuous memory formulation, derived from simple hypotheses on the behavior of partially observed dynamical systems. The final corrected models can hence be simulated using most classical time stepping schemes. The capabilities of the CD-ROM approach are demonstrated on two classical examples from Computational Fluid Dynamics, as well as a parametric case, the Kuramoto-Sivashinsky equation.

**Keywords:** Reduced Order Models, Deep Learning, Neural ODE, Computational Fluid Dynamics

## 1. Introduction

### 1.1. General context

The simulation of complex physical processes requires solving high dimensional partial differential equations (PDE), often at large numerical cost. This constitutes a major limitation on the applicability of accurate simulation methods to engineering problems, such as model predictive control or iterative design optimisation. Despite this complexity, physical problems can often be suitably

<footnote><sup>a</sup>TAU, Inria / LISN, Université Paris-Saclay / CNRS, 1 Rue Raimond Castaing, 91190, Gif-sur-Yvette, France</footnote>

<footnote><sup>b</sup>IRT - SystemX, 8 Avenue de la Vauve, 91120 Palaiseau, France</footnote>

<footnote><sup>c</sup>LISN, Université Paris-Saclay & CNRS, Campus Universitaire bât.507, Rue du Belvédère, 91405, Orsay, France</footnote>

<footnote>*Corresponding author
Email address: emmanuel.menier@inria.fr (Emmanuel Menier)</footnote>
approximated by a reduced number of dominant structures, suggesting the ex-
istence of low-dimensional manifolds supporting the underlying dynamics and
fueling the development of reduction methods.

The construction of representative reduced order models often relies on the identification of a suitable modal basis to project a description of the system onto. These modes are most often extracted from data, through different methods. For example, the well-known Proper Orthogonal Decomposition (POD)-Galerkin method ([1]) has been applied to the reduction of linear systems with great results ([2]). It was also applied to more complex problems such as the Navier Stokes equations with varying degrees of success ([3, 4]). A three dimensional Galerkin model was shown in [5] to capture the transient and oscillation regimes of the cylinder flow, provided the POD basis included a mode tailored to support the transition. Models based on Galerkin POD formulations have the advantage of reducing computational cost while preserving a high level of interpretability in terms of coherent structures and involved dynamics. Nevertheless, despite these promising results, the predictive capabilities of Galerkin models were shown to be limited and often lead to inaccurate reduced models. Indeed, while the POD is optimal for the reconstruction error in a suitable norm, the Galerkin projection might not be ([2]).

Because of these limitations, some research efforts have been focused on
deriving purely data-driven models via sparse regression methods ([6, 7]), cluster
modeling ([8]) or Koopman theory ([9]) to cite just a few. The representation
of the dynamics with the Koopman operator has garnered a lot of interest in
recent years, both through the DMD algorithm ([10]) and its extensions (e.g.,
EDMD and KDMD [11, 12]) which leverage nonlinear state transformations to
simplify the modeling problem.

These purely data-driven works are part of a wider trend towards hybridisa-
tion of simulation methods with machine learning techniques. Specifically, neu-
ral networks which have been applied to numerous physical modeling problems,
from black box resolutions through graphs networks allowing the resolution of
problems on unstructured meshes ([13]), to more interpretable approaches like
the physics informed neural networks proposed in [14]. Over the last few years,
several works have proposed unsupervised learning techniques to avoid the mesh
construction (mesh-free methods). These algorithms train deep neural networks
to approximate PDE solutions by satisfying the differential operator, initial con-
ditions and boundary conditions for a specific PDE. [15] used fully connected
layers to approximate the solution on complex geometries. [16, 17] discussed the
possibility to solve high dimensional problems through neural networks. The
authors of [17] proved that the neural network converges to the PDE as the num-
ber of hidden units increases. A more detailed review on machine learning/CFD
hybridization can be found in [18]. At the same time, the deep learning com-
munity has developed powerful tools by adapting physical modeling concepts to
deep learning, for instance the Neural ODE approach was proposed for the mod-
eling of continuous transformations ([19]), stable neural architectures were also
developed by embedding neural networks with invariant structures ([20, 21, 22]).
Neural networks have also been applied in the context of model order reduction,
notably to learn useful representations of a system’s state ([23, 24, 25, 26]).

As part of this interest in hybridization, we propose a novel architecture for the development of closures for Galerkin reduced models. Despite their limited performance with nonlinear systems, the physical information embedded in their structure can be retained and used as a basis for the construction of a hybrid reduction approach. Applications of this idea have already been proposed in the literature ([27]) and are discussed in the following paragraph. We refer the reader to [28] for a thorough review on the topic.

## 1.2. Related Work

Reduced order modeling methods offer powerful and straightforward industrial applications, thus, their improvement has been the subject of a large body of literature. The proposed solutions can be divided into *intrusive* and *non-intrusive* methods. *Intrusive* approaches aim at learning a closure model for the available ROM such that high fidelity data are fitted, while *non intrusive* approaches propose to completely replace pre-existing models with learnt black-box forecasting methods. Several works have proposed to leverage deep learning methods to develop non intrusive forecasting models on low dimensional spaces and represent physical simulation problems ([29, 30, 31, 32, 33]). Specifically, the necessity of exploiting temporal information to reconstruct accurate dynamics has been underlined in [34, 35, 36, 37]. More theoretically grounded works have also proposed to develop intrusive non-Markovian closure models for existing Reduced Order Modeling approaches ([27, 38]). They are motivated by the Mori-Zwanzig formalism ([39]), which provides a theoretical framework for the modeling of partially observed systems. Similar to these efforts, our work aims at developing a theoretically grounded reduced modeling method for the forecasting of physical systems, preserving the physical insights provided by the projected PDE describing the problem. Contrary to most previous works, which are based on the use of recurrent neural networks, our proposal is constructed around a time-continuous memory formulation with numerous advantages over discrete time models. Moreover, we underline below the higher degree of interpretability of our solution. The main novelties of this work are listed below and Table 1 also provides a summary of the potential crossovers between existing methods and the present work.

*   **Intrusivity:** Closure modeling has been a topic of interest since the early days of numerical simulation, thus, developing *intrusive* correction models which combine with existing physical models is not a novel approach in and of itself. However, most deep learning approaches do not take this route and propose to learn forecasting models from scratch, ignoring the underlying physical laws. In this work, we show that the Galerkin ROM can be used inside the training loop to optimise the closure model in an *a-posteriori* fashion, that is, by simulating the whole model and assessing its performance. Embedding the existing ROM in the training strategy allows us to leverage pre-existing physical information rather than replace it with a physics-agnostic model.
<table><caption>Table 1: Comparison of closure approaches from the recent literature.</caption><thead><tr><th></th><th>Non-Markovian</th><th>Non Intrusive</th><th>Continuous</th><th>Memory interpretability</th></tr></thead><tbody><tr><td>Wang et al. ([27])</td><td>✓</td><td></td><td></td><td>✓</td></tr><tr><td>Pawar et al. ([38])</td><td>✓</td><td></td><td></td><td></td></tr><tr><td>Maulik et al. ([30])</td><td>✓</td><td>✓</td><td>✓</td><td></td></tr><tr><td>Vlachas et al. ([40])</td><td>✓</td><td>✓</td><td></td><td>✓</td></tr><tr><td>Wu et al. ([37])</td><td>✓</td><td>✓</td><td></td><td>✓</td></tr><tr><td>Maulik et al. ([35])</td><td>✓</td><td>✓</td><td></td><td></td></tr><tr><td>Pawar et al. ([31])</td><td>✓</td><td>✓</td><td></td><td></td></tr><tr><td>Our work</td><td>✓</td><td></td><td>✓</td><td>✓</td></tr></tbody></table>

* **Continuity:** The proposed model is embedded with a novel time continuous memory formulation. This increases the applicability of the model as it can be plugged in any initial value problem solver without being biased against specific time-step and/or numerical scheme choices made during training. Moreover, this flexibility implies that this work can be used to model arbitrarily stiff problems through the use of adaptive time-stepping schemes. Finally, the continuous structure allows the model to be used in combination with irregularly spaced data, often encountered in real-life problems, with no additional considerations to interpolation between samples.

* **Memory Interpretability:** Contrary to classical recurrent neural networks such as the Long Short Term Memory (LSTM) used in the state of the art, the proposed memory formulation was specifically designed for numerical simulation purposes. Particularly, the time evolution of the memory has a closed form solution, which means that the memory term can be initialised to any desired degree of precision. Our proposed formulation also allows for the evaluation of the *time persistence* of information in memory, increasing the overall interpretability of the model.

* **End to End Training:** In contrast with other works which propose to learn the closure model in an *a-priori* fashion, i.e. by learning the dynamics correction as a standalone regression problem, we integrate the imperfect model directly within the training strategy. Thus the correction model accurately learns to compensate for the sensitivity of the Galerkin ROM and account for the long term effects of unstable, low energy modes in the system. Indeed, it has been observed in the literature that this strategy lead to more stable and accurate models ([41, 24, 42]).

The outline of the paper is as follows: the POD-Galerkin reduction method is introduced in Section 2 and its limitations are discussed. Section 3 details our main contribution on the derivation of a continuous in time correction architecture for reduced order models. Motivation for our work through comparisons
with existing approaches is provided in Section 4, while the selected test cases and results are respectively discussed in Section 5 and 6. Section 7 concludes the paper.

## 2. Model reduction approach

In this Section we aim at describing the POD-Galerkin model reduction strategy, here illustrated in the context of fluid dynamics. In particular we focus on the incompressible Navier-Stokes equations owing to their relevance in engineering applications and since they represent a challenging test-bed for model reduction techniques. The dynamics of the velocity field $\mathbf{u}(x, t)$ and pressure field $p(x, t)$ are governed by:

$$ \frac{\partial \mathbf{u}}{\partial t} + (\mathbf{u} \cdot \nabla)\mathbf{u} = -\nabla p + \frac{1}{Re} \nabla^2 \mathbf{u} \quad (1) \\ \nabla \cdot \mathbf{u} = 0 $$

with $Re$ the Reynolds number. Computationally demanding operations are required to numerically solve Eq. (1). The typical approach of POD Galerkin Reduced Order Models (ROMs) is then employed to identify a low-dimensional approximation of the Navier-Stokes equations based on time invariant spatial features. The method, and its derivation for the Navier-Stokes equations, is presented in Section 2.1, and its limitations are discussed in Section 2.2. For a deeper introduction to POD Galerkin reduced order models we refer the interested reader to [43, 44].

### 2.1. POD Galerkin models

After gathering data from experiments or the numerical simulation of Eq. (1), one can use the Proper Orthogonal Decomposition, also known as Principal Component Analysis (PCA), to identify dominant structures in the data. The POD method relies on a matrix, referred to as the snapshot matrix, $\mathbf{U} = \{\mathbf{u}(t_i)|i = 1, \dots, n_t\}$ whose columns $\mathbf{u}(t) \in \mathbb{R}^{d_{n_x}}$ correspond to the $d$-dimensional state of the system (for example, the three components of the fluid velocity) indexed upon the $n_x$ points of a spatial grid at time $t$. In the sequel, it is supposed that $n_x \ge n_t$, which is common in numerical simulation applications, but not a requirement of the POD method.

Let us consider the thin singular value decomposition (SVD) of this real-valued matrix: $\mathbf{U} = \mathbf{V}\mathbf{\Sigma}\mathbf{W}^T$, where the matrix $\mathbf{V} \in \mathbb{R}^{d_{n_x} \times n_t}$ and $\mathbf{W} \in \mathbb{R}^{n_t \times n_x}$ respectively hold the left and right singular vectors of $\mathbf{U}$. The diagonal matrix $\mathbf{\Sigma} \in \mathbb{R}^{n_t \times n_t}$ holds the singular values of $\mathbf{U}$, usually arranged such that $\sigma_1 \ge \sigma_2 \ge \dots \ge \sigma_{n_t} \ge 0$.

Columns of $\mathbf{V}$ are the time-invariant spatial modes defining a POD basis. These POD modes are orthonormal $\langle \mathbf{v}_i, \mathbf{v}_j \rangle = \delta_{ij}$, with $\delta_{ij}$ the Kronecker delta and $\langle \cdot, \cdot \rangle$ is here Euclidean: $\langle \mathbf{v}_i, \mathbf{v}_j \rangle = \mathbf{v}_i^T \mathbf{v}_j$. These modes are useful for dimensionality reduction because, for any $r < n_t$, the subspace spanned by the basis $\mathbf{V}_r = \{\mathbf{v}_1, \mathbf{v}_2, \dots, \mathbf{v}_r\}$ optimally approximates the data $\mathbf{U}$ over the set
of $dn_x \times r$ matrices in the sense that it minimizes the reconstruction error $E_r$ defined as:

$$E_r = \|\mathbf{U} - \mathbf{V}_r \mathbf{V}_r^\mathrm{T} \mathbf{U}\|_F, \quad \forall r \in [1, \dots, n_t]. \qquad (2)$$

This error can be related to the sum of the discarded singular values: $E_r^2 = \sum_{k=r+1}^{n_t} \sigma_k^2$. This implies that the information captured by the first $r$ modes in the POD basis can be quantified by looking at the following ratio:

$$R(r) = \frac{\sum_{k=1}^{r} \sigma_k^2}{\sum_{k=1}^{n_t} \sigma_k^2}. \qquad (3)$$

Introducing a projector $\mathbf{P} = \mathbf{V}_r \mathbf{V}_r^\mathrm{T}$ of $\mathbb{R}^{dn_x}$ onto the invariant subspace spanned by the basis $\mathbf{V}_r$, and $\mathbf{Q} = (\mathbf{I} - \mathbf{P})$ its orthogonal complement, the full state vector $\mathbf{u}(t)$ can be decomposed as:

$$\mathbf{u}(t) = \mathbf{P}\mathbf{u}(t) + \mathbf{Q}\mathbf{u}(t). \qquad (4)$$

The key idea is to select a suitable $r$ in Eq. (3) based on an energy criterion such that $\mathbf{Q}\mathbf{u}(t)$ can be neglected for the problem under consideration. The full state $\mathbf{u}(t)$ can then be approximated as:

$$\mathbf{u}(t) \approx \tilde{\mathbf{u}}(t) = \mathbf{V}_r \mathbf{a}(t) \qquad (5)$$

where $\mathbf{a}(t) := \mathbf{V}_r^\mathrm{T} \mathbf{u}(t)$ is the vector containing temporal coefficients associated with the spatial modes in $\mathbf{V}_r$. From the approximation defined in Eq. (5), model reduction arises whenever $r \ll n_x$. Under this modal decomposition ansatz, the original spatio-temporal problem of solving the PDE is reduced to the computation of the temporal coefficients $\mathbf{a}(t)$ as the spatial information is expressed in the invariant POD basis $\mathbf{V}_r$. However, to achieve a gain in computational cost, one needs a way to compute the reduced coordinates vector without accessing the full order solution. This is where the Galerkin projection method can be used to obtain an equation for the reduced coordinates. The method is based on the reduction of the original model equation:

$$\frac{d\mathbf{u}(t)}{dt} = \mathbf{f}(\mathbf{u}(t)) \qquad (6)$$

where $\mathbf{f}$ corresponds to the full order dynamics discretised in space. By projecting this equation onto the time-independent POD basis, a dynamical model for the reduced coordinates is obtained:

$$\frac{d\mathbf{a}(t)}{dt} = \mathbf{V}_r^\mathrm{T} \mathbf{f}(\mathbf{u}(t)). \qquad (7)$$

With this Galerkin projection, the dimension of the problem is reduced but the model still requires knowledge of the full state $\mathbf{u}(t)$ to be resolved. To obtain a model that solely depends on the coordinates $\mathbf{a}(t)$, the dynamics $\mathbf{f}$ of
the model have to be evaluated from the approximate form of the state (see Eq. (5)). This leads to the following system of r ODEs:

$$ \frac{da(t)}{dt} \approx V_r^T f(V_r a(t)). \qquad (8) $$

If $f$ is the right hand side of the incompressible Navier-Stokes in Eq. (1) and the POD basis is retrieved from the velocity field snapshot matrix, the approximation in (8) leads to the following ROM:

$$ \frac{da_i(t)}{dt} = \sum_{j=1}^{r} \tilde{\mathcal{L}}_{i,j} a_i(t) + \sum_{j=1}^{r} \sum_{k=1}^{r} \tilde{Q}_{i,j,k} a_j(t) a_k(t), \quad \forall i = 1, \dots, r \quad (9) $$

where tensors $\tilde{\mathcal{L}}$ and $\tilde{Q}$ are reduced versions of the dissipative and advective parts of the original Navier Stokes equations. Because operators $\mathcal{L} = \Delta \mathbf{u}(x,t)$ and $\tilde{Q} = \nabla(\mathbf{u}(x,t)\mathbf{u}(x,t))$ are spatial operators (here $\mathbf{u}(x,t) \in \mathbb{R}^d$ refers to the velocity vector), they can be pre-assembled from the POD modes. In contrast with this favorable situation, other PDE problems might necessitate that the dynamics be evaluated from the reconstructed solution $\mathbf{u}(t) \approx V_r \mathbf{a}(t)$, which can imply a significant cost and require additional reduction work, [45]. It is worth noting that the pressure term is not present in Eq. (9). This is due to the fact that, in the full order model, the pressure term accounts for the incompressibility condition $(\nabla \cdot \mathbf{u} = 0)$ and can be interpreted as a Lagrange multiplier used to ensure this constraint during the simulation.

Because the reconstructed solution is a linear combination of the POD modes which are all divergence-free, the incompressibility condition is always satisfied and the pressure term can then be ignored without significant accuracy loss ([2, 5]).

## 2.2. Limitations of POD-Galerkin models

Because POD-Galerkin models rely on the projection of the original model equations, they retain part of the structure from the full order physical model. This makes them relatively robust to new initial conditions comparatively to purely data-based approaches which can struggle to model conditions unseen during training, as well as suffer from overfitting depending on the way they were trained. It should be noted however that approximation errors are embedded by construction in their structure. Equation (7) shows that an exact expression for the dynamics of the reduced coordinates is obtained by projection of the model equation on the POD basis. If the orthogonal complement to $V_r$ is not neglected, equation (7) reads:

$$ \frac{da(t)}{dt} = V_r^T f(\mathbf{u}(t)) = V_r^T f(\mathbf{P}\mathbf{u}(t) + \mathbf{Q}\mathbf{u}(t)), \qquad (10) $$

with $\|\mathbf{Q}\mathbf{u}(t)\|_F \ll \|\mathbf{P}\mathbf{u}(t)\|_F$ if $r$ is large enough. Considering a first order approximation of Eq. (10), the temporal dynamics of the POD coefficients can be evaluated from:
$$ \frac{da(t)}{dt} \approx V_r^T f(V_r a(t)) + \underbrace{V_r^T J_{Pu(t)} Q u(t)}_{\mathcal{R}} \quad (11) $$

with $J_{Pu(t)}$ being the Jacobian of the discretised nonlinear system $f$, evaluated at the state $Pu(t)$. The residual term $\mathcal{R}$ accounts for the dynamics of the complement of the POD basis in the space spanned by the POD basis. Classical POD-Galerkin reduced models ignore this error term and directly compute the dynamics of the model from the approximate solution as presented in Section 2.1.

This approximation means that small errors on the dynamics will compound over time and lead to significant discrepancies between the true and simulated trajectories. This is especially true in the case of nonlinear dynamical systems where orthogonal projection on the POD basis can suppress an important part of the dynamics. In the case of the Navier Stokes equations, POD-Galerkin models have been shown to fail to reproduce the dynamics even in simple cases like the flow over a cylinder [5]. It should be noted that different projection schemes have been proposed to alleviate this issue ([46]), however, these often require additional considerations on the construction of the reduced model.

The aim of this paper is to retain the simplicity of the POD-Galerkin method, and learn the complement of the ROM using deep learning methods, thus, we call our proposed method **Complemented Deep - Reduced Order Model**.

## 2.3. Non-Markovianity and Takens' theorem

The residual depends on information from a subspace orthogonal to the span of the POD basis. This means that an accurate correction model cannot be directly computed from the reduced state $a(t)$. However, we leverage the fact that the information lost by projection of the full order state can be retrieved by considering past states of the system. This hypothesis is formalized by the Takens' theorem ([47]), which states that, under mild conditions, the dynamics of a state vector can be reproduced by constructing a time-embedding from time-lagged observables: $z(t) = (z(t), z(t-\tau), z(t-2\tau), \dots, z(t-k\tau))$, with $k$ large enough.<fnref n="1" /> This is illustrated in Figure 1, which shows the Lorenz attractor observed via its embedded $X$-component. By constructing a 3-dimensional embedding of the obtained time-series, an attractor is obtained, which preserves the topology of the true attractor (e.g., symmetries, correlation dimension, etc.).

This suggests that the correction term $\mathcal{R}$ in POD-Galerkin reduced order models is non-Markovian and should consider past states of the system ($a(t)$ -

<footnote n="1">¹We tacitly assumed here that $z$ is a suitable observable. Observability analysis goes beyond the scope of this paper. Nevertheless, it is well known that higher harmonic POD modes are enslaved to dominant POD modes [7, 48]. Therefore, in the rest of the paper, we will assume that the dynamics of the unobserved space spanned by the columns of $Q$ can be retrieved from past observations of the dynamics in the space defined by $P$.</footnote>
<figure><img src="image_3.png" /><figcaption>Figure 1: Illustration of the Takens’ delay embedding theorem. Left: Original Lorenz attractor; middle: Delayed time series of the X coordinate; right: Reconstructed attractor.</figcaption></figure>

τ), a(t − 2 τ ) , . . . ). However, such a discrete time-embedding is not well suited in the context of continuous time models such as equation (9). Indeed, using discrete time steps in combination with adaptive step time marching schemes, such as the Runge-Kutta method [49], would require specific considerations about interpolation between the simulation steps and the required embedding steps. To address this time-continuity issue while retaining a non-Markovian correction structure, we propose to use delay differential equations (DDE) with a continuous embedding of the past information:

$$
\frac{d}{dt} \mathbf{a}(t) = \mathbf{g}(\mathbf{a}(t), \mathbf{y}(t)) , \quad \mathbf{y}(t) = \int_{-\infty}^{t} e^{(\tau-t)\lambda} \mathbf{a}(\tau) d\tau \quad (12)
$$

with λ ∈ R + sufficiently large for the integral above to be bounded.

These equations retain information from past states of the system in a time-
continuous manner and are used for a number of modeling applications such
as epidemiology or population dynamics ([50, 51]). In this formulation, the
dynamics g depend both on the partially observed states of the system a(t)
and a memory variable y(t), corresponding to the integral of past observables,
damped in time by an exponential decay. In fact, the memory term y can
be defined in many ways, depending on the problem at hand. However, the
exponential decay formulation was chosen because i) it provides the model with
the ability to consider recent states of the system while older observations are
discarded and ii) it can be solved by directly augmenting the original system
with a second ODE for the memory:

$$
\begin{align}
\frac{d}{dt} \mathbf{a}(t) &= \mathbf{g}(\mathbf{a}(t), \mathbf{y}(t)), \notag \\
\frac{d}{dt} \mathbf{y}(t) &= \mathbf{a}(t) - \lambda \mathbf{y}(t). \tag{13}
\end{align}
$$

The exponential decay e ( τ − t ) λ acts as a filter of width 1 /λ on the observables evolution and we show in a later section that the value of the decay rate ( λ )
can be learned in a data driven setting. Other applications of the same augmented ODE exist in literature but with different purposes such as modelling of subgrid-scales in Large-Eddy Simulations (LES) [52] or to find unstable-steady solutions of the Navier-Stokes equations [53]. As mentioned earlier, we have made the convenient choice of the exponential kernel since it can be described by a linear ODE and easily constrained to model dissipative dynamics, but alternative kernels exist in the DDE literature which could be considered if the exponential kernel became too constraining for certain cases.

It is worth noting that the memory variables $\mathbf{y}$ have the same dimension as the observations $\mathbf{a}$ in (13). This limitation might introduce a significant information bottleneck in the model. Indeed, Takens' theorem states that the dimension required to obtain a satisfactory embedding can go as high as twice the intrinsic dimension of the true attractor. Although there is no similar result for the continuous case, the limited dimension of the memory $\mathbf{y}$ may prevent deriving an accurate correction model. As a result, we define an encoding map $\mathbf{E}: \mathbb{R}^r \to \mathbb{R}^{n_E}$, used to lift the observations $\mathbf{a}$ to a higher dimensional space to increase the dimension of the memory:

$$ \mathbf{y}(t) = \int_{-\infty}^{t} e^{(\tau-t)\lambda} \mathbf{E}(\mathbf{a}(\tau)) d\tau. \quad (14) $$

In fact, the use of such an encoding map to unfold non-linear dynamics and recover a linear ODE is rooted in Koopman theory as each encoded coordinate can be considered as an observable of the original state. While approaches such as dictionary learning ([11]) and kernel methods ([12]) have been proposed to learn these observables, we use neural networks to avoid additional optimization considerations and retain flexibility, similar to the works of [23, 25]. Using the modified DDE architecture (Eq.(14)) to close reduced order models, the correction operator $\mathcal{R}$ acting on the memory $\mathbf{y}(t)$ becomes an application from memory space to phase space: $\mathbb{R}^{n_E} \to \mathbb{R}^r$. Finally, the following augmented reduced order model is obtained:

$$ \begin{aligned} \frac{d}{dt} \mathbf{a}(t) &= \mathbf{V}_r^\top \mathbf{f}(\mathbf{V}_r \mathbf{a}) + \mathcal{R}(\mathbf{y}), \\ \frac{d}{dt} \mathbf{y}(t) &= \mathbf{E}(\mathbf{a}) - \lambda \mathbf{y}. \end{aligned} \quad (15) $$

This proposed augmented ROM architecture has a similar form to the Mori-Zwanzig formalism [39] which derives an equation for the dynamics of partially observed systems. Although the parallel with the CD-ROM approach is not formally established, we discuss these similarities in section 4.1 as it provides additional motivation for the above choices. Finally, the augmented ROM formulation is summarised in Figure 2 to help illustrate the idea.

### 3. Data driven learning of the residual

Although the motivation for the existence of the residual $\mathcal{R}$ and encoding map $\mathbf{E}$ was outlined above, little information about their form can be derived
<figure><img src="image_4.png" /><figcaption>Figure 2: Visualisation of the CD-ROM approach. The full order solution snapshots are projected on the POD basis to obtain time series of the reduced coordinates (a). These reduced coordinates are then augmented with a memory y(t). Finally, the augmented dynamical model is presented on the right hand side of the image.</figcaption></figure>

from the previous expressions. In this context, we leverage the universal function approximator property of artificial neural networks [54] to model the missing terms.

## 3.1. Neural Networks

Our work relies on feed-forward neural networks with a very simple structure. These networks can be used to learn any smooth nonlinear continuous application $\Psi: \mathbb{R}^r \to \mathbb{R}^k$ by optimising the weights of a sequence of L layers:

$$ \Psi(\mathbf{a}) = \eta_L \circ \eta_{L-1} \circ \cdots \circ \eta_1(\mathbf{a}) \qquad (16) $$

$$ \eta_l = \sigma_l(W_l \eta_{l-1} + b_l), \quad \forall l \in \{2, \dots, L\} \qquad (17) $$

where $\sigma_l$ is a nonlinear activation function, and the dimension of $\eta_l$ corresponds to the number of neurons in the layer $l$. It has been shown that, provided the dimension of the layer is high enough, the trainable parameters $W_l$ and $b_l$ can be optimised to approximate any function, [54]. The encoder $E(\mathbf{a})$ and the residual $\mathcal{R}(\mathbf{y})$ are both approximated with neural networks with parameters $\theta_E = (\{W_l, b_l\}, l = 1, \dots, L_E)$ and $\theta_R$ respectively.

## 3.2. Memory time scales

Physical systems often involve a variety of phenomena each evolving at different time scales. Capturing these phenomena can be critical to accurately model these systems, which is why the memory should be able to retain information at different rates. The time scales accounted for by the memory are driven by the parameter $\lambda$ of the exponential decay in Eq. (12), which acts as a low-pass filter on the encoded trajectory (see Figure 3 for an illustration).

To retain information at different rates, $\lambda$ can be adjusted for each observable $E_i(\mathbf{a})$. Equation (15) is then modified accordingly with the single $\lambda$ parameter replaced with a diagonal matrix $\Lambda \in \mathbb{R}_+^{n_E}$ whose entries can be learned to select time scales relevant to the system at hand. This finally results in the CD-ROM architecture:
<figure><img src="image_5.png" /><figcaption>Figure 3: Superposition of sinusoidal signals of frequencies 1, 10 and 100 Hz, filtered by an exponential decay with a decay rate $\lambda = 50$. The higher frequency (100 Hz) is filtered out while the other two are reproduced in the memory signal. This implies that recent events in memory due to the 100 Hz frequency are filtered out so that the memory is mostly driven by older events associated with lower frequencies. Note that the memory signal has been scaled by a factor $\lambda$ for clarity.</figcaption></figure>

$$ \begin{aligned} \frac{d}{dt} \mathbf{a}(t) &= \mathbf{V}_r^\top \mathbf{f}(\mathbf{V}_r \mathbf{a}) + \mathcal{R}(\mathbf{y}; \theta_R), \\ \frac{d}{dt} \mathbf{y}(t) &= \mathbf{E}(\mathbf{a}; \theta_E) - \mathbf{\Lambda} \mathbf{y}. \end{aligned} \qquad (18) $$

### 3.3. Training strategy

Optimising the parameters of the model (18) requires consideration of past states of the system as well as their impact on the dynamics in the future. In some sense, the problem is similar to the optimisation of classical recurrent neural networks for the simulation of dynamical systems [40]. The major difference is that the aim of the present work is to derive a continuous time dynamical model, while recurrent networks have traditionally been used to model transitions between discrete time instants. As a result, we cannot use standard backpropagation through time to optimise the model. Instead, the NeuralODE approach [19] is used. Initially proposed as the continuous equivalent of residual networks ([55]), Neural ODEs are based on the well established adjoint backpropagation method ([56]). It is briefly presented here while a more detailed derivation is given in Appendix A. The adjoint backpropagation method can be used to solve constrained optimization problems of the form:

$$ \min_{\theta} J(\mathbf{x}(T; \theta)) \qquad (19) $$

$$ \text{s.t.} \quad \frac{d}{dt} \mathbf{x}(t) = \mathbf{f}(\mathbf{x}; \theta) \qquad (20) $$

$$ \mathbf{x}(0) = \mathbf{x}_0. \qquad (21) $$

The method allows for the computation of the gradient $\frac{d}{d\theta} J$, at roughly twice the cost of integrating forward in time from $t = 0$ to $t = T$, through the following expression:
$$ \frac{d}{d\theta} J = \int_{0}^{T} -\mu^{\mathrm{T}} \frac{\partial f}{\partial \theta} dt \qquad (22) $$

where the adjoint state $\mu$ can be computed by integrating the adjoint ODE:

$$ \frac{d}{dt} \mu^{\mathrm{T}}(t) = -\mu^{\mathrm{T}} \frac{\partial f}{\partial \mathbf{x}}|_{t}, \qquad (23) $$

$$ \mu^{\mathrm{T}}(T) = \frac{\partial J}{\partial \mathbf{x}}|_{T}. \qquad (24) $$

Computing the gradient of $J$ with respect to the parameters thus requires integration of the adjoint ODE backwards in time, from $t = T$ to $t = 0$. If the dynamics $f$ are expressed by a Neural Network, the required vector-Jacobian products $\mu^{\mathrm{T}} \frac{\partial f}{\partial \mathbf{x}}|_{t}$ and $\mu^{\mathrm{T}} \frac{\partial f}{\partial \theta}|_{t}$ can seamlessly be evaluated through automatic differentiation by any classical deep learning framework.

As explained in Section 2.1, reduced models of the incompressible Navier-Stokes equations can be expressed directly in terms of the reduced coordinates $\mathbf{a}$ and evaluated through simple tensorial expressions (see equation (9)). This means that the ROM dynamics can be directly computed and back-propagated through once the reduced operators $\tilde{\mathcal{L}}$ and $\tilde{\mathcal{Q}}$ are assembled (see equation 9). The Deep ROM architecture (18) can then be optimised within the Neural ODE framework by concatenating the reduced state and the memory into a single state vector $\mathbf{z}(t) = [\mathbf{a}(t), \mathbf{y}(t)]$, with dynamics $f(\mathbf{z}; \theta_{\mathbf{E}}, \theta_{\mathbf{R}}, \mathbf{\Lambda})$:

$$ f(\mathbf{z}; \theta_{\mathbf{E}}, \theta_{\mathbf{R}}, \mathbf{\Lambda}) = \begin{pmatrix} \tilde{\mathcal{L}}\mathbf{a} + \mathbf{a}^{\mathrm{T}}\tilde{\mathcal{Q}}\mathbf{a} & + & \mathcal{R}(\mathbf{y}; \theta_{\mathbf{R}}) \\ \mathbf{E}(\mathbf{a}; \theta_{\mathbf{E}}) & - & \mathbf{\Lambda}\mathbf{y} \end{pmatrix}. \qquad (25) $$

Since the last modes in the POD basis $\mathbf{V}$ are typically associated to low-energy dissipative scales [28], if $r$ is not large enough, the ROM model retrieved by $\mathbf{V}_r$ might be unstable and diverge after a few time integration steps. This potentially unstable primal model is however embedded in the optimization framework so that the residual is trained accordingly, resulting in a stable, accurate, model.

### 3.4. Training Data and Memory Initialisation

To optimise the parameters of equation (25) through the Neural ODE approach, one only needs knowledge of the true trajectory of the reduced coordinates $\mathbf{a}^*(t)$. As presented in Section 2.1, this trajectory data can be obtained by projecting the solutions computed with the full order solver on the POD basis:

$$ \mathbf{a}^*(t) = \mathbf{V}_r^{\mathrm{T}} \mathbf{u}(t). \qquad (26) $$

If the snapshots are sampled with a time interval $\Delta_t$, simulating the corrected ROM (25) for $n_t$ time steps leads to the following $L^2$ mean squared error:
$$J = \frac{1}{n_t} \sum_{i=1}^{n_t} \| \mathbf{a}(i \Delta_t) - \mathbf{a}^*(i \Delta_t) \|_2^2 \quad (27)$$

for which a gradient can be computed by integrating the adjoint equation (23). It should be noted that the memory term needs to be initialised by evaluating the memory integral at time $t = 0$ which requires knowledge of the past of the true trajectory:

$$\mathbf{y}(0) = \int_{-\infty}^{0} \mathbf{E}(\mathbf{a}^*(\tau); \theta_{\mathbf{E}})e^{\mathbf{\Lambda}\tau}d\tau. \quad (28)$$

To be able to compute this integral, the infinite horizon of the memory integral can be relaxed by defining a finite $\tau_{\min}$ from the longest time scale $\lambda_{\min}$ of the matrix $\mathbf{\Lambda}$ and a threshold $\epsilon \ll 1$, chosen to be small enough such that the relative error made on the initial memory is sufficiently small:

$$\tau_{\min} = -\frac{\log \epsilon}{\lambda_{\min}}, \quad (29)$$

$$\mathbf{y}(0) \approx \int_{-\tau_{\min}}^{0} \mathbf{E}(\mathbf{a}^*(\tau); \theta_{\mathbf{E}})e^{\mathbf{\Lambda}\tau}d\tau. \quad (30)$$

Because this initial memory term directly depends on the parameters of the encoder and the matrix $\mathbf{\Lambda}$, it needs to be re-computed at each training epoch. This can be done very efficiently on a GPU through a simple trapezoidal approximation of the integral in Eq. (30). It should be noted that, while the initialization of the memory is necessary to obtain an exact model of the system at hand, it could limit the applicability of the method to certain real life settings. However, there are ways to make it less critical, such as initializing the memory with white noise during training. This would of course impact the accuracy of the model, depending on the system at hand.

To conclude this section, algorithm 1 summarises the NeuralODE training procedure.

## 4. Interpretation of the model

Before presenting the selected test cases and results, further justification and insights into the model are discussed in this Section. First, the Mori-Zwanzig formalism is introduced to frame our model in the context of dynamical systems theory. Then, we study how it can be compared to purely data driven approaches, such as reservoir computing.

### 4.1. The Mori-Zwanzig formalism

The Mori-Zwanzig formalism ([57, 58, 39]) provides a closed form for the dynamics of partially observed systems by distinguishing three separate terms:
<caption>**Algorithm 1** Training the CD-ROM as a Neural ODE</caption>

**Require:** $f(\mathbf{z}; \theta)$ the dynamics of the CD-ROM, $J$ the loss, $\eta$ the learning rate
**for** $i \leftarrow 1$ to $N$ **do**
&emsp;Sample a batch of trajectories $\mathbf{a}_{-\tau_{\min}}^{*} \rightarrow T$
&emsp;Compute the initial memory $\mathbf{y}_0$
&emsp;$\mathbf{z}_0 \leftarrow [\mathbf{a}_0^*, \mathbf{y}_0]$
&emsp;$[\mathbf{a}_t, \mathbf{y}_t] \leftarrow \mathbf{z}_0 + \int_0^t f(\mathbf{z}; \theta)dt$
&emsp;Compute the loss $J(\mathbf{a}_0 \rightarrow T, \mathbf{a}_0^*)$
&emsp;$\mu_T \leftarrow \frac{dJ}{dz_T}$
&emsp;$\mu_t \leftarrow \mu_T - \int_T^t \mu^T \frac{\partial f}{\partial z_t} dt$
&emsp;$\frac{d}{d\theta} J \leftarrow \int_0^T -\mu_t^T \frac{\partial f}{\partial \theta} dt$
&emsp;$\theta \leftarrow \theta - \eta \frac{d}{d\theta} J$
**end for**

$$ \frac{d}{dt} \mathbf{a}(t) = \Omega(\mathbf{a}(t)) + \int_{0}^{t} K(\mathbf{a}(s))ds + F(t) \quad (31) $$

where $\mathbf{a}(t)$ are observables of a system defined as a projection of the full order state onto the observable space and $\Omega$ is the projected part of the original dynamics. These two quantities can respectively be identified as the reduced coordinates and dynamics of classical POD-Galerkin models. The remaining terms account for the impact of the non-observed coordinates of the system on the resolved dynamics. $K$ represents the dynamical exchanges between resolved and unresolved dynamics during the simulation, while $F$ accounts for the incomplete knowledge of the initial condition and system dynamics. Under the condition that the unresolved dynamics be dissipative, which is reasonable when the unobserved coordinates correspond to the small scales of a dynamical system, and that the boundary of the integral term in Eq. (31) be modified to $-\infty$, the last term in the Mori-Zwanzig formalism vanishes, leading to the following formulation:

$$ \frac{d}{dt} \mathbf{a}(t) = \Omega(\mathbf{a}(t)) + \int_{-\infty}^{t} K(\mathbf{a}(s))ds. \quad (32) $$

The residual model $\mathcal{R}$ proposed above can then be identified with the memory integral defined by the time-convolution kernel $K$, providing a strong connection with our non-Markovian correction hypothesis. Framing the CD-ROM architecture in the context of the Mori-Zwanzig formalism further justifies our modeling choices. Yet, it does not provide additional insights into the form of the correction model since little information is known about the convolution operator which can be infinite dimensional in certain cases. Choices about the structure of $K$ need to be made. In this work, assumptions are made about the vanishing impact of past states of the system on the residual model, accounted for by the matrix $\Lambda$ in our approach.
## 4.2. Deep Learning interpretation of the model

In this paragraph, the link between the CD-ROM architecture and classical deep learning models is discussed. Numerous methods have been proposed to model sequential data, each relying on a specific mechanism to extract and retain meaningful information from the past states of the system. The most popular architectures, GRU and LSTM ([59, 60, 61]), both use a combination of *gating* mechanisms to learn long term dependencies in a sequence. Other approaches like reservoir computing rely on an underlying dynamical system forced by the sequence data to predict the required output. In fact, strong similarities can be identified between our approach and *echo state networks* (ESN [62, 63]), a widespread reservoir computing architecture. ESNs are based on the simulation of random dynamics described by matrices $W_{in}$ and $W_{R}$. The matrix $W_{in}$ is used to encode some data $x_{t}$ in higher dimension, while the matrix $W_{R}$ holds the weights of the reservoir used to advance the state (memory) $y_{t}$ in time:

$$y_{t+1} = \sigma(W_{in}x_{t} + W_{R}y_{t}). \qquad (33)$$

For the sake of comparison, the nonlinear activation function $\sigma$ is dropped, and we consider that equation (33) results from the Euler integration with time-step $\Delta_{t}$ of a continuous system describing the dynamics of the memory $\frac{dy}{dt}$:

$$y_{t+1} = W_{in}x_{t} + W_{R}y_{t} \qquad (34)$$

$$= y_t + \Delta_t \frac{dy}{dt} \qquad (35)$$

which leads to the following expression:

$$\frac{dy}{dt} = \frac{W_{in}}{\Delta_t} x_t + \frac{(W_R - I)}{\Delta_t} y_t. \qquad (36)$$

Because the spectral radius of the matrix $W_{R}$ is constrained to be less than unity, all the eigenvalues of the operator $\bar{W}_{R} = \frac{(W_{R}-I)}{\Delta_{t}}$ have negative real parts, leading the ESN to have memory dissipation properties similar to those of our model. To underline this similarity, we can express the state of the ESN at time $t$ by diagonalising the operator $\bar{W}_{R} = P\Lambda_{R}P^{-1}$:

$$y(t) = \int_{0}^{T} P e^{\Lambda_{R}(t-\tau)} P^{-1} \frac{W_{in}}{\Delta_{t}} x(\tau) d\tau. \qquad (37)$$

Thus, a parallel between our model and the ESN is outlined. Major differences remain in the fact that the dynamics ($W_{R}$) and encoding matrix ($W_{in}$) are not optimised during the ESN training. It should also be noted that non-linearity is introduced in the ESN dynamics through the $\sigma$ activation function, while our model is based on linear memory dynamics and on a nonlinear encoding operator. These comparisons between direct deep learning methods and the continuous correction approach help build intuition about the role of each term in the model. The encoder can be compared with the input gate of an
<figure><img src="image_6.png" /><figcaption>Figure 4: Computational domain for the cylinder case.</figcaption></figure>

LSTM, or the $W_{in}$ matrix of the ESN, the $\Lambda$ matrix provides a tunable *forget* mechanism, while the residual term $\mathcal{R}$ plays the role of the output operator in our “continuous recurrent network”.

## 5. Case presentation and reduced models

In this Section, we introduce the simulation cases selected to demonstrate the ability of the CD-ROM approach to improve the performance of POD-Galerkin models. We first present two flow problems to illustrate the benefits of the CD-ROM architecture in the context of fluid mechanics. The first case is the standard configuration of the flow over a cylinder, often used as a benchmark for reduction methods. The second case is the fluidic pinball flow, introduced in [64] for the development of new control strategies. Finally, we introduce the case of the 1D Kuramoto-Sivashinsky which we use to demonstrate the ability of the CD-ROM approach to extend to parametric simulation problems.

### 5.1. Flow over a cylinder

The two-dimensional incompressible flow over a cylinder has been extensively studied in the context of reduced order modeling [4, 5] and model identification [6, 7] which makes this test case a good initial benchmark for the proposed correction method. The flow is simulated at a Reynolds number of $Re = 100$ based on the cylinder diameter and the velocity of the incoming flow. In this regime, the flow is laminar and exhibits vortex shedding in the wake of the cylinder.

The flow is governed by the incompressible Navier-Stokes equations, here solved using the FEniCs finite elements solver [65, 66]. The retained mesh is shown in Figure 4. It is rectangle-shaped, spanning from $x = -5$ to $x = 15$ in the streamwise direction, and from $y = -5$ to $y = 5$ in the transverse direction. The cylinder has a diameter $D = 1$, centered around the origin. The inflow is modeled as a uniform axial flow ($u(-5, y) = [U_\infty, 0]$) and a free-slip condition is used for the lower and upper boundaries of the rectangle while a no-slip condition is enforced at the cylinder surface. Finally, a stress-free condition is used for the outlet.
<figure><img src="image_7.png" /><figcaption>Figure 5: Vorticity fields of the modes selected for the reduced order modeling of the cylinder flow.</figcaption></figure>

A reduced order modeling strategy is employed to obtain a baseline model for the correction approach. The main results are summarised here but we refer to [5] for details on the reduction strategy. The vortex shedding regime of the cylinder flow is simulated with the FE niCs solver to obtain snapshot data and compute the POD modes. The first two modes, accounting for more than 95% of the Frobenius norm of the snapshot matrix are selected. The steady solution of the system is computed with a Newton method to construct a so- called shift mode (a_Δ). This mode is computed as a vector orthogonal to the plane described by the first two modes, pointing to the base flow solution u_b, and serves as a support for the simulation of the transition of the system from its steady state to the vortex shedding limit cycle. A three-dimensional POD basis is thus finally obtained. They are shown in Fig. 5 in terms of the vorticity field.

The training data corresponds to the simulation of the transition of the system to the limit cycle of oscillations, starting from an initial condition u_0. This initial condition is taken as a point close to the base flow u_b which is the fixed point of the incompressible Navier Stokes equations. Following the procedure of Loiseau et al. ([7]), we chose

$$u_0 = u_b + \epsilon v_1, \qquad (38)$$

where v_1 is the first POD mode, and ϵ > 0 is a small coefficient used to perturb the unstable base flow.

Through Galerkin projection of the discretised Navier-Stokes equations, a system of 3 coupled ODEs is obtained, describing the dynamics of the reduced coordinates vector. The results of the simulation of the transition using both the Finite Elements model and the Galerkin ROM are displayed in Figure 6. Even though the three equation model is able to simulate the transient dynamics, its trajectory strongly diverges from the projected snapshot data. The transition starts much later than in the full order simulation, due to a growth rate of the ROM’s transition lower than what it should be. Another significant issue with the model is its stabilisation around the limit cycle (a_Δ = 0) where an overshoot can be observed before the vortex shedding regime is established, which is not observed in the snapshot data. These discrepancies between the two trajectories can be attributed to the ignored residual term in the dynamics, making this model a good baseline for the correction approach illustrated in Section 6.
<figure><img src="image_8.png" /><figcaption>Figure 6: True trajectory and simulation of an uncorrected Galerkin ROM. Plots (1) & (2) show the phase space trajectories projected on the $a_1 - a_Δ$ and $a_1 - a_2$ planes, while (3) shows the time evolution of the shift mode coefficient.</figcaption></figure>

A test trajectory is also simulated in FEniCs by taking a random initial condition in the phase space spanned by the 3 selected modes, such that:

$$ \mathbf{u}_{0,test} = \sum_{i=0}^{3} \mathbf{v}_i a_i, \quad a_i \sim \mathcal{N}(0, 1), \qquad (39) $$

where $\mathbf{v}_i$ are the POD modes, and $a_i$ are random reduced coordinates sampled from a normal distribution. Starting from this initial condition, the finite element model is simulated in FEniCs for 700 seconds to ensure the system reaches the oscillation regime. The performance of the various models on this data trajectory are presented in section 6.

## 5.2. Fluidic pinball

The second case used to demonstrate the approach is the so-called fluidic pinball. Initially proposed as a challenging test bed for the development of control laws [64], the fluidic pinball case offers a good trade-off between complexity of its dynamics and interpretability [67, 68]. The simulation domain (Figure 7) is composed of three equidistant cylinders, each generating vortices in the wake which interact to create rich dynamics. The mesh used for the simulation was provided by the authors of [69]. Displayed in Figure 7, the domain is a rectangle spanning from $x = -6$ to $x = 20$ in the streamwise direction, and $y = -6$ to $y = 6$ in the transverse direction. Three identical cylinders with diameter $D = 1$ are arranged in an equilateral triangle, with centers' coordinates $(0, 0.75)$, $(0, -0.75)$ and $(-1.25, 0)$ respectively. The boundary conditions are identical to those of the cylinder case in Sec. 5.1. The inflow is modeled as a uniform axial flow, the upper and lower boundaries of the computational domain are modeled as free-slip, while a no-slip condition is used for the walls of the three cylinders and the outlet is modeled as stress-free.

The flow is simulated at a Reynolds number of $Re = 130$. At this Reynolds number, the flow dynamics have been shown to be chaotic [68], which makes it a challenging problem for any correction model as the smallest error on the dynamics will make the simulated trajectory diverge exponentially fast in time from the truth. These kind of chaotic problems are starting to get traction as interesting benchmarks for forecasting tasks and modeling problems [70].
<figure><img src="image_9.png" /><figcaption>Figure 7: Computational domain of the fluidic pinball case.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Figure 8: Flow complexity of the pinball case. Left: Vorticity field of the flow in the chaotic regime. Right: Spectrum of the snapshot matrix.</figcaption></figure>

The flow is simulated for 1800 seconds in the chaotic regime, which, based on the ergodic property of the system, yields a trajectory long enough to be representative of it's attractor.

It should be noted that the pinball flow is much more complex than the cylinder, as evidenced from Figure 8 where it is seen that many POD modes are required to account for a significant part of the energy. Almost a thousand modes would be required to capture 99% of the Frobenius norm of the snapshot data, while only 8 are required to achieve the same accuracy in the cylinder case.

We chose to build a POD-Galerkin model of this flow using only the first 10 POD modes. Although this choice is somewhat arbitrary, it was made to challenge the correction method, as the mean reconstruction error of about 60% has an important impact on the approximated dynamics. Indeed, the obtained reduced model quickly separates from the original trajectory, as expected from a chaotic system. More problematic is the fact that the reduced model is very unstable and diverges after 70s of simulation, as shown in Figure 9. Application of the correction method is presented in Section 6, a discussion on the impact of the number of modes is provided in Appendix B, while additional information on the computational cost of the approach is given in Appendix C.
<figure><img src="image_11.png" /><figcaption>Figure 9: From top to bottom: time evolution of the first 3 POD coefficients. Black line: true value of the pinball’s POD coefficients, orange line: trajectories obtained by integrating an uncorrected Galerkin reduced order model constructed using the first 10 POD modes.</figcaption></figure>

### 5.3. Parametric Kuramoto-Sivashinsky equation

Finally, to illustrate the ability of the proposed method to extend to parametric problems, we introduce the case of the 1D parametric Kuramoto-Sivashinsky (KS) equation. This case is often used to validate physical modeling methods as it is fairly inexpensive to simulate and presents non linear dynamics. Moreover, depending on the parameters used for the simulation, the dynamics become chaotic, making it significantly more challenging to forecasting approaches. The KS equation is formulated as follows:

$$ \frac{\partial u}{\partial t} = -\frac{1}{2} \nabla \cdot u^2 - \Delta u - \nu \Delta^2 u, \qquad (40) $$

$$ u(x + L, t) = u(x, t), \qquad (41) $$

$$ u(x, 0) = g(x), \qquad (42) $$

where L is the length of the 1D simulation domain, g is the initial condition and $\nu$ is a parameter that controls the degree of dissipativity of the system. Taking inspiration from [27], we propose to learn a corrected ROM for this problem under varying $\nu$ values. As in [27], we chose L = 22 and g is computed as the sum of the four leading Fourier modes with coefficients 0.06. The problem is discretised spatially on a basis of N = 513 Fourier modes, and integrated in time through the semi-implicit third order scheme from [71]. This choice of discretisation implies that the periodicity condition (equation 41) is satisfied by construction. The simulation is carried out for a duration of 50 seconds, using a time step ($\Delta t$) of 0.025s.

As in [27], the parameter $\nu$ is varied in the range [0.3, 1.5]. As mentioned in the previous paragraph, this parameter controls the degree of dissipation in the system, thus, low $\nu$ values lead to more chaotic dynamics and a harder model reduction task. This is represented on figure 10, which displays the differences between simulations carried out at the limits of the chosen parameter range.

To create a reduced model of the system, the simulation problem is solved for 25 parameter values selected within a log-linear range from $\nu = 0.3$ to $\nu = 1.5$. With this initial data, the proper orthogonal decomposition of the snapshot
<figure><img src="image_12.png" /><figcaption>Figure 10: Left: Simulations of the KS equation carried out under two different parameter values. Right: Relative snapshot reconstruction accuracy against the number of POD modes used for reconstruction.</figcaption></figure>

matrix is computed. The evolution of the relative snapshot reconstruction error depending on the number of selected modes is shown in Figure 10. To assemble the Galerkin reduced model, we select the 25 leading POD modes, which account for more than 99.5% of the information in the snapshot data. The computed POD modes form a basis $V$, which can be used to approximate the solution field $u(t, \nu)$ computed for a given time and parameter value as follows: $u(t, \nu) \approx \bar{u}(t, \nu) = V\mathbf{a}(t, \nu)$.

As in the non-parametric case, computing the vector of reduced coordinates $\mathbf{a}(t, \nu)$ is sufficient to fully determine the approximate solution $\bar{u}(t, \nu)$. Finally, the Galerkin projection method described in section 2 is applied to the KS equation, yielding the following reduced model:

$$ \frac{d\mathbf{a}}{dt} = -\frac{1}{2}\mathbf{a}^\top \tilde{\mathcal{Q}}\mathbf{a} - \tilde{\mathcal{L}}\mathbf{a} - \nu\tilde{\mathcal{L}}^2\mathbf{a}, \quad (43) $$

$$ \tilde{\mathcal{Q}}_{i,j,k} = \langle v_i, \nabla(v_j v_k) \rangle, \quad i, j, k = 1, \dots, r \quad (44) $$

$$ \tilde{\mathcal{L}}_{i,j} = \langle v_i, \Delta v_j \rangle, \quad i, j = 1, \dots, r \quad (45) $$

$$ \tilde{\mathcal{L}}_{i,j}^2 = \langle v_i, \Delta^2 v_j \rangle, \quad i, j = 1, \dots, r \quad (46) $$

where $v_i$ are the POD modes, and $\langle \cdot, \cdot \rangle$ is an inner product defined over the computational domain. It can be noted that the $\tilde{\mathcal{Q}}$ and $\tilde{\mathcal{L}}$ operators are the one dimensional equivalent of the reduced Navier-Stokes operators introduced in section 2. Similarly, the operator $\tilde{\mathcal{L}}^2$ is a linear operator corresponding to the fourth order derivative in equation 40. To test the model, 62 test parameter values are selected randomly in the range $\nu \in [0.3, 1.5]$ following a log-uniform distribution.

The uncorrected Galerkin ROM is simulated using the semi-implicit scheme from [71]. Figure 11 presents the results obtained by simulating the Galerkin model under different parameter values. The figure clearly underlines the difficulty of modeling lower $\nu$ values, as we observe that the Galerkin model diverges
<figure><img src="image_13.png" /><figcaption>Figure 11: Coefficients of the first 4 POD modes simulated with the uncorrected Galerkin ROM under different parameter values. Dashed line: Projected DNS data, full line: Simulation of uncorrected the Galerkin ROM.</figcaption></figure>

quickly from the true trajectory. Results obtained by augmenting the Galerkin ROM with the CD-ROM architecture are presented in section 7.

# 6. Results and discussion

In this section, we present the results obtained by applying the CD-ROM method to the imperfect Galerkin ROMs presented in the previous section. While some of the design choices regarding each specific cases are discussed in the following paragraphs, we refer the reader to Appendix D for a description of the various training details and hyper-parameter choices.

## 6.1. *Cylinder case*

The reduced model derived in Section 5 above was shown to be efficient for the simulation of the vortex shedding regime, but not suited to the simulation of transient dynamics. To apply the proposed correction procedure to this model, the modeling terms introduced in Sec. 3 are added to the reduced model. The dimension of the memory is chosen to be ten times the dimension of the reduced state. The residual $\mathcal{R}$ and encoder $\mathbf{E}$ are defined as multi layer perceptrons, using the Rectified Linear Unit activation function. Finally, the diagonal of the memory matrix $\mathbf{\Lambda}$ is initialised at random from a normal distribution. The model is trained with the Adam optimizer, and the trajectory loss introduced earlier (Eq. (27)) is used in combination with additional regularization terms. Details on the regularization of the loss are provided in Appendix E.

As described in section 5, the model is trained on the trajectory data obtained by simulating the transition from an initial condition close to the base flow. Since this target trajectory starts close to the base flow of the system, which is stationary, the initial memory can be computed with minimal error through the following integral:
<figure><img src="image_14.png" /><figcaption>Figure 12: CD-ROM trajectory obtained at different levels of training. The first mode's amplitude is displayed in the top panel, while the third mode's amplitude is shown in the bottom panel.</figcaption></figure>

$$y(0) = E(a_b^*) \int_{-\tau_{\min}}^{0} e^{\Lambda \tau} d\tau$$

where $a_b^* = V_r^T u_b$ are the reduced coordinates of the base flow and $\tau_{\min}$ is the longest time horizon defined by the $\Lambda$ matrix as in Equation (29). Finally, the parameters of the models are progressively optimised to reproduce the true transition trajectory as shown in Figure 12.

The CD-ROM model is integrated in time using an adaptive RK-45 scheme. Simulation results on the training trajectory are presented on Figure 13. It can be seen that the corrected model follows the training trajectory perfectly, triggering the transition at the right time, and correcting the oscillations of the original ROM during the stabilisation on the limit cycle. Moreover, the graph shows that the correction applied by the residual model is strong during the transition, where the original ROM struggles, and becomes minimal during the rest of the trajectory.

Finally, we present results of the performance of the model on the test trajectory. The first 8 seconds of DNS simulation are used to initialise the memory following Equation (30). Figure 14 presents the performance of each model on this trajectory. Because the starting point is not close to the base flow, the uncorrected ROM instantly exhibits transient dynamics, however, the growth rate of the transition is still inaccurate and the shift mode's trajectory presents the same non-physical oscillations around the limit cycle observed in the training trajectory. The corrected model does much better than the original ROM, simulating a more accurate transition, and stabilising on the limit cycle almost perfectly.

## 6.2. Fluidic pinball results

While the cylinder case discussed above offers a simple test bed for the presentation of the approach and its potential, it has already been shown that very
<figure><img src="image_15.png" /><figcaption>Figure 13: Results obtained with the corrected ROM for the simulation of the cylinder flow's transient dynamics. Left: phase space trajectory from base flow to limit cycle simulated with the finite element solver, the original Galerkin ROM and our corrected model. Right: time evolution of the shift mode's coefficient simulated with the same models, as well as the norm of the correction applied by our model.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Figure 14: Simulation results obtained on a test trajectory. As in Figure 13, the first graph represents the phase space trajectory, while the time evolution of the shift mode's coefficient is presented on the second graph.</figcaption></figure>
parsimonious models could be used to model its dynamics [7], which suggests
this configuration might not require a high dimensional neural network to learn
a correction term. The case of the fluidic pinball is more challenging and can
better underline the ability of our method to handle complex physics. As pre-
sented in the previous section, the POD-Galerkin approach is not well suited to
the reduction of the pinball case. The number of modes required to reconstruct
the snapshot data with a satisfying accuracy is very large and using a small
number of modes leads to an unstable model.

To apply the correction approach to the pinball case, a correction model
is built. The encoder and residual models are multi layer perceptrons and
the Sigmoid Linear Unit activation function is used as it leads to smoother
integration. The weights of the residual model are initialised to be close to 0 so
that the ROM is initially almost uncorrected. Diagonal entries of the memory
matrix are initialised as a log-linear range of time horizons, ranging from 0.6 to
3.84 seconds. The training data consists of 1800 seconds (15000 snapshots) of
DNS simulation in the chaotic regime. The leading two thirds of the simulated
DNS trajectory are used for training while the remaining third is set aside for
testing.

As presented in the previous section, the uncorrected 10-mode ROM is un-
stable and diverges after 70 s of simulation. As a result, trying to optimise
the correction model for long trajectories directly would lead to a very unstable
training process. To address this issue and stabilise the training, the corrected
ROM is trained on sub-trajectories of only 10 seconds This length is chosen as it
is short enough for the model to remain stable and long enough for the impact of
the encoder on the memory to be accounted for. Once a good correction model
has been trained for 10 seconds long trajectories, the length is progressively
augmented to attain the target horizon of 120 s. Besides the stabilisation of the
training, using sub-trajectories also allows for parallel training. Multiple sub-
trajectories can be sampled from the snapshot data and integrated in parallel on
a GPU to dramatically speed-up training, more information about the training
strategy as well as the training parameters used is provided in Appendix D.

The pinball correction model was trained to follow true trajectories up to
120 seconds. Results of the simulated flow fields are presented on Figure 15.
It can be seen that the projection on the 10-dimensional POD basis effectively
filters part of the spatial structures, leading to the divergence of the uncor-
rected reduced model. In contrast, the corrected ROM is able to reproduce the
projected flow field accurately.

As with the cylinder case, the CD-ROM is simulated with an adaptive RK-
45 scheme. Trajectory results simulated from a condition in the training basis
are presented on Figure 16. The model starts quickly diverging from the true
trajectory after the training horizon (120 s) as is expected from the chaotic
nature of the problem. Deriving a model to perfectly follow the DNS trajectory
would here make little modeling sense.

More interesting is the fact that, despite leaving the training trajectory, the
corrected model does not become unstable, even when integrated for over 1000
seconds with an initial condition outside of the training data. This suggests
<figure><img src="image_17.png" /><figcaption>Figure 15: Vorticity fields obtained after 110 s of simulation in the chaotic regime.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Figure 16: Results of the correction approach applied to the 10-mode ROM of the fluidic pinball. The initial condition is a random point selected in the training trajectory. Plots (a), (b), and (c) describe the coefficients of the first three POD modes simulated with different models. Plot (d) presents the relative Euclidean distance between the true data and the trajectory simulated with the CD-ROM.</figcaption></figure>

that the dynamics correction learned by the model has some physical consistency, dissipating the necessary energy which would otherwise have made the simulation diverge in the uncorrected case.

The intuition that the model was able to learn a Physics-compatible correction is confirmed when looking at the statistics of the attractor spanned by the CD-ROM's trajectory. Using the *nolitsa* library [72], the correlation dimension ([73]) was estimated, as well as the maximum Lyapunov exponent ([74]) of the corrected and true trajectories. The results are shown in Figure 17 where the model is seen to reproduce well the characteristics of the true attractor. One can also look at the probability distributions of the mode's amplitudes, presented on Figure 18. Once again, the simulated trajectory reproduces the results of the true simulation.

Note that long trajectories (several hundreds of seconds) were simulated with the corrected model to obtain these statistics. Due to the chaotic nature of the problem and the length of the integration period, the model has visited parts of the attractor different from those seen in the training trajectory. Despite this,
<figure><img src="image_19.png" /><figcaption>Figure 17: Statistics of the 10-dimensional attractors described by the projected DNS data (black) and the CD-ROM (blue). Left: Average number of points $g$ in a sphere of radius $r$, the data follows the law $g = r^{C_d}$ where $C_d$ corresponds to the correlation dimension of the attractor. Right: Time evolution of the average distance $d$ between trajectories starting from arbitrarily close initial conditions on the attractor. The data follows the law $d = e^{tl_e}$ where $l_e$ corresponds to the maximum lyapunov exponent of the system at hand.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>Figure 18: Estimated probability density functions for the coefficient of each mode. Plain blue line: Statistics of the trajectory simulated with the corrected model; dotted black line: Statistics of the projected DNS data. Labels refer to the mode index.</figcaption></figure>

the CD-ROM remains stable and describes an attractor with statistics that are very similar to those of the true attractor, further supporting the approach for the reduced modeling of complex dynamics.

### 6.3. Parametric KS equation results

In this paragraph, we present the results obtained by applying the CD-ROM architecture to a parametric case, the KS equation presented in section 5.3. To do so, the reduced Galerkin model of equation 43 is augmented with the proposed CD-ROM architecture (equation 18). To account for the parametric nature of the problem, the coefficient $\nu$ is passed to both the residual ($\mathcal{R}$) and encoder ($\mathbf{E}$) models, yielding the following CD-ROM system:

$$ \begin{aligned} \frac{d}{dt}\mathbf{a} &= -\frac{1}{2}\mathbf{a}^\top\tilde{\mathcal{Q}}\mathbf{a} - \tilde{\mathcal{L}}\mathbf{a} - \nu\tilde{\mathcal{L}}^2\mathbf{a} + \mathcal{R}(\mathbf{y}, \nu; \boldsymbol{\theta}_R), \\ \frac{d}{dt}\mathbf{y} &= \mathbf{E}(\mathbf{a}, \nu; \boldsymbol{\theta}_E) - \mathbf{\Lambda}\mathbf{y}. \end{aligned} \quad (47) $$
<figure><img src="image_21.png" /><figcaption>Figure 19: Coefficients of the first 4 POD modes simulated with the CD-ROM as well as the uncorrected Galerkin ROM under different parameter values in the test dataset. Dashed line: Projected DNS data, full orange line: Simulation of the Galerkin ROM, full blue line: Simulation of the CD-ROM.</figcaption></figure>

The residual and encoder models are both expressed as multi layer perceptrons, using the SiLU activation function. The weights of both neural networks, as well as the memory matrix $\Lambda$, are optimised using the Adam optimizer. As in the previous fluidic pinball case, we start by optimising the model on small sub-trajectories, then gradually increase the length of the sub-trajectories as the model reaches the desired accuracy. The model is trained on the data generated to compute the POD modes and assemble the Galerkin ROM (see Section 5.3). This training data corresponds to simulations carried out under 25 different parameter values in the range [0.3,1.5]. After training, the model is tested on 62 new simulations carried out under different parameter values selected randomly following a log-uniform distribution, as described in Section 5.3.

Contrary to the two previous flow cases, it is more efficient to simulate the KS equations using the semi implicit time-stepping scheme of [71], thus we use this scheme to integrate the CD-ROM model in time. Figure 19 presents the results obtained by simulating the CD-ROM using $\nu$ values not included in the training data. On the two cases presented, the CD-ROM performs better than the Baseline Galerkin ROM. It can be seen that the model diverges from the true trajectory data earlier in the case where the value of $\nu$ is lower. This is expected as we showed in section 5.3 that lower $\nu$ values lead to more complex dynamics.

To assess the performance of the CD-ROM over the whole test set, the Euclidean distance between the simulated reduced coordinates vector $\mathbf{a}$ and the projected DNS data $\mathbf{a}^*$ is computed over time, and for every test parameter values:

$$d(t, \nu) = \|\mathbf{a}(t, \nu) - \mathbf{a}^*(t, \nu)\|_2. \qquad (48)$$

Figure 20 presents the values of the error metric (equation 48) for every test
<figure><img src="image_22.png" /><figcaption>Figure 20: Error metric (equation 48) computed for each test parameter value.</figcaption></figure>

parameter value at select time steps. The figure shows that the CD-ROM is able to remain significantly closer to the true trajectory than the baseline Galerkin ROM for almost 20 seconds. The model then behaves differently depending on the parameter value. Cases in the $v \in [0.3, 0.5]$ range presenting the more chaotic dynamics quickly diverge from the true trajectory, while the CD-ROM is able to beat the baseline on the rest of the test cases for up to 50 seconds.

These results demonstrate the ability of the CD-ROM architecture to improve the Galerkin model order reduction approach in a parametric setting. The trained CD-ROM model is able to reproduce the dynamics of the full order system better than its uncorrected counterpart, even when using parameter values different from the training conditions.

## 7. Conclusion

This paper was concerned with an augmented reduced order modeling strategy based on the hybridisation of the classical Galerkin projection method and simple neural networks. By studying the limitations of the Galerkin projection, we established links with the field of partially observed systems, which led us to use past observables of the studied system as a critical ingredient for the correction of Galerkin models. Building on this result, the CD-ROM architecture was proposed to extract and exploit useful information from the system trajectory, by embedding the model with a delay differential equation structure. Moreover, the training strategy based on adjoint optimisation ensures *a-posteriori* performance of the model on the training trajectory.

The CD-ROM approach was demonstrated on two simple CFD test cases, namely the flow over a cylinder and the fluidic pinball. Numerical experiments have shown that the corrected models were able to capture the true dynamics with a high degree of accuracy, reproducing the true transition in the case of the
cylinder flow, and following the training trajectory for multiple Lyapunov times
in the fluidic pinball case. Moreover, these experiments outlined the reliability of
the corrected model as it performed better than its uncorrected counterpart even
outside the training conditions. The results obtained with the fluidic pinball
are particularly encouraging. We showed that the correction model was able
to stabilise the original Galerkin model in a consistent physical manner, as the
attractor simulated with the CD-ROM approach presents statistics similar to
the original attractor.

The ability of the proposed approach to extend to parametric problems
was also demonstrated. The CD-ROM method was applied to the case of the
Kuramoto-Sivashinsky equation with varying viscosity. After training the model
on a small number of parameter values in a selected range, we showed that the
CD-ROM approach improved the performance of the baseline Galerkin ROM
over the whole parameter range, even when simulating using parameter values
outside of the training data. This is of interest for many real-world situations,
e.g. industrial applications where a low cost parametric model is a key-enabler.

On-going work investigates the application of the CD-ROM approach on a
more concrete industrial use case, to demonstrate the capabilities of the method
on a real-world application. We are also looking into the best way to build para-
metric models of dynamical systems using the CD-ROM, as recent approaches
leveraging hyper-networks to dynamically modify the model weights have shown
promising results on data driven forecasting tasks.
# Appendix A. Neural ODEs and Adaptive Checkpoint Adjoint

## Appendix A.1. Adjoint Backpropagation

The adjoint backpropagation algorithm on which Neural ODEs are based can be used to solve problems of the form<fnref n="2" />:

$$ 
\begin{array}{ll}
\min_{\theta} & J(x(T)) \\
\text{s.t.} & \dot{x} - g(x; \theta) = 0 \qquad (A.1) \\
& x(0) = x_0.
\end{array}
$$

Generally, these kind of problems are solved with a gradient descent method, which requires evaluation of the following gradient:

$$ \frac{dJ}{d\theta} = \left. \frac{\partial J}{\partial x} \frac{\partial x}{\partial \theta} \right|_T. \qquad (A.2) $$

Because the derivative of $\dot{x}(t) = g(x; \theta)$ is parameterised by $\theta$, there is an implicit relation $x = x(t; \theta)$. The estimation of the sensitivity $\frac{\partial x}{\partial \theta}|_T$ requires the consideration of the impact of $\theta$ on the whole time-integration, from $t = 0$ to $t = T$. The adjoint backpropagation is used to avoid computing this term, which is done by evaluating the sensitivity of the following Lagrangian:

$$ \mathcal{L} = J + \int_{0}^{T} \mu(t)(\dot{x}(t) - g(x; \theta))dt \qquad (A.3) $$

$$ \dot{x}(t) = g(x; \theta) \implies \frac{dJ}{d\theta} = \frac{d\mathcal{L}}{d\theta}. \qquad (A.4) $$

The vector of Lagrangian multiplier $\mu$ is a function of time. Distributing the product in the integral and integrating the first term by parts leads to an expression where the sensitivity $\frac{\partial x}{\partial \theta}|_T$ can be isolated:

$$ \mathcal{L} = J + [\mu x]_0^T - \int_0^T \dot{\mu}x(t) + \mu(t)g(x; \theta)dt \qquad (A.5) $$

$$ \Rightarrow \frac{d\mathcal{L}}{d\theta} = \left(\frac{\partial J}{\partial x} + \mu(T)\right) \left. \frac{\partial x}{\partial \theta} \right|_T - \mu(0) \left. \frac{\partial x}{\partial \theta} \right|_0^0 - \int_0^T \left(\dot{\mu} + \mu(t) \frac{\partial g}{\partial x}\right) \left. \frac{\partial x}{\partial \theta} \right|_t + \mu(t) \left. \frac{\partial g}{\partial \theta} \right|_{x=x(t)} dt. \qquad (A.6) $$

A so-called adjoint equation can be derived to avoid having to solve for the sensitivity $\frac{\partial x}{\partial \theta}$. Enforcing a vanishing variation of the Lagrangian wrt $x$ at optimality yields:

<footnote n="2">²This specific criterion and constraints choice are well suited to our approach. The adjoint backpropagation algorithm can be used to solve more complex problems; however, we chose to restrict the scope to this formulation to clarify the derivation.</footnote>
$$ \frac{d\mu}{dt} = -\mu \frac{\partial g}{\partial x} \Big|_{t} \qquad (A.7) \\ \mu(T) = -\frac{\partial J}{\partial x} \Big|_{T}. $$

Solving the adjoint equation (A.7) for the values of $\mu(t)$ leads most of the terms in Equation (A.6) to vanish, so that:

$$ \frac{d\mathcal{L}}{d\theta} = - \int_{0}^{T} \mu(t) \frac{\partial g}{\partial \theta} \Big|_{x=x(t)} dt. \qquad (A.8) $$

As mentioned earlier, both this integral and the adjoint equation can be easily evaluated if $\dot{x}(t)$ is approximated by a neural network as the required vector-Jacobian products $\mu \frac{\partial g}{\partial x}$ and $\mu(t) \frac{\partial g}{\partial \theta}$ can be easily computed in a deep learning framework.

## Appendix A.2. Adaptive checkpoint adjoint

The original NeuralODE paper ([19]) proposed integrating forward in time to obtain the initial condition of the adjoint equation evaluated from $x(T)$, while discarding intermediate values $x(t)$. This choice was made with the goal of reducing the memory footprint of the method. However, discarding the intermediate time-steps $x(t)$ means that they have to be recomputed during backpropagation to evaluate the adjoint ODE (A.7), as well as the gradient integral (A.8). This can be done as a single backward in time integration by concatenating the different quantities $(x, \mu, \frac{d\mathcal{L}}{d\theta})$ in a single state vector:

$$ z = \left[ x, \mu, \frac{d\mathcal{L}}{d\theta} \right] \qquad (A.9) $$

$$ \frac{dz}{dt} = \left[ g(x; \theta), -\mu \frac{\partial g}{\partial x}, -\mu \frac{\partial g}{\partial \theta} \right] \qquad (A.10) $$

$$ z(T) = \left[ x(T), -\frac{\partial J}{\partial x} \Big|_{T}, \frac{\partial J}{\partial x} \frac{\partial g}{\partial \theta} \Big|_{T} \right]. \qquad (A.11) $$

Not only does this increase the computational cost of the method, but it can also lead to erroneous gradients, as numerical errors during integration can lead to differences between the forward and backward trajectories for $x(t)$. This was observed in [75] and is illustrated on Figure A.21, which shows that, despite using the same parameters, the forward and backward time trajectories diverge due to numerical errors. To address this issue, one can use the Adaptive Checkpoint Adjoint method. This method retains the intermediate integration steps $(t, x(t))$ and simply evaluates the integrals for the adjoint $\mu$ and the gradient $\frac{dJ}{d\theta}$ at the time steps selected during the forward integration, using the forward trajectory $x(t)$.

Although this increases the memory footprint of the method, it limits the computational cost of the backward pass, as the time steps are already selected,
<figure><img src="image_23.png" /><figcaption>Figure A.21: Forward and backward time integration of the Van der Pol oscillator. Numerical errors during the integration lead to the divergence of the two trajectories which should be identical.</figcaption></figure>

and the trajectory $x(t)$ does not have to be integrated another time. Most importantly, this limits the numerical errors introduced by the integration schemes, which can have a significant impact on the training, especially in the case of chaotic time-derivatives.

# Appendix B. Number of modes for the fluidic pinball reduction

## Appendix B.1. Uncorrected models comparison

In Section 5, we present a 10-mode reduced order model of the fluidic pinball problem. The number of modes is chosen somewhat arbitrarily to challenge the method. The more modes are used to model the flow, the better the model will be at reproducing the true dynamics, reducing the complexity of the required residual term. In this Appendix, we provide more insights into the impact of the number of modes on the modeling problem to clarify the choice of using 10 modes to model the pinball flow.

Figure B.22 presents the performance of different reduced models of the pinball flow. It is clearly seen that increasing the number of modes is beneficial for the performance of the reduced models. The magnitude of the closure term, as well as the speed at which the reduced model diverges from the true trajectory, are reduced when the number of modes increases. This Figure also shows that a higher number of modes leads to more stable reduced models. However, further experiments showed that even well resolved models such as the one using 173 modes were not stable and would diverge in certain conditions.

## Appendix B.2. Training convergence

To add to the argument, a comparison of the training convergence between two models using different number of modes is discussed. Two models are built, using respectively 10 and 55 POD modes, and the correction architecture is trained using the same parameters, presented in Table B.2.
<figure><img src="image_24.png" /><figcaption>Figure B.22: Comparison of the performance of POD-Galerkin models at different degrees of reduction. Top: Simulated value for the amplitude of the first POD mode; center: Value of the closure term for the first POD mode computed on the training trajectory; bottom: Relative distance between simulated and true trajectories.</figcaption></figure>

<table><caption>Table B.2: Parameters, e.g. learning rates (lr), used for the training of the compared CD-ROM architectures.</caption><thead><tr><th>Memory Dimension</th><th>5 × POD dimension</th></tr></thead><tbody><tr><td>Residual lr</td><td>10<sup>-3</sup></td></tr><tr><td>Encoder lr</td><td>10<sup>-3</sup></td></tr><tr><td>Λ lr</td><td>2 10<sup>-4</sup></td></tr><tr><td>Optimiser</td><td>AdamW</td></tr></tbody></table>

In both cases, the same loss is used, combining the optimisation of the distance between simulated and true trajectory with the residual regularization introduced in Appendix E:

$$J = \frac{1}{r} \left( \frac{1}{n_t} \sum_{i=1}^{n_t} \| \mathbf{a}(i\Delta_t) - \mathbf{a}^*(i\Delta_t) \|_2^2 + c \frac{1}{n_t} \sum_{i=1}^{n_t} \| \mathcal{R}(\mathbf{y}(i\Delta_t)) - \dot{\mathbf{a}}^*(i\Delta_t) \|_2^2 \right). \quad (B.1)$$

where $c = 0.1$ is a constant weighting the importance of the stochastic residual regularization term w.r.t. the trajectory loss. Notice that, to ease the comparison between the two models, the loss is scaled by the number of modes used in the ROM. The two models are trained in the same fashion, sub-trajectories of a hundred time steps are sampled in the training base and the loss (B.1) is optimized until a chosen threshold ($5 \times 10^{-4}$) is reached, at which point the length of the sub-trajectories is increased by fifty time steps. This process is repeated until the model is able to reproduce sub-trajectories of a thousand time steps.

Figure B.23 presents the evolution of the loss for the two models, as expected, the 55-mode model is quicker to train as it first reaches the threshold in 1400 epochs, while it takes the 10-mode model more than twice the number of epochs to achieve the same performance. Similarly, the 55-mode model achieves the required precision on trajectories of a 1000 time steps in only 4700 epochs, which is again more than twice as fast as the 10-mode model.
<figure><img src="image_25.png" /><figcaption>Figure B.23: Loss evolution of two CD-ROM architectures using a different number of modes.</figcaption></figure>

These results confirm the interest of only using the first 10 modes to challenge our correction approach. The relatively high magnitude of the residual term to be learned, the instability embedded in the model and the low degree of resolution of the problem (only 42% of the snapshot information) all constitute significant complexities which could arise in real world applications. While we showed that using a higher number of modes would simplify the modeling problems, this constraining choice helps demonstrate the applicability of the CD-ROM method to challenging modeling problems.

## Appendix C. Computational cost

While the computational cost of the overall CD-ROM approach and the way it compares to classical full order methods will strongly vary with the nature of the problem it is applied to, the simulation software used as well as the available hardware, we provide some elements of comparison with the POD Galerkin method in this appendix. Focusing on the case of the fluidic pinball (see section 5.2), we distinguish several components of the computational costs entailed by the CD-ROM method:

### Appendix C.1. ROM assembly

Because most POD Galerkin models are often restricted to a very low number of modes, the cost of assembling the reduced model can often be overlooked. However, some problems might require a high number of POD modes to achieve a satisfactory resolution. For example, the fluidic pinball case requires up to a thousand modes to capture 99% of the snapshot information, which directly impacts the cost of assembling the reduced operators $\tilde{\mathcal{L}}$ and $\tilde{\mathcal{Q}}$ in equation 9. Specifically, the reduced advective operator $\tilde{\mathcal{Q}}$ requires the computation of $O(n^3)$ inner products, $n$ being the number of selected POD modes.

This leads to exploding ROM assembly costs as the number of POD modes grows higher. While the assembly remains a one-time, parallelisable operation, we observed that assembling a 250 modes ROM on a 50 cpu machine took more than a full day of computation. This underlines the interest of representing the
<figure><img src="image_26.png" /><figcaption>Figure C.24: Performance and computational cost of simulating different models starting from an initial condition outside the CD-ROM training basis.</figcaption></figure>

dynamics on a low dimensional basis of modes, as assembling a thousand mode ROM would become prohibitively expensive.

## Appendix C.2. ROM Simulation

Figure C.24 presents the comparative simulating costs and performance of different reduced models on the fluidic pinball case. The figure shows that, although the CD-ROM does diverge from the true trajectory after some time, it performs better than its uncorrected counterparts. Specifically, the uncorrected 173 modes reduced model which captures more than 90% of the snapshot information, diverges earlier than the CD-ROM, while being more expensive to simulate.

It can also be seen that the CD-ROM is significantly more expensive to simulate than the simple galerkin model. This can be explained by the cost of evaluating the neural networks embedded in the CD-ROM architecture. Indeed, neural networks require the evaluation of matrix vector products of relatively high dimension. In the pinball case, we use two hidden layers of 250 neurons for the correction model, which explains the computational cost increase. However, it is interesting to note that the cost of evaluating a multi layer perceptron scales quadratically with its width (number of neurons per layer), while the cost of evaluating the advection term in the galerkin ROM scales cubically with the number of modes. This explains the fact that the 173 modes Galerkin ROM is more expensive to simulate than the 10 dimensional CD-ROM, while having a lower accuracy. Similar to the previous paragraph, this shows the interest of correcting a low dimensional model, rather than simply increasing the number of modes.

## Appendix C.3. Notes on Full Order Models

Providing a detailed comparison with full order methods is outside of the scope of this study, as the computational cost of a full order model will depend on a large number of choices, ranging from simulation software implementation
to numerical integration choices. We can however state that in fluid mechanics examples, the full order models were extremely expensive to simulate when compared to the reduced models studied in this paper. For example, generating the snapshot and test data for the fluidic pinball case took more than a day on a 50 cores machine. By comparison, the simulation of the trained CD-ROM for the same duration is of the order of the second on a normal computer.

## Appendix D. Hyper-Parameters and training

This appendix details the various design and training choices made for the different models presented in the results section 6. As explained earlier in the paper, the models are trained using progressively longer prediction horizons. The main advantage of this method, as opposed to directly training with the target prediction horizon, is that potentially unstable systems in the earlier learning stages will struggle to reach long term prediction horizons without diverging, making the training extremely inefficient. Moreover, using this strategy means that a single, long DNS trajectory can be separated in numerous sub trajectories which can be batched together and simulated in parallel, yielding a very efficient training process. The procedure is presented in algorithm 2.

**Algorithm 2** *Batching and Sub-trajectories*

**Require:** $T_f, T_i$ the initial and final prediction horizons, $T_{inc}$ the horizon increment, $\mathcal{L}_t$ the loss threshold, $f(\mathbf{z}, \theta)$ a CD-ROM model

$$T \leftarrow T_i$$

$$\mathcal{L} \leftarrow +\infty$$

while $T < T_f$ do

while $\mathcal{L} > \mathcal{L}_t$ do

Sample a batch of trajectories $\mathbf{a}^*_{-\tau_{min} \to T}$

Compute the initial memory $\mathbf{y}_0$

$\mathbf{z}_0 \leftarrow [\mathbf{a}^*_0, \mathbf{y}_0]$

$[\mathbf{a}_t, \mathbf{y}_t] \leftarrow \mathbf{z}_0 + \int_0^t f(\mathbf{z}; \theta) d\mathbf{t}$

▷ CD-ROM simulation

$\mathcal{L} \leftarrow \mathcal{L}(\mathbf{a}_0 \to T, \mathbf{a}^*_{0 \to T})$

Backpropagation & Gradient Step

end while

$$T \leftarrow T + T_{inc}$$

Reorganise batches according to the new $T$

end while

Finally, the hyper-parameters values and training details for the different models trained using the above strategy are presented below. In an effort to improve readability, the various values are organised in table D.3.
<table><caption>Table D.3: Hyper parameters used in the training of the different models presented in section 6. The various parts of the CD-ROM are designated as follows, $\mathcal{R}$ the neural network predicting the residual of the Galerkin model, $\mathbf{E}$ the memory encoder model, $\mathbf{A}$ the diagonal memory matrix. The final line presents the order of magnitude of the training times of each model on single RTX 2080 gpu.</caption><thead><tr><th></th><th>Cylinder</th><th>Pinball</th><th>KS</th></tr></thead><tbody><tr><td>Memory Size</td><td>30</td><td>50</td><td>30</td></tr><tr><td>Corrector Neurons</td><td>(30,30,30,30,3)</td><td>(50,250,250,250,10)</td><td>(56,150,150,150,25)</td></tr><tr><td>Encoder Neurons</td><td>(3,9,15,21,27)</td><td>(10,17,24,31,40)</td><td>(26,21,16,11,5)</td></tr><tr><td>Activation</td><td>SiLU</td><td>SiLU</td><td>SiLU</td></tr><tr><td>Optimizer</td><td>Adam</td><td>Adam</td><td>Adam</td></tr><tr><td>Learning Rate (E & R)</td><td>10<sup>-3</sup></td><td>10<sup>-3</sup></td><td>5 &times; 10<sup>-4</sup></td></tr><tr><td>Weight Decay (E & R)</td><td>10<sup>-4</sup></td><td>10<sup>-2</sup></td><td>10<sup>-3</sup></td></tr><tr><td>Learning Rate (A)</td><td>10<sup>-4</sup></td><td>2 &times; 10<sup>-4</sup></td><td>5 &times; 10<sup>-4</sup></td></tr><tr><td>Weight Decay (A)</td><td>0</td><td>0</td><td>0</td></tr><tr><td>Time Integrator</td><td>Scipy RK-45</td><td>Scipy RK-45</td><td>Semi Implicit 3rd order [71]</td></tr><tr><td>Training Time</td><td>1/2 day</td><td>1 day</td><td>1/2 day</td></tr></tbody></table>

## Appendix E. Regularization

### Appendix E.1. Residual regularization

We observed that, only training the model to follow the trajectory data can lead to poor local optima with large corrections applied to the original model. Since the magnitude of the correction is meant to be small when $r$ is sufficiently large, this tends to indicate over-fitting. This is an issue as such a model does not capture the true dynamics, and will quickly diverge when evaluated on conditions different from the training trajectory. To address this, a regularization term can be added in the model, to limit the magnitude of the corrections applied to the ROM:

$$J = \frac{1}{n_t} \sum_{i=1}^{n_t} \left( \|\mathbf{a}(i \Delta_t) - \mathbf{a}^*(i \Delta_t)\|_2^2 + \rho \|\mathcal{R}(\mathbf{y}(i \Delta_t))\|_2^2 \right) \quad (\text{E.1})$$

where $n_t$ is the number of time steps of length $\Delta_t$ in the optimised trajectory, $\mathbf{a}^*$ is the training trajectory data and $\rho$ is a hyper-parameter chosen to balance the importance of the regularization w.r.t. the rest of the loss. Although this fairly simple approach already helped guiding the training, we obtained better results by adding a custom regularization loss. This loss is based on the approximate value of the residual on the true trajectory, which can be obtained by computing the time-derivative of the true reduced coordinates $\dot{\mathbf{a}}^*$ through finite differences, and computing the difference with the derivative defined by the uncorrected ROM:

$$\mathcal{R}^* = \dot{\mathbf{a}}^* - \tilde{\mathcal{L}}\mathbf{a}^* - \mathbf{a}^* \tilde{\mathcal{Q}}\mathbf{a}^*. \quad (\text{E.2})$$

With this correction data, the memory can be evaluated from the true trajectory integrating in time the following ODE:
$$ \frac{d}{dt} \mathbf{y} = \mathbf{E}(\mathbf{a}^*; \theta_E) - \mathbf{\Lambda} \mathbf{y}. \qquad (E.3) $$

One can then define a regularization term $\mathcal{L}_{corr}$ for the loss:

$$ \mathcal{L}_{corr} = \frac{1}{n_t} \sum_{i=1}^{n_t} \| \mathcal{R}(\mathbf{y}(i\Delta_t)) - \dot{\mathbf{a}}^*(i\Delta_t) \|_2^2. \qquad (E.4) $$

The regularized loss definition was observed to lead to models with better generalization properties. Note that $\mathcal{L}_{corr}$ can be computed for a small random subset of the training batch at each epoch, to accelerate training while keeping a “stochastic” regularization for the residual.

## Appendix E.2. Encoded space regularization

Regularizing the encoded space can both help smooth the training process and increase the robustness of the model to unseen conditions. Taking inspiration from existing work ([23]), we propose to add the identity function to the encoder model. This means that useful information is embedded in memory in the form of a time convolution of the past resolved states of the system. The encoded state $\mathbf{E}(\mathbf{a}_t)$ is then constituted of both the reduced state $\mathbf{a}_t$ and a learnable nonlinear transformation $\mathcal{MLP}(\cdot, \theta_E)$ of it:

$$ \mathbf{E}(\mathbf{a}_t, \theta_E) = [\mathbf{a}_t, \mathcal{MLP}(\mathbf{a}, \theta_E)] \qquad (E.5) $$

With this structure, the encoder will only be learning additional nonlinear transformations of the state, simplifying the training and introducing a level of structure in the encoded space as we ensure its first dimensions are coherent with the phase space.

## Appendix F. Memory Dimension

The dimension of the memory in the CD-ROM formulation is a hyper parameter that should be chosen depending on the case and the dimension of the reduced state. Choosing an excessively low memory dimension will lead to poor prediction performance, while using too high of a dimension will negatively impact the computational cost of the corrected model and might lead to overfitting and poor numerical conditioning.

In the present examples, good results were obtained using memory dimensions ranging from 2x to 10x the dimension of the reduced state. However, in cases where performance becomes critical and the memory dimension should be as small as possible, a simple mask can be added to the output of the encoder to isolate necessary memory dimensions:

$$ \tilde{\mathbf{E}}(x) = \mathbf{m}_{\theta} \odot \mathbf{E}(x) \qquad (F.1) $$

where $\odot$ is the Hadamard (pointwise) product and $\mathbf{m}_{\theta} \in \mathbb{R}^{n_E}$ a vector with entries $m_i \in [0, 1]$. By adding the $L_1$ norm ($\|\mathbf{m}_{\theta}\|_1$) of this vector to the training loss, the model learns a sparse mask $m_{\theta}$. This mask isolates redundant memory dimensions, which can then be pruned after training.
Acknowledgements

This work was supported by IRT SystemX, in the context of the *Intelligence Artificielle et Ingénierie Augmentée* program.

References

[1] J. L. Lumley, Stochastic tools in turbulence. volume 12. applied mathematics and mechanics, Tech. rep., Pennsylvania state university, dept of aerospace engineering (1970).

[2] C. W. Rowley, S. T. Dawson, Model reduction for flow analysis and control, Annual Review of Fluid Mechanics 49 (1) (2017) 387–417.
doi:10.1146/annurev-fluid-010816-060042.
URL https://doi.org/10.1146/annurev-fluid-010816-060042

[3] G. Berkooz, P. Holmes, J. Lumley, The proper orthogonal decomposition in the analysis of turbulent flows, Annual Review of Fluid Mechanics 25 (2003) 539–575. doi:10.1146/annurev.fl.25.010193.002543.

[4] A. E. Deane, I. G. Kevrekidis, G. E. Karniadakis, S. A. Orszag, Low-dimensional models for complex geometry flows: Application to grooved channels and circular cylinders, Physics of Fluids 3 (1991) 2337–2354.

[5] B. Noack, K. Afanasiev, M. Morzynski, G. Tadmor, F. Thiele, A hierarchy of low-dimensional models for the transient and post-transient cylinder wake, Journal of Fluid Mechanics 497 (2003) 335–363. doi:10.1017/S0022112003006694.

[6] S. L. Brunton, J. L. Proctor, J. N. Kutz, Discovering governing equations from data by sparse identification of nonlinear dynamical systems, Proceedings of the National Academy of Sciences 113 (15) (2016) 3932–3937.
doi:10.1073/pnas.1517384113.
URL https://www.pnas.org/content/113/15/3932

[7] J.-C. Loiseau, S. L. Brunton, Constrained sparse galerkin regression, Journal of Fluid Mechanics 838 (2018) 42–67.

[8] E. Kaiser, B. R. Noack, L. Cordier, A. Spohn, M. Segond, M. Abel, G. Daviller, J. Östh, S. Krajnović, R. K. Niven, Cluster-based reduced-order modelling of a mixing layer, Journal of Fluid Mechanics 754 (2014) 365–414. doi:10.1017/jfm.2014.355.

[9] I. Mezić, Analysis of fluid flows via spectral properties of the koopman operator, Annual Review of Fluid Mechanics 45 (1) (2013) 357–378. doi:10.1146/annurev-fluid-011212-140652.

[10] P. Schmid, J. Sesterhenn, Dynamic mode decomposition of numerical and experimental data, Journal of Fluid Mechanics 656 (11 2008). doi:10.1017/S0022112010001217.
[11] M. O. Williams, I. G. Kevrekidis, C. W. Rowley, A data-driven approximation of the koopman operator: Extending dynamic mode decomposition, Journal of Nonlinear Science 25 (6) (2015) 1307-1346. doi:10.1007/s00332-015-9258-5.

[12] I. G. K. Matthew O. Williams, Clarence W. Rowley, A kernel-based method for data-driven koopman spectral analysis, Journal of Computational Dynamics 2 (2) (2015) 247-265.

[13] T. Pfaff, M. Fortunato, A. Sanchez-Gonzalez, P. Battaglia, Learning mesh-based simulation with graph networks, in: International Conference on Learning Representations, 2021.

[14] M. Raissi, P. Perdikaris, G. 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 (2019) 686-707. doi:https://doi.org/10.1016/j.jcp.2018.10.045.

[15] J. Berg, K. Nyström, A unified deep ANN approach to PDEs in complex geometries, Neurocomputing 317 (2018) 28-41.

[16] M. Raissi, Deep hidden physics models: Deep learning of nonlinear partial differential equations, JMLR 19 (1) (2018) 932-955.

[17] J. Sirignano, K. Spiliopoulos, DGM: A deep learning algorithm for solving PDEs, Journal of computational physics 375 (2018) 1339-1364.

[18] S. L. Brunton, B. R. Noack, P. Koumoutsakos, Machine learning for fluid mechanics, Annual Review of Fluid Mechanics 52 (1) (2020) 477-508. doi:10.1146/annurev-fluid-010719-060214.

[19] R. T. Q. Chen, Y. Rubanova, J. Bettencourt, D. K. Duvenaud, Neural ordinary differential equations, in: S. Bengio, H. Wallach, H. Larochelle, K. Grauman, N. Cesa-Bianchi, R. Garnett (Eds.), Advances in Neural Information Processing Systems, Vol. 31, Curran Associates, Inc., 2018.

[20] E. Haber, L. Ruthotto, Stable architectures for deep neural networks, Inverse Problems 34 (1) (2017) 014004. doi:10.1088/1361-6420/aa9a90.

[21] S. Greydanus, M. Dzamba, J. Yosinski, Hamiltonian neural networks, in: H. Wallach, H. Larochelle, A. Beygelzimer, F. d'Alché-Buc, E. Fox, R. Garnett (Eds.), Advances in Neural Information Processing Systems, Vol. 32, Curran Associates, Inc., 2019.

[22] J. Behrmann, D. K. Duvenaud, J.-H. Jacobsen, Invertible residual networks, in: ICML, 2019.

[23] Q. Li, F. Dietrich, E. M. Bollt, 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 27 (10) (2017) 103111. doi:10.1063/1.4993854.
[24] S. Otto, C. Rowley, Linearly-recurrent autoencoder networks for learning dynamics, SIAM Journal on Applied Dynamical Systems 18 (12 2017). doi:10.1137/18M1177846.

[25] S. Pan, K. Duraisamy, Physics-informed probabilistic learning of linear embeddings of nonlinear dynamics with guaranteed stability, SIAM Journal on Applied Dynamical Systems 19 (1) (2020) 480–509. doi:10.1137/19m1267246.

[26] M. M. N. B. Erichson, M. Mahoney, Physics-informed autoencoders for lyapunov-stable fluid flow prediction, in: Machine Learning and the Physical Sciences Workshop, Conference on Neural Information Processing Systems, 2019. URL https://arxiv.org/abs/1905.10866

[27] Q. Wang, N. Ripamonti, J. S. Hesthaven, Recurrent neural network closure of parametric pod-galerkin reduced-order models based on the mori-zwanzig formalism, Journal of Computational Physics 410 (2020) 109402.

[28] S. E. Ahmed, S. Pawar, O. San, A. Rasheed, T. Iliescu, B. R. Noack, On closures for reduced order models—a spectrum of first-principle to machine-learned avenues, Physics of Fluids 33 (9) (2021) 091301.

[29] P. R. Vlachas, G. Arampatzis, C. Uhler, P. Koumoutsakos, Multiscale simulations of complex systems by learning their effective dynamics, Nature Machine Intelligence (2022). doi:10.1038/s42256-022-00464-w.

[30] R. Maulik, A. Mohan, B. Lusch, S. Madireddy, P. Balaprakash, D. Livescu, Time-series learning of latent-space dynamics for reduced-order model closure, Physica D: Nonlinear Phenomena 405 (2020) 132368. doi:10.1016/j.physd.2020.132368.

[31] S. Pawar, S. M. Rahman, H. Vaddireddy, O. San, A. Rasheed, P. Vedula, A deep learning enabler for nonintrusive reduced order modeling of fluid flows, Physics of Fluids 31 (8) (2019) 085101. doi:10.1063/1.5113494.

[32] S. Fresca, A. Manzoni, Real-time simulation of parameter-dependent fluid flows through deep learning-based reduced order models, Fluids 6 (7) (2021). doi:10.3390/fluids6070259.

[33] H. Gao, J.-X. Wang, M. Zahr, Non-intrusive model reduction of large-scale, nonlinear dynamical systems using deep learning, Physica D: Nonlinear Phenomena 412 (2020) 132614. doi:10.1016/j.physd.2020.132614.

[34] S. R. Bukka, R. Gupta, A. R. Magee, R. K. Jaiman, Assessment of unsteady flow predictions using hybrid deep learning based reduced-order models, Physics of Fluids 33 (1) (2021) 013601. doi:10.1063/5.0030137.
[35] R. Maulik, B. Lusch, P. Balaprakash, Reduced-order modeling of advection-dominated systems with recurrent neural networks and convolutional autoencoders, Physics of Fluids 33 (3) (2021) 037106. doi:10.1063/5.0039986.

[36] H. Eivazi, H. Veisi, M. H. Naderi, V. Esfahanian, Deep neural networks for nonlinear model order reduction of unsteady flows, Physics of Fluids 32 (10) (2020) 105104. doi:10.1063/5.0020526.

[37] P. Wu, J. Sun, X. Chang, W. Zhang, R. Arcucci, Y. Guo, C. C. Pain, Data-driven reduced order model with temporal convolutional neural network, Computer Methods in Applied Mechanics and Engineering 360 (2020) 112766. doi:https://doi.org/10.1016/j.cma.2019.112766.

[38] S. Pawar, S. E. Ahmed, O. San, A. Rasheed, Data-driven recovery of hidden physics in reduced order modeling of fluid flows, Physics of Fluids 32 (3) (2020) 036602. doi:10.1063/5.0002051.

[39] R. Zwanzig, Nonequilibrium statistical mechanics, Oxford university press, 2001.

[40] P. Vlachas, W. Byeon, Z. Y. Wan, T. Sapsis, P. Koumoutsakos, Data-driven forecasting of high-dimensional chaotic systems with long short-term memory networks, Proceedings of the Royal Society A 474. doi:https://doi.org/10.1098/rspa.2017.0844.

[41] C. Michelén Ströfer, H. Xiao, End-to-end differentiable learning of turbulence models from indirect observations (04 2021).

[42] K. Um, R. Brand, Y. R. Fei, P. Holl, N. Thuerey, Solver-in-the-loop: Learning from differentiable physics to interact with iterative pde-solvers, in: H. Larochelle, M. Ranzato, R. Hadsell, M. Balcan, H. Lin (Eds.), Advances in Neural Information Processing Systems, Vol. 33, Curran Associates, Inc., 2020, pp. 6111–6122.

[43] P. Holmes, J. L. Lumley, G. Berkooz, C. W. Rowley, Turbulence, coherent structures, dynamical systems and symmetry, Cambridge university press, 2012.

[44] T. Lassila, A. Manzoni, A. Quarteroni, G. Rozza, Model order reduction in fluid dynamics: challenges and perspectives, Reduced Order Methods for modeling and computational reduction (2014) 235–273.

[45] S. Chaturantabut, D. C. Sorensen, Discrete empirical interpolation for nonlinear model reduction (2009) 4316–4321 doi:10.1109/CDC.2009.5400045.

[46] K. Carlberg, M. Barone, H. Antil, Galerkin v. least-squares petrov-galerkin projection in nonlinear model reduction, Journal of Computational Physics 330 (2017) 693–734. doi:https://doi.org/10.1016/j.jcp.2016.10.033.
[47] F. Takens, Detecting strange attractors in turbulence, in: D. Rand, L.-S. Young (Eds.), Dynamical Systems and Turbulence, Warwick 1980, Springer Berlin Heidelberg, Berlin, Heidelberg, 1981, pp. 366–381.

[48] J. L. Callaham, S. L. Brunton, J.-C. Loiseau, On the role of nonlinear correlations in reduced-order modeling, arXiv preprint arXiv:2106.02409 (2021).

[49] J. Dormand, P. Prince, A family of embedded runge-kutta formulae, Journal of Computational and Applied Mathematics 6 (1) (1980) 19–26. doi:https://doi.org/10.1016/0771-050X(80)90013-3.

[50] K. Gopalsamy, Stability and Oscillations in Delay Differential Equations of Population Dynamics, Kluwer Academic Publishers, 1992.

[51] E. E. Salpeter, S. R. Salpeter, Mathematical model for the epidemiology of tuberculosis, with estimates of the reproductive number and infection-delay function, American Journal of Epidemiology 147 (4) (1998) 398–406. doi:10.1093/oxfordjournals.aje.a009463.

[52] C. Pruett, T. Gatski, C. E. Grosch, W. Thacker, The temporally filtered navier–stokes equations: properties of the residual stress, Physics of Fluids 15 (8) (2003) 2127–2140.

[53] E. Åkervik, L. Brandt, D. S. Henningson, J. Hœpffner, O. Marxen, P. Schlatter, Steady solutions of the navier-stokes equations by selective frequency damping, Physics of fluids 18 (6) (2006) 068102.

[54] M. Leshno, V. Y. Lin, A. Pinkus, S. Schocken, Multilayer feed-forward networks with a nonpolynomial activation function can approximate any function, Neural Networks 6 (6) (1993) 861–867. doi:https://doi.org/10.1016/S0893-6080(05)80131-5.

[55] K. He, X. Zhang, S. Ren, J. Sun, Deep residual learning for image recognition, in: 2016 IEEE Conference on Computer Vision and Pattern Recognition (CVPR), 2016, pp. 770–778. doi:10.1109/CVPR.2016.90.

[56] L. Pontriagin, G. Boltjanskij, V. Boltanskij, K. M. R. Collection, V. Gamkrelidze, K. Trirogoff, L. Neustadt, R. Gamkrelidze, W. Neustadt, E. Misenko, et al., The Mathematical Theory of Optimal Processes, Interscience publishers, Interscience Publishers, 1962.

[57] H. Mori, Transport, Collective Motion, and Brownian Motion, Progress of Theoretical Physics 33 (3) (1965) 423–455. doi:10.1143/PTP.33.423.

[58] R. Zwanzig, Nonlinear generalized langevin equations, Journal of Statistical Physics (1973). URL https://doi.org/10.1007/BF01008729
[59] K. Cho, B. Merrienboer, C. Gulcehre, F. Bougares, H. Schwenk, Y. Bengio, Learning phrase representations using rnn encoder-decoder for statistical machine translation (06 2014). doi:10.3115/v1/D14-1179.

[60] S. Hochreiter, J. Schmidhuber, Long short-term memory, Neural computation 9 (1997) 1735-80. doi:10.1162/neco.1997.9.8.1735.

[61] J. Chung, C. Gulcehre, K. Cho, Y. Bengio, Empirical evaluation of gated recurrent neural networks on sequence modeling (12 2014).

[62] M. Lukosevicius, H. Jaeger, Reservoir computing approaches to recurrent neural network training, Comput. Sci. Rev. 3 (2009) 127-149.

[63] N. Doan, W. Polifke, L. Magri, Physics-informed echo state networks, Journal of Computational Science 47 (2020) 101237. doi:10.1016/j.jocs.2020.101237.

[64] B. R. Noack, M. Morzyński, The fluidic pinball — a toolkit for multiple-input multiple-output flow control, Tech. rep., Chair of Virtual Engineering, Poznan University of Technology, Poland (2017).

[65] A. Logg, K. Mardal, G. Wells, Automated Solution of Differential Equations by the Finite Element Method, 2012. doi:https://doi.org/10.1007/978-3-642-23099-8.

[66] M. Alnaes, J. Blechta, J. Hake, A. Johansson, K. B., A. Logg, C. Richardson, J. Ring, M. Rognes, G. Wells, The FEniCS Project, 2015. doi:https://doi.org/10.11588/ans.2015.100.20553.

[67] N. Deng, B. R. Noack, M. Morzyński, L. R. Pastur, Low-order model for successive bifurcations of the fluidic pinball, Journal of Fluid Mechanics 884 (12 2019). doi:10.1017/jfm.2019.959.

[68] N. Deng, L. Pastur, M. Morzyński, B. Noack, Route to chaos in the fluidic pinball (04 2021).

[69] G. Y. Cornejo Maceda, B. R. Noack, F. Lusseyran, M. Morzynski, L. Pastur, N. Deng, Taming the fluidic pinball with artificial intelligence control, in: European Fluid Mechanics Conference, Vienne, Austria, 2018.

[70] W. Gilpin, Chaos as an interpretable benchmark for forecasting and data-driven modelling (2021). arXiv:2110.05266.

[71] S. K. Kar, A semi-implicit runge-kutta time-difference scheme for the two-dimensional shallow-water equations, Monthly Weather Review 134 (10) (2006) 2916-2926. doi:10.1175/MWR3214.1.

[72] M. Mannattil, Nolitsa, https://github.com/charlespwd/project-title (2018).
[73] P. Grassberger, I. Procaccia, Measuring the strangeness of strange attractors, Physica D: Nonlinear Phenomena 9 (1) (1983) 189-208.
doi:https://doi.org/10.1016/0167-2789(83)90298-1.

[74] A. Wolf, J. B. Swift, H. L. Swinney, J. A. Vastano, Determining lyapunov exponents from a time series, Physica D: Nonlinear Phenomena 16 (3) (1985) 285-317.
doi:https://doi.org/10.1016/0167-2789(85)90011-9.

[75] J. Zhuang, N. Dvornek, X. Li, S. Tatikonda, X. Papademetris, J. Duncan, Adaptive checkpoint adjoint method for gradient estimation in neural ODE, in: H. D. III, A. Singh (Eds.), Proceedings of the 37th International Conference on Machine Learning, Vol. 119 of Proceedings of Machine Learning Research, PMLR, 2020, pp. 11639-11649.
**Declaration of interests**

[x] The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

[ ] The authors declare the following financial interests/personal relationships which may be considered as potential competing interests: