# Augmenting Physical Models with Deep Networks for Complex Dynamics
Forecasting

Yuan Yin1 ∗Vincent Le Guen2,3 ∗Jérémie Dona1 ∗Emmanuel de Bézenac1
∗Ibrahim Ayed1,4 Nicolas Thome2 Patrick Gallinari1,5
1 Sorbonne Université, CNRS, LIP6, Paris, France
2 Conservatoire National des Arts et Métiers, CEDRIC, Paris, France
3 EDF R&D, Chatou, France
4 Theresis Lab, Thales
5 Criteo AI Lab, Paris, France Equal contribution, authors sorted by reverse
alphabetical order.

###### Abstract

Forecasting complex dynamical phenomena in settings where only partial
knowledge of their dynamics is available is a prevalent problem across various
scientific fields. While purely data-driven approaches are arguably
insufficient in this context, standard physical modeling based approaches tend
to be over-simplistic, inducing non-negligible errors. In this work, we
introduce the APHYNITY framework, a principled approach for augmenting
incomplete physical dynamics described by differential equations with deep
data-driven models. It consists in decomposing the dynamics into two
components: a physical component accounting for the dynamics for which we have
some prior knowledge, and a data-driven component accounting for errors of the
physical model. The learning problem is carefully formulated such that the
physical model explains as much of the data as possible, while the data-driven
component only describes information that cannot be captured by the physical
model, no more, no less. This not only provides the existence and uniqueness
for this decomposition, but also ensures interpretability and benefits
generalization. Experiments made on three important use cases, each
representative of a different family of phenomena, i.e. reaction-diffusion
equations, wave equations and the non-linear damped pendulum, show that
APHYNITY can efficiently leverage approximate physical models to accurately
forecast the evolution of the system and correctly identify relevant physical
parameters.

##  1 Introduction

Modeling and forecasting complex dynamical systems is a major challenge in
domains such as environment and climate (Rolnick et al., 2019), health science
(Choi et al., 2016), and in many industrial applications (Toubeau et al.,
2018). Model Based (MB) approaches typically rely on partial or ordinary
differential equations (PDE/ODE) and stem from a deep understanding of the
underlying physical phenomena. Machine learning (ML) and deep learning methods
are more prior agnostic yet have become state-of-the-art for several spatio-
temporal prediction tasks (Shi et al., 2015; Wang et al., 2018; Oreshkin et
al., 2020; Donà et al., 2020), and connections have been drawn between deep
architectures and numerical ODE solvers, e.g. neural ODEs (Chen et al., 2018;
Ayed et al., 2019b). However, modeling complex physical dynamics is still
beyond the scope of pure ML methods, which often cannot properly extrapolate
to new conditions as MB approaches do.

Combining the MB and ML paradigms is an emerging trend to develop the
interplay between the two paradigms. For example, Brunton et al. (2016); Long
et al. (2018b) learn the explicit form of PDEs directly from data, Raissi et
al. (2019); Sirignano & Spiliopoulos (2018) use NNs as implicit methods for
solving PDEs, Seo et al. (2020) learn spatial differences with a graph
network, Ummenhofer et al. (2020) introduce continuous convolutions for fluid
simulations, de Bézenac et al. (2018) learn the velocity field of an
advection-diffusion system, Greydanus et al. (2019); Chen et al. (2020)
enforce conservation laws in the network architecture or in the loss function.

The large majority of aforementioned MB/ML hybrid approaches assume that the
physical model adequately describes the observed dynamics. This assumption is,
however, commonly violated in practice. This may be due to various factors,
e.g. idealized assumptions and difficulty to explain processes from first
principles (Gentine et al., 2018), computational constraints prescribing a
fine grain modeling of the system (Ayed et al., 2019a), unknown external
factors, forces and sources which are present (Large & Yeager, 2004). In this
paper, we aim at leveraging prior dynamical ODE/PDE knowledge in situations
where this physical model is incomplete, i.e. unable to represent the whole
complexity of observed data. To handle this case, we introduce a principled
learning framework to Augment incomplete PHYsical models for ideNtIfying and
forecasTing complex dYnamics (APHYNITY). The rationale of APHYNITY,
illustrated in Figure 1 on the pendulum problem, is to augment the physical
model when—and only when—it falls short.

Designing a general method for combining MB and ML approaches is still a
widely open problem, and a clear problem formulation for the latter is lacking
(Reichstein et al., 2019). Our contributions towards these goals are the
following:

  * •

We introduce a simple yet principled framework for combining both approaches.
We decompose the data into a physical and a data-driven term such that the
data-driven component only models information that cannot be captured by the
physical model. We provide existence and uniqueness guarantees (Section 3.1)
for the decomposition given mild conditions, and show that this formulation
ensures interpretability and benefits generalization.

  * •

We propose a trajectory-based training formulation (Section 3.2) along with an
adaptive optimization scheme (Section 3.3) enabling end-to-end learning for
both physical and deep learning components. This allows APHYNITY to
automatically adjust the complexity of the neural network to different
approximation levels of the physical model, paving the way to flexible learned
hybrid models.

  * •

We demonstrate the generality of the approach on three use cases (reaction-
diffusion, wave equations and the pendulum) representative of different PDE
families (parabolic, hyperbolic), having a wide spectrum of application
domains, e.g. acoustics, electromagnetism, chemistry, biology, physics
(Section 4). We show that APHYNITY is able to achieve performances close to
complete physical models by augmenting incomplete ones, both in terms of
forecasting accuracy and physical parameter identification. Moreover, APHYNITY
can also be successfully extended to the partially observable setting (see
discussion in Section 5).

|  |
---|---|---
(a) Data-driven Neural ODE | (b) Simple physical model | (c) Our APHYNITY framework
Figure 1: Predicted dynamics for the damped pendulum vs. ground truth (GT)
trajectories
$\nicefrac{{\mathop{}\\!\mathrm{d}^{2}\theta}}{{\mathop{}\\!\mathrm{d}t^{2}}}+\omega_{0}^{2}\sin\theta+\alpha\nicefrac{{\mathop{}\\!\mathrm{d}\theta}}{{\mathop{}\\!\mathrm{d}t}}=0$.
We show that in (a) the data-driven approach (Chen et al., 2018) fails to
properly learn the dynamics due to the lack of training data, while in (b) an
ideal pendulum cannot take friction into account. The proposed APHYNITY shown
in (c) augments the over-simplified physical model in (b) with a data-driven
component. APHYNITY improves both forecasting (MSE) and parameter
identification (Error $T_{0}$) compared to (b).

##  2 Related work

#### Correction in data assimilation

Prediction under approximate physical models has been tackled by traditional
statistical calibration techniques, which often rely on Bayesian methods
(Pernot & Cailliez, 2017). Data assimilation techniques, e.g. the Kalman
filter (Kalman, 1960; Becker et al., 2019), 4D-var (Courtier et al., 1994),
prediction errors are modeled probabilistically and a correction using
observed data is applied after each prediction step. Similar residual
correction procedures are commonly used in robotics and optimal control (Chen,
2004; Li et al., 2014). However, these sequential (two-stage) procedures
prevent the cooperation between prediction and correction. Besides, in model-
based reinforcement learning, model deficiencies are typically handled by
considering only short-term rollouts (Janner et al., 2019) or by model
predictive control (Nagabandi et al., 2018). The originality of APHYNITY is to
leverage model-based prior knowledge by augmenting it with neurally
parametrized dynamics. It does so while ensuring optimal cooperation between
the prior model and the augmentation.

#### Augmented physical models

Combining physical models with machine learning (gray-box or hybrid modeling)
was first explored from the 1990’s: Psichogios & Ungar (1992); Thompson &
Kramer (1994); Rico-Martinez et al. (1994) use neural networks to predict the
unknown parameters of physical models. The challenge of proper MB/ML
cooperation was already raised as a limitation of gray-box approaches but not
addressed. Moreover these methods were evaluated on specific applications with
a residual targeted to the form of the equation. In the last few years, there
has been a renewed interest in deep hybrid models bridging data assimilation
techniques and machine learning to identify complex PDE parameters using
cautiously constrained forward model (Long et al., 2018b; de Bézenac et al.,
2018), as discussed in introduction. Recently, some approaches have
specifically targetted the MB/ML cooperation. HybridNet (Long et al., 2018a)
and PhICNet (Saha et al., 2020) both use data-driven networks to learn
additive perturbations or source terms to a given PDE. The former considers
the favorable context where the perturbations can be accessed, and the latter
the special case of additive noise on the input. Wang et al. (2019); Mehta et
al. (2020) propose several empirical fusion strategies with deep neural
networks but lack theoretical groundings. PhyDNet (Le Guen & Thome, 2020)
tackles augmentation in partially-observed settings, but with specific
recurrent architectures dedicated to video prediction. Crucially, all the
aforementioned approaches do not address the issues of uniqueness of the
decomposition or of proper cooperation for correct parameter identification.
Besides, we found experimentally that this vanilla cooperation is inferior to
the APHYNITY learning scheme in terms of forecasting and parameter
identification performances (see experiments in Section 4.2).

##  3 The APHYNITY Model

In the following, we study dynamics driven by an equation of the form:

$\frac{\mathop{}\\!\mathrm{d}X_{t}}{\mathop{}\\!\mathrm{d}t}=F(X_{t})$ (1)

defined over a finite time interval $[0,T]$, where the state $X$ is either
vector-valued, i.e. we have $X_{t}\in\mathbb{R}^{d}$ for every $t$, (pendulum
equations in Section 4), or $X_{t}$ is a $d$-dimensional vector field over a
spatial domain $\Omega\subset\mathbb{R}^{k}$, with $k\in\\{2,3\\}$, i.e.
$X_{t}(x)\in\mathbb{R}^{d}$ for every $(t,x)\in[0,T]\times\Omega$ (reaction-
diffusion and wave equations in Section 4). We suppose that we have access to
a set of observed trajectories
${\mathcal{D}}=\\{X_{\cdot}:[0,T]\rightarrow{\mathcal{A}}\ |\ \forall
t\in[0,T],\nicefrac{{\mathop{}\\!\mathrm{d}X_{t}}}{{\mathop{}\\!\mathrm{d}t}}=F(X_{t})\\}$,
where ${\mathcal{A}}$ is the set of $X$ values (either $\mathbb{R}^{d}$ or
vector field). In our case, the unknown $F$ has ${\mathcal{A}}$ as domain and
we only assume that $F\in{\mathcal{F}}$, with $({\mathcal{F}},\|\cdot\|)$ a
normed vector space.

###  3.1 Decomposing dynamics into physical and augmented terms

As introduced in Section 1, we consider the common situation where incomplete
information is available on the dynamics, under the form of a family of ODEs
or PDEs characterized by their temporal evolution
$F_{p}\in{\mathcal{F}}_{p}\subset{\mathcal{F}}$. The APHYNITY framework
leverages the knowledge of ${\mathcal{F}}_{p}$ while mitigating the
approximations induced by this simplified model through the combination of
physical and data-driven components. ${\mathcal{F}}$ being a vector space, we
can write:

$F=F_{p}+F_{a}$

where $F_{p}\in{\mathcal{F}}_{p}$ encodes the incomplete physical knowledge
and $F_{a}\in{\mathcal{F}}$ is the data-driven augmentation term complementing
$F_{p}$. The incomplete physical prior is supposed to belong to a known
family, but the physical parameters (e.g. propagation speed for the wave
equation) are unknown and need to be estimated from data. Both $F_{p}$ and
$F_{a}$ parameters are estimated by fitting the trajectories from
${\mathcal{D}}$.

The decomposition $F=F_{p}+F_{a}$ is in general not unique. For example, all
the dynamics could be captured by the $F_{a}$ component. This decomposition is
thus ill-defined, which hampers the interpretability and the extrapolation
abilities of the model. In other words, one wants the estimated parameters of
$F_{p}$ to be as close as possible to the true parameter values of the
physical model and $F_{a}$ to play only a complementary role w.r.t $F_{p}$, so
as to model only the information that cannot be captured by the physical
prior. For example, when $F\in{\mathcal{F}}_{p}$, the data can be fully
described by the physical model, and in this case it is sensible to desire
$F_{a}$ to be nullified; this is of central importance in a setting where one
wishes to identify physical quantities, and for the model to generalize and
extrapolate to new conditions. In a more general setting where the physical
model is incomplete, the action of $F_{a}$ on the dynamics, as measured
through its norm, should be as small as possible.

This general idea is embedded in the following optimization problem:

$\underset{F_{p}\in{\mathcal{F}}_{p},F_{a}\in{\mathcal{F}}}{\min}~{}~{}~{}\left\|F_{a}\right\|~{}~{}~{}\mathrm{subject~{}to}~{}~{}~{}~{}\forall
X\in{\mathcal{D}},\forall
t,\frac{\mathop{}\\!\mathrm{d}X_{t}}{\mathop{}\\!\mathrm{d}t}=(F_{p}+F_{a})(X_{t})$
(2)

The originality of APHYNITY is to leverage model-based prior knowledge by
augmenting it with neurally parametrized dynamics. It does so while ensuring
optimal cooperation between the prior model and the augmentation.

A first key question is whether the minimum in Eq. (2) is indeed well-defined,
in other words whether there exists indeed a decomposition with a minimal norm
$F_{a}$. The answer actually depends on the geometry of ${\mathcal{F}}_{p}$,
and is formulated in the following proposition proven in Appendix B:

######  Proposition 1 (Existence of a minimizing pair).

If ${\mathcal{F}}_{p}$ is a proximinal set111A proximinal set is one from
which every point of the space has at least one nearest point. A Chebyshev set
is one from which every point of the space has a unique nearest point. More
details in Appendix A., there exists a decomposition minimizing Eq. (2).

Proximinality is a mild condition which, as shown through the proof of the
proposition, cannot be weakened. It is a property verified by any boundedly
compact set. In particular, it is true for closed subsets of finite
dimensional spaces. However, if only existence is guaranteed, while forecasts
would be expected to be accurate, non-uniqueness of the decomposition would
hamper the interpretability of $F_{p}$ and this would mean that the identified
physical parameters are not uniquely determined.

It is then natural to ask under which conditions solving problem Eq. (2) leads
to a unique decomposition into a physical and a data-driven component. The
following result provides guarantees on the existence and uniqueness of the
decomposition under mild conditions. The proof is given in Appendix B:

######  Proposition 2 (Uniqueness of the minimizing pair).

If ${\mathcal{F}}_{p}$ is a Chebyshev set11footnotemark: 1, Eq. (2) admits a
unique minimizer. The $F_{p}$ in this minimizer pair is the metric projection
of the unknown $F$ onto ${\mathcal{F}}_{p}$.

The Chebyshev assumption condition is strictly stronger than proximinality but
is still quite mild and necessary. Indeed, in practice, many sets of interest
are Chebyshev, including all closed convex spaces in strict normed spaces and,
if ${\mathcal{F}}=L^{2}$, ${\mathcal{F}}_{p}$ can be any closed convex set,
including all finite dimensional subspaces. In particular, all examples
considered in the experiments are Chebyshev sets.

Propositions 1 and 2 provide, under mild conditions, the theoretical
guarantees for the APHYNITY formulation to infer the correct MB/ML
decomposition, thus enabling both recovering the proper physical parameters
and accurate forecasting.

###  3.2 Solving APHYNITY with deep neural networks

In the following, both terms of the decomposition are parametrized and are
denoted as $F_{p}^{\theta_{p}}$ and $F_{p}^{\theta_{a}}$. Solving APHYNITY
then consists in estimating the parameters $\theta_{p}$ and $\theta_{a}$.
$\theta_{p}$ are the physical parameters and are typically low-dimensional,
e.g. 2 or 3 in our experiments for the considered physical models. For
$F_{a}$, we need sufficiently expressive models able to optimize over all
${\mathcal{F}}$: we thus use deep neural networks, which have shown promising
performances for the approximation of differential equations (Raissi et al.,
2019; Ayed et al., 2019b).

When learning the parameters of $F_{p}^{\theta_{p}}$ and $F_{a}^{\theta_{a}}$,
we have access to a finite dataset of trajectories discretized with a given
temporal resolution $\Delta t$:
${\mathcal{D}}_{\text{train}}=\\{(X^{(i)}_{k\Delta t})_{0\leq
k\leq\left\lfloor{\nicefrac{{T}}{{\Delta t}}}\right\rfloor}\\}_{1\leq i\leq
N}$. Solving Eq. (2) requires estimating the state derivative
$\nicefrac{{\mathop{}\\!\mathrm{d}X_{t}}}{{\mathop{}\\!\mathrm{d}t}}$
appearing in the constraint term. One solution is to approximate this
derivative using e.g. finite differences as in (Brunton et al., 2016;
Greydanus et al., 2019; Cranmer et al., 2020). This numerical scheme requires
high space and time resolutions in the observation space in order to get
reliable gradient estimates. Furthermore it is often unstable, leading to
explosive numerical errors as discussed in Appendix D. We propose instead to
solve Eq. (2) using an integral trajectory-based approach: we compute
$\widetilde{X}^{i}_{k\Delta t,X_{0}}$ from an initial state $X^{(i)}_{0}$
using the current $F_{p}^{\theta_{p}}+F_{a}^{\theta_{a}}$ dynamics, then
enforce the constraint $\widetilde{X}^{i}_{k\Delta t,X_{0}}=X^{i}_{k\Delta
t}$. This leads to our final objective function on $(\theta_{p},\theta_{a})$:

$\underset{\theta_{p},\theta_{a}}{\min}~{}~{}~{}\left\|F_{a}^{\theta_{a}}\right\|~{}~{}~{}\mathrm{subject~{}to}~{}~{}~{}~{}\forall
i,\forall k,\widetilde{X}^{(i)}_{k\Delta t}=X^{(i)}_{k\Delta t}$ (3)

where $\widetilde{X}^{(i)}_{k\Delta t}$ is the approximate solution of the
integral $\int_{X_{0}^{(i)}}^{X_{0}^{(i)}+k\Delta
t}(F_{p}^{\theta_{p}}+F_{a}^{\theta_{a}})(X_{s})\mathop{}\\!\mathrm{d}X_{s}$
obtained by a differentiable ODE solver.

In our setting, where we consider situations for which $F^{\theta_{p}}_{p}$
only partially describes the physical phenomenon, this coupled MB + ML
formulation leads to different parameter estimates than using the MB
formulation alone, as analyzed more thoroughly in Appendix C. Interestingly,
our experiments show that using this formulation also leads to a better
identification of the physical parameters $\theta_{p}$ than when fitting the
simplified physical model $F^{\theta_{p}}_{p}$ alone (Section 4). With only an
incomplete knowledge on the physics, $\theta_{p}$ estimator will be biased by
the additional dynamics which needs to be fitted in the data. Appendix F also
confirms that the integral formulation gives better forecasting results and a
more stable behavior than supervising over finite difference approximations of
the derivatives.

###  3.3 Adaptively constrained optimization

The formulation in Eq. (3) involves constraints which are difficult to enforce
exactly in practice. We considered a variant of the method of multipliers
(Bertsekas, 1996) which uses a sequence of Lagrangian relaxations
$\mathcal{L}_{\lambda_{j}}(\theta_{p},\theta_{a})$:

$\mathcal{L}_{\lambda_{j}}(\theta_{p},\theta_{a})=\|F_{a}^{\theta_{a}}\|+\lambda_{j}\cdot\mathcal{L}_{traj}(\theta_{p},\theta_{a})$
(4)

where
$\mathcal{L}_{traj}(\theta_{p},\theta_{a})=\sum_{i=1}^{N}\sum_{h=1}^{T/\Delta
t}\|X^{(i)}_{h\Delta t}-\widetilde{X}^{(i)}_{h\Delta t}\|$.

Initialization: $\lambda_{0}\geq 0,\tau_{1}>0,\tau_{2}>0$;

for  _epoch = $1:N_{epochs}$_ do

for  _iter in $1:N_{iter}$_ do

for  _batch in $1:B$_ do

$\theta_{j+1}=\theta_{j}-\tau_{1}\nabla\left[\lambda_{j}\mathcal{L}_{traj}(\theta_{j})+\left\|F_{a}\right\|\right]\;$

$\lambda_{j+1}=\lambda_{j}+$ $\tau_{2}\mathcal{L}_{traj}(\theta_{j+1})\;$

Algorithm 1 APHYNITY

This method needs an increasing sequence $(\lambda_{j})_{j}$ such that the
successive minima of $\mathcal{L}_{\lambda_{j}}$ converge to a solution (at
least a local one) of the constrained problem Eq. (3). We select
$(\lambda_{j})_{j}$ by using an iterative strategy: starting from a value
$\lambda_{0}$, we iterate, minimizing $\mathcal{L}_{\lambda_{j}}$ by gradient
descent222Convergence to a local minimum isn’t necessary, a few steps are
often sufficient for a successful optimization., then update $\lambda_{j}$
with: $\lambda_{j+1}=\lambda_{j}+\tau_{2}\mathcal{L}_{traj}(\theta_{j+1})$,
where $\tau_{2}$ is a chosen hyper-parameter and
$\theta=(\theta_{p},\theta_{a})$. This procedure is summarized in Algorithm 1.
This adaptive iterative procedure allows us to obtain stable and robust
results, in a reproducible fashion, as shown in the experiments.

##  4 Experimental validation

We validate our approach on 3 classes of challenging physical dynamics:
reaction-diffusion, wave propagation, and the damped pendulum, representative
of various application domains such as chemistry, biology or ecology (for
reaction-diffusion) and earth physic, acoustic, electromagnetism or even
neuro-biology (for waves equations). The two first dynamics are described by
PDEs and thus in practice should be learned from very high-dimensional
vectors, discretized from the original compact domain. This makes the learning
much more difficult than from the one-dimensional pendulum case. For each
problem, we investigate the cooperation between physical models of increasing
complexity encoding incomplete knowledge of the dynamics (denoted Incomplete
physics in the following) and data-driven models. We show the relevance of
APHYNITY (denoted APHYNITY models) both in terms of forecasting accuracy and
physical parameter identification.

###  4.1 Experimental setting

We describe the three families of equations studied in the experiments. In all
experiments, ${\mathcal{F}}=\mathcal{L}^{2}({\mathcal{A}})$ where
${\mathcal{A}}$ is the set of all admissible states for each problem, and the
$\mathcal{L}^{2}$ norm is computed on ${\mathcal{D}}_{train}$ by:
$\|F\|^{2}\approx\sum_{i,k}\|F(X^{(i)}_{k\Delta t})\|^{2}$. All considered
sets of physical functionals ${\mathcal{F}}_{p}$ are closed and convex in
${\mathcal{F}}$ and thus are Chebyshev. In order to enable the evaluation on
both prediction and parameter identification, all our experiments are
conducted on simulated datasets with known model parameters. Each dataset has
been simulated using an appropriate high-precision integration scheme for the
corresponding equation. All solver-based models take the first state $X_{0}$
as input and predict the remaining time-steps by integrating $F$ through the
same differentiable generic and common ODE solver (4th order Runge-
Kutta)333This integration scheme is then different from the one used for data
generation, the rationale for this choice being that when training a model one
does not know how exactly the data has been generated.. Implementation details
and architectures are given in Appendix E.

#### Reaction-diffusion equations

We consider a 2D FitzHugh-Nagumo type model (Klaasen & Troy, 1984). The system
is driven by the PDE $\frac{\partial u}{\partial t}=a\Delta
u+R_{u}(u,v;k),\frac{\partial v}{\partial t}=b\Delta v+R_{v}(u,v)$ where $a$
and $b$ are respectively the diffusion coefficients of $u$ and $v$, $\Delta$
is the Laplace operator. The local reaction terms are
$R_{u}(u,v;k)=u-u^{3}-k-v,R_{v}(u,v)=u-v$. The state is $X=(u,v)$ and is
defined over a compact rectangular domain $\Omega$ with periodic boundary
conditions. The considered physical models are:  • Param PDE ($a,b$), with
unknown ($a,b$) diffusion terms and without reaction terms:
${\mathcal{F}}_{p}=\\{F_{p}^{a,b}:(u,v)\mapsto(a\Delta u,b\Delta v)\ |\ a\geq
a_{\min}>0,b\geq b_{\min}>0\\}$;  • Param PDE ($a,b,k$), the full PDE with
unknown parameters:
${\mathcal{F}}_{p}=\\{F_{p}^{a,b,k}:(u,v)\mapsto\linebreak(a\Delta
u+R_{u}(u,v;k),b\Delta v+R_{v}(u,v)\ |\ a\geq a_{\min}>0,b\geq
b_{\min}>0,k\geq k_{\min}>0\\}$.

#### Damped wave equations

We investigate the damped-wave PDE: $\frac{\partial^{2}w}{\partial
t^{2}}-c^{2}\Delta w+k\frac{\partial w}{\partial t}=0$ where $k$ is the
damping coefficient. The state is $X=(w,\frac{\partial w}{\partial t})$ and we
consider a compact spatial domain $\Omega$ with Neumann homogeneous boundary
conditions. Note that this damping differs from the pendulum, as its effect is
global. Our physical models are:  • Param PDE ($c$), without damping term:
${\mathcal{F}}_{p}=\\{F^{c}_{p}:(u,v)\mapsto(v,c^{2}\Delta u)\ |\
c\in[\epsilon,+\infty)\text{ with }\epsilon>0\\}$;  • Param PDE ($c,k$):
${\mathcal{F}}_{p}=\\{F^{c,k}_{p}:(u,v)\mapsto(v,c^{2}\Delta u-kv)\ |\
c,k\in[\epsilon,+\infty)\text{ with }\epsilon>0\\}$.

#### Damped pendulum

The evolution follows the ODE
$\nicefrac{{\mathop{}\\!\mathrm{d}^{2}\theta}}{{\mathop{}\\!\mathrm{d}t^{2}}}+\omega_{0}^{2}\sin\theta+\alpha\nicefrac{{\mathop{}\\!\mathrm{d}\theta}}{{\mathop{}\\!\mathrm{d}t}}=0$,
where $\theta(t)$ is the angle, $\omega_{0}$ the proper pulsation ($T_{0}$ the
period) and $\alpha$ the damping coefficient. With state
$X=(\theta,\nicefrac{{\mathop{}\\!\mathrm{d}\theta}}{{\mathop{}\\!\mathrm{d}t}})$,
the ODE is
$F_{p}^{\omega_{0},\alpha}:X\mapsto(\nicefrac{{\mathop{}\\!\mathrm{d}\theta}}{{\mathop{}\\!\mathrm{d}t}},-\omega_{0}^{2}\sin\theta-\alpha\nicefrac{{\mathop{}\\!\mathrm{d}\theta}}{{\mathop{}\\!\mathrm{d}t}})$.
Our physical models are:  • Hamiltonian (Greydanus et al., 2019), a
conservative approximation, with
${\mathcal{F}}_{p}=\\{F^{\mathcal{H}}_{p}:(u,v)\mapsto(\partial_{y}\mathcal{H}(u,v),-\partial_{x}\mathcal{H}(u,v))\
|\ \mathcal{H}\in H^{1}(\mathbb{R}^{2})\\}$, $H^{1}(\mathbb{R}^{2})$ is the
first order Sobolev space.  • Param ODE ($\omega_{0}$), the frictionless
pendulum: ${\mathcal{F}}_{p}=\\{F^{\omega_{0},\alpha=0}_{p}|\
\omega_{0}\in[\epsilon,+\infty)\text{ with }\epsilon>0\\}$ • Param ODE
($\omega_{0},\alpha$), the full pendulum equation:
${\mathcal{F}}_{p}=\\{F^{\omega_{0},\alpha}_{p}|\
\omega_{0},\alpha\in[\epsilon,+\infty)\text{ with }\epsilon>0\\}$.

#### Baselines

As purely data-driven baselines, we use Neural ODE (Chen et al., 2018) for the
three problems and PredRNN++ (Wang et al., 2018, for reaction-diffusion only)
which are competitive models for datasets generated by differential equations
and for spatio-temporal data. As MB/ML methods, in the ablations studies (see
Appendix F), we compare for all problems, to the vanilla MB/ML cooperation
scheme found in (Wang et al., 2019; Mehta et al., 2020). We also show results
for True PDE/ODE, which corresponds to the equation for data simulation (which
do not lead to zero error due to the difference between simulation and
training integration schemes). For the pendulum, we compare to Hamiltonian
neural networks (Greydanus et al., 2019; Toth et al., 2020) and to the the
deep Galerkin method (DGM, Sirignano & Spiliopoulos, 2018). See additional
details in Appendix E.

###  4.2 Results

Table 1: Forecasting and identification results on the (a) reaction-diffusion,
(b) wave equation, and (c) damped pendulum datasets. We set for (a) $a=1\times
10^{-3},b=5\times 10^{-3},k=5\times 10^{-3}$, for (b) $c=330$, $k=50$ and for
(c) $T_{0}=6$, $\alpha=0.2$ as true parameters. $\log$ MSEs are computed
respectively over 25, 25, and 40 predicted time-steps. %Err param. averages
the results when several physical parameters are present. For each level of
incorporated physical knowledge, equivalent best results according to a
Student t-test are shown in bold. n/a corresponds to non-applicable cases.

Dataset |  | Method |  $\log$ MSE | %Err param. | $\|F_{a}\|^{2}$
---|---|---|---|---|---
(a) Reaction-diffusion |  Data-driven | Neural ODE | -3.76$\pm$0.02 | n/a | n/a
PredRNN++ | -4.60$\pm$0.01 | n/a | n/a
Incomplete physics | Param PDE ($a,b$) | -1.26$\pm$0.02 | 67.6 | n/a
APHYNITY Param PDE ($a,b$) | -5.10$\pm$0.21 | 2.3 | 67
Complete physics | Param PDE ($a,b,k$) | -9.34$\pm$0.20 | 0.17 | n/a
APHYNITY Param PDE ($a,b,k$) | -9.35$\pm$0.02 | 0.096 | 1.5e-6
3-6 | True PDE | -8.81$\pm$0.05 | n/a | n/a
| APHYNITY True PDE | -9.17$\pm$0.02 | n/a | 1.4e-7
(b) Wave equation | Data-driven | Neural ODE | -2.51$\pm$0.29 | n/a | n/a
Incomplete physics | Param PDE ($c$) | 0.51$\pm$0.07 | 10.4 | n/a
APHYNITY Param PDE ($c$) | -4.64$\pm$0.25 | 0.31 | 71.
Complete physics | Param PDE $(c,k)$ | -4.68$\pm$0.55 | 1.38 | n/a
APHYNITY Param PDE $(c,k)$ | -6.09$\pm$0.28 | 0.70 | 4.54
3-6 | True PDE | -4.66$\pm$0.30 | n/a | n/a
| APHYNITY True PDE | -5.24$\pm$0.45 | n/a | 0.14
(c) Damped pendulum | Data-driven | Neural ODE | -2.84$\pm$0.70 | n/a | n/a
Incomplete physics | Hamiltonian | -0.35$\pm$0.10 | n/a | n/a
APHYNITY Hamiltonian | -3.97$\pm$1.20 | n/a | 623
3-6 | Param ODE ($\omega_{0}$) | -0.14$\pm$0.10 | 13.2 | n/a
| Deep Galerkin Method ($\omega_{0}$) | -3.10$\pm$0.40 | 22.1 | n/a
| APHYNITY Param ODE ($\omega_{0}$) | -7.86$\pm$0.60 | 4.0 | 132
|  Complete physics | Param ODE ($\omega_{0},\alpha$) | -8.28$\pm$0.40 | 0.45 | n/a
| Deep Galerkin Method ($\omega_{0},\alpha$) | -3.14$\pm$0.40 | 7.1 | n/a
| APHYNITY Param ODE ($\omega_{0},\alpha$) | -8.31$\pm$0.30 | 0.39 | 8.5
3-6 | True ODE | -8.58$\pm$0.20 | n/a | n/a
| APHYNITY True ODE | -8.44$\pm$0.20 | n/a | 2.3

We analyze and discuss below the results obtained for the three kind of
dynamics. We successively examine different evaluation or quality criteria.
The conclusions are consistent for the three problems, which allows us to
highlight clear trends for all of them.

#### Forecasting accuracy

The data-driven models do not perform well compared to True PDE/ODE (all
values are test errors expressed as $\log$ MSE): -4.6 for PredRNN++ vs. -9.17
for reaction-diffusion, -2.51 vs. -5.24 for wave equation, and -2.84 vs. -8.44
for the pendulum in Table 4.2. The Deep Galerkin method for the pendulum in
complete physics DGM ($\omega_{0},\alpha$), being constrained by the equation,
outperforms Neural ODE but is far inferior to APHYNITY models. In the
incomplete physics case, DGM ($\omega_{0}$) fails to compensate for the
missing information. The incomplete physical models, Param PDE ($a,b$) for the
reaction-diffusion, Param PDE ($c$) for the wave equation, and Param ODE
($\omega_{0}$) and Hamiltonian models for the damped pendulum, have even
poorer performances than purely data-driven ones, as can be expected since
they ignore important dynamical components, e.g. friction in the pendulum
case. Using APHYNITY with these imperfect physical models greatly improves
forecasting accuracy in all cases, significantly outperforming purely data-
driven models, and reaching results often close to the accuracy of the true
ODE, when APHYNITY and the true ODE models are integrated with the same
numerical scheme (which is different from the one used for data generation,
hence the non-null errors even for the true equations), e.g. -5.92 vs. -5.24
for wave equation in Table 4.2. This clearly highlights the capacity of our
approach to augment incomplete physical models with a learned data-driven
component.

#### Physical parameter estimation

Confirming the phenomenon mentioned in the introduction and detailed in
Appendix C, incomplete physical models can lead to bad estimates for the
relevant physical parameters: an error respectively up to 67.6% and 10.4% for
parameters in the reaction-diffusion and wave equations, and an error of more
than 13% for parameters for the pendulum in Table 4.2. APHYNITY is able to
significantly improve physical parameters identification: 2.3% error for the
reaction-diffusion, 0.3% for the wave equation, and 4% for the pendulum. This
validates the fact that augmenting a simple physical model to compensate its
approximations is not only beneficial for prediction, but also helps to limit
errors for parameter identification when dynamical models do not fit data
well. This is crucial for interpretability and explainability of the
estimates.

#### Ablation study

We conduct ablation studies to validate the importance of the APHYNITY
augmentation compared to a naive strategy consisting in learning
$F=F_{p}+F_{a}$ without taking care on the quality of the decomposition, as
done in (Wang et al., 2019; Mehta et al., 2020). Results shown in Table 4.2 of
Appendix F show a consistent gain of APHYNITY for the three use cases and for
all physical models: for instance for Param ODE ($a,b$) in reaction-diffusion,
both forecasting performances ($\log\text{MSE}=$-5.10 vs. -4.56) and
identification parameter (Error$=$ 2.33% vs. 6.39%) improve. Other ablation
results are provided in Appendix F showing the relevance of the the
trajectory-based approach described in Section 3.2 (vs supervising over finite
difference approximations of the derivative $F$).

#### Flexibility

When applied to complete physical models, APHYNITY does not degrade accuracy,
contrary to a vanilla cooperation scheme (see ablations in Appendix F). This
is due to the least action principle of our approach: when the physical
knowledge is sufficient for properly predicting the observed dynamics, the
model learns to ignore the data-driven augmentation. This is shown by the norm
of the trained neural net component $F_{a}$, which is reported in Table 4.2
last column: as expected, $\|F_{a}\|^{2}$ diminishes as the complexity of the
corresponding physical model increases, and, relative to incomplete models,
the norm becomes very small for complete physical models (for example in the
pendulum experiments, we have $\|F_{a}\|=8.5$ for the APHYNITY model to be
compared with 132 and 623 for the incomplete models). Thus, we see that the
norm of $F_{a}$ is a good indication of how imperfect the physical models
${\mathcal{F}}_{p}$ are. It highlights the flexibility of APHYNITY to
successfully adapt to very different levels of prior knowledge. Note also that
APHYNITY sometimes slightly improves over the true ODE, as it compensates the
error introduced by different numerical integration methods for data
simulation and training (see Appendix E).

(a) Param PDE ($a$, $b$), diffusion-only

(b) APHYNITY Param PDE ($a$, $b$)

(c) Ground truth simulation

Figure 2: Comparison of predictions of two components $u$ (top) and $v$
(bottom) of the reaction-diffusion system. Note that $t=4$ is largely beyond
the dataset horizon ($t=2.5$).

(a) Neural ODE

(b) APHYNITY Param PDE ($c$)

(c) Ground truth simulation

Figure 3: Comparison between the prediction of APHYNITY when $c$ is estimated
and Neural ODE for the damped wave equation. Note that $t+32$, last column for
(a, b, c) is already beyond the training time horizon ($t+25$), showing the
consistency of APHYNITY method.

#### Qualitative visualizations

Results in Figure 2 for reaction-diffusion show that the incomplete diffusion
parametric PDE in Figure 2(a) is unable to properly match ground truth
simulations: the behavior of the two components in Figure 2(a) is reduced to
simple independent diffusions due to the lack of interaction terms between $u$
and $v$. By using APHYNITY in Figure 2(b), the correlation between the two
components appears together with the formation of Turing patterns, which is
very similar to the ground truth. This confirms that $F_{a}$ can learn the
reaction terms and improve prediction quality. In Figure 3, we see for the
wave equation that the data-driven Neural ODE model fails at approximating
$\nicefrac{{\mathop{}\\!\mathrm{d}w}}{{\mathop{}\\!\mathrm{d}t}}$ as the
forecast horizon increases: it misses crucial details for the second component
$\nicefrac{{\mathop{}\\!\mathrm{d}w}}{{\mathop{}\\!\mathrm{d}t}}$ which makes
the forecast diverge from the ground truth. APHYNITY incorporates a Laplacian
term as well as the data-driven $F_{a}$ thus capturing the damping phenomenon
and succeeding in maintaining physically sound results for long term
forecasts, unlike Neural ODE.

#### Extension to non-stationary dynamics

We provide additional results in Appendix G to tackle datasets where physical
parameters of the equations vary in each sequence. To this end, we design an
encoder able to perform parameter estimation for each sequence. Results show
that APHYNITY accommodates well to this setting, with similar trends as those
reported in this section.

#### Additional illustrations

We give further visual illustrations to demonstrate how the estimation of
parameters in incomplete physical models is improved with APHYNITY. For the
reaction-diffusion equation, we show that the incomplete parametric PDE
underestimates both diffusion coefficients. The difference is visually
recognizable between the poorly estimated diffusion (Figure 4(a)) and the true
one (Figure 4(c)) while APHYNITY gives a fairly good estimation of those
diffusion parameters as shown in Figure 4(b).

(a) $a=0.33\times 10^{-3},b=0.94\times 10^{-3}$, diffusion estimated with
Param PDE $(a,b)$

(b) $a=0.97\times 10^{-3},b=4.75\times 10^{-3}$, diffusion estimated with
APHYNITY Param PDE $(a,b)$

(c) $a=1.0\times 10^{-3},b=5.0\times 10^{-3}$, true diffusion

Figure 4: Diffusion predictions using coefficient learned with (a) incomplete
physical model Param PDE $(a,b)$ and (b) APHYNITY-augmented Param PDE$(a,b)$,
compared with the (c) true diffusion

##  5 Conclusion

In this work, we introduce the APHYNITY framework that can efficiently augment
approximate physical models with deep data-driven networks, performing
similarly to models for which the underlying dynamics are entirely known. We
exhibit the superiority of APHYNITY over data-driven, incomplete physics, and
state-of-the-art approaches combining ML and MB methods, both in terms of
forecasting and parameter identification on three various classes of physical
systems. Besides, APHYNITY is flexible enough to adapt to different
approximation levels of prior physical knowledge.

An appealing perspective is the applicability of APHYNITY on partially-
observable settings, such as video prediction. Besides, we hope that the
APHYNITY framework will open up the way to the design of a wide range of more
flexible MB/ML models, e.g. in climate science, robotics or reinforcement
learning. In particular, analyzing the theoretical decomposition properties in
a partially-observed setting is an important direction for future work.

## Acknowledgements:

Funding (P. Gallinari), Chaires de recherche et d’enseignement en intelligence
artificielle (Chaires IA), DL4Clim project.

## References

  * Ayed et al. (2019a) Ibrahim Ayed, Nicolas Cedilnik, Patrick Gallinari, and Maxime Sermesant.  Ep-net: Learning cardiac electrophysiology models for physiology-based constraints in data-driven predictions.  In Yves Coudière, Valéry Ozenne, Edward J. Vigmond, and Nejib Zemzemi (eds.), _Functional Imaging and Modeling of the Heart - 10th International Conference, FIMH 2019, Bordeaux, France, June 6-8, 2019, Proceedings_ , volume 11504 of _Lecture Notes in Computer Science_ , pp. 55–63. Springer, 2019a.
  * Ayed et al. (2019b) Ibrahim Ayed, Emmanuel de Bézenac, Arthur Pajot, Julien Brajard, and Patrick Gallinari.  Learning dynamical systems from partial observations.  _arXiv preprint arXiv:1902.11136_ , 2019b.
  * Becker et al. (2019) Philipp Becker, Harit Pandya, Gregor Gebhardt, Cheng Zhao, James Taylor, and Gerhard Neumann.  Recurrent kalman networks: Factorized inference in high-dimensional deep feature spaces.  _International Conference on Machine Learning (ICML)_ , 2019.
  * Bertsekas (1996) Dimitri P. Bertsekas.  _Constrained Optimization and Lagrange Multiplier Methods (Optimization and Neural Computation Series)_.  Athena Scientific, 1 edition, 1996.
  * Brunton et al. (2016) Steven L. Brunton, Joshua L. Proctor, and J. Nathan Kutz.  Discovering governing equations from data by sparse identification of nonlinear dynamical systems.  _Proceedings of the National Academy of Sciences_ , 113(15):3932–3937, 2016.
  * Chen et al. (2018) Tian Qi Chen, Yulia Rubanova, Jesse Bettencourt, and David K. Duvenaud.  Neural ordinary differential equations.  In _Advances in neural information processing systems (NeurIPS)_ , pp. 6571–6583, 2018.
  * Chen (2004) Wen-Hua Chen.  Disturbance observer based control for nonlinear systems.  _IEEE/ASME transactions on mechatronics_ , 9(4):706–710, 2004.
  * Chen et al. (2020) Zhengdao Chen, Jianyu Zhang, Martin Arjovsky, and Léon Bottou.  Symplectic recurrent neural networks.  _International Conference on Learning Representations (ICLR)_ , 2020\.
  * Choi et al. (2016) Edward Choi, Mohammad Taha Bahadori, Jimeng Sun, Joshua Kulas, Andy Schuetz, and Walter Stewart.  RETAIN: An interpretable predictive model for healthcare using reverse time attention mechanism.  In _Advances in Neural Information Processing Systems (NeurIPS)_ , pp. 3504–3512, 2016.
  * Courtier et al. (1994) Philippe Courtier, J-N Thépaut, and Anthony Hollingsworth.  A strategy for operational implementation of 4d-var, using an incremental approach.  _Quarterly Journal of the Royal Meteorological Society_ , 120(519):1367–1387, 1994.
  * Cranmer et al. (2020) Miles Cranmer, Sam Greydanus, Stephan Hoyer, Peter Battaglia, David Spergel, and Shirley Ho.  Lagrangian neural networks.  _ICLR 2020 Deep Differential Equations Workshop_ , 2020.
  * de Bézenac et al. (2018) Emmanuel de Bézenac, Arthur Pajot, and Patrick Gallinari.  Deep learning for physical processes: Incorporating prior scientific knowledge.  _International Conference on Learning Representations (ICLR)_ , 2018\.
  * Donà et al. (2020) Jérémie Donà, Jean-Yves Franceschi, Sylvain Lamprier, and Patrick Gallinari.  Pde-driven spatiotemporal disentanglement.  _International Conference on Learning Representations (ICLR)_ , 2020\.
  * Dormand & Prince (1980) John R Dormand and Peter J Prince.  A family of embedded runge-kutta formulae.  _Journal of computational and applied mathematics_ , 6(1):19–26, 1980.
  * Fletcher & Moors (2014) James Fletcher and Warren Moors.  Chebyshev sets.  _Journal of the Australian Mathematical Society_ , 98:161–231, 04 2014.  doi: 10.1017/S1446788714000561.
  * Gentine et al. (2018) P. Gentine, M. Pritchard, S. Rasp, G. Reinaudi, and G. Yacalis.  Could machine learning break the convection parameterization deadlock?  _Geophysical Research Letters_ , 45(11):5742–5751, 2018.
  * Greydanus et al. (2019) Samuel Greydanus, Misko Dzamba, and Jason Yosinski.  Hamiltonian neural networks.  In _Advances in Neural Information Processing Systems (NeurIPS)_ , pp. 15353–15363, 2019.
  * Janner et al. (2019) Michael Janner, Justin Fu, Marvin Zhang, and Sergey Levine.  When to trust your model: Model-based policy optimization.  In _Advances in Neural Information Processing Systems (NeurIPS)_ , pp. 12519–12530, 2019.
  * Johnson (1987) Gordon G Johnson.  A nonconvex set which has the unique nearest point property.  _Journal of Approximation Theory_ , 51(4):289 – 332, 1987.
  * Kalman (1960) Rudolph Emil Kalman.  A new approach to linear filtering and prediction problems.  1960\.
  * Klaasen & Troy (1984) Gene A. Klaasen and William C. Troy.  Stationary wave solutions of a system of reaction-diffusion equations derived from the fitzhugh–nagumo equations.  _SIAM Journal on Applied Mathematics_ , 44(1):96–110, 1984.  doi: 10.1137/0144008.
  * Large & Yeager (2004) William Large and Stephen Yeager.  Diurnal to decadal global forcing for ocean and sea-ice models: The data sets and flux climatologies, 05 2004.
  * Le Guen & Thome (2020) Vincent Le Guen and Nicolas Thome.  Disentangling physical dynamics from unknown factors for unsupervised video prediction.  In _Computer Vision and Pattern Recognition (CVPR)_. 2020.
  * Li et al. (2014) Shihua Li, Jun Yang, Wen-Hua Chen, and Xisong Chen.  _Disturbance observer-based control: methods and applications_.  CRC press, 2014.
  * Long et al. (2018a) Yun Long, Xueyuan She, and Saibal Mukhopadhyay.  Hybridnet: integrating model-based and data-driven learning to predict evolution of dynamical systems.  _Conference on Robot Learning (CoRL)_ , 2018a.
  * Long et al. (2018b) Zichao Long, Yiping Lu, Xianzhong Ma, and Bin Dong.  PDE-Net: Learning PDEs from data.  In _International Conference on Machine Learning (ICML)_ , 2018b.
  * Mehta et al. (2020) Viraj Mehta, Ian Char, Willie Neiswanger, Youngseog Chung, and Jeff Schneider.  Neural dynamical systems.  _ICLR 2020 Deep Differential Equations Workshop_ , 2020.
  * Nagabandi et al. (2018) Anusha Nagabandi, Gregory Kahn, Ronald S Fearing, and Sergey Levine.  Neural network dynamics for model-based deep reinforcement learning with model-free fine-tuning.  In _2018 IEEE International Conference on Robotics and Automation (ICRA)_ , pp. 7559–7566. IEEE, 2018.
  * Oreshkin et al. (2020) Boris N. Oreshkin, Dmitri Carpov, Nicolas Chapados, and Yoshua Bengio.  N-BEATS: Neural basis expansion analysis for interpretable time series forecasting.  _International Conference on Learning Representations (ICLR)_ , 2020\.
  * Paszke et al. (2019) Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, Alban Desmaison, Andreas Kopf, Edward Yang, Zachary DeVito, Martin Raison, Alykhan Tejani, Sasank Chilamkurthy, Benoit Steiner, Lu Fang, Junjie Bai, and Soumith Chintala.  Pytorch: An imperative style, high-performance deep learning library.  In H. Wallach, H. Larochelle, A. Beygelzimer, F. d’Alché Buc, E. Fox, and R. Garnett (eds.), _Advances in Neural Information Processing Systems 32_ , pp. 8024–8035. 2019.
  * Pernot & Cailliez (2017) Pascal Pernot and Fabien Cailliez.  A critical review of statistical calibration/prediction models handling data inconsistency and model inadequacy.  _AIChE Journal_ , 63(10):4642–4665, 2017.
  * Psichogios & Ungar (1992) Dimitris C Psichogios and Lyle H Ungar.  A hybrid neural network-first principles approach to process modeling.  _AIChE Journal_ , 38(10):1499–1511, 1992.
  * Raissi et al. (2019) Maziar Raissi, Paris Perdikaris, and George Em Karniadakis.  Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations.  _Journal of Computational Physics_ , 473:686–707, 2019\.
  * Reichstein et al. (2019) Markus Reichstein, Gustau Camps-Valls, Bjorn Stevens, Martin Jung, Joachim Denzler, Nuno Carvalhais, and & Prabhat.  Deep learning and process understanding for data-driven Earth system science.  _Nature_ , 566:195–204, 2019.
  * Rico-Martinez et al. (1994) R Rico-Martinez, JS Anderson, and IG Kevrekidis.  Continuous-time nonlinear signal processing: a neural network based approach for gray box identification.  In _Proceedings of IEEE Workshop on Neural Networks for Signal Processing_ , pp. 596–605. IEEE, 1994.
  * Rolnick et al. (2019) David Rolnick, Priya L Donti, Lynn H Kaack, Kelly Kochanski, Alexandre Lacoste, Kris Sankaran, Andrew Slavin Ross, Nikola Milojevic-Dupont, Natasha Jaques, Anna Waldman-Brown, et al.  Tackling climate change with machine learning.  In _NeurIPS 2019 workshop on Climate Change with Machine Learning_ , 2019.
  * Saha et al. (2020) Priyabrata Saha, Saurabh Dash, and Saibal Mukhopadhyay.  PhICNet: Physics-incorporated convolutional recurrent neural networks for modeling dynamical systems.  _arXiv preprint arXiv:2004.06243_ , 2020.
  * Seo et al. (2020) Sungyong Seo, Chuizheng Meng, and Yan Liu.  Physics-aware difference graph networks for sparsely-observed dynamics.  _International Conference on Learning Representations (ICLR)_ , 2020\.
  * Shi et al. (2015) Xingjian Shi, Zhourong Chen, Hao Wang, Dit-Yan Yeung, Wai-Kin Wong, and Wang-chun Woo.  Convolutional LSTM network: A machine learning approach for precipitation nowcasting.  In _Advances in neural information processing systems (NeurIPS)_ , pp. 802–810, 2015.
  * Sirignano & Spiliopoulos (2018) Justin Sirignano and Konstantinos Spiliopoulos.  Dgm: A deep learning algorithm for solving partial differential equations.  _Journal of computational physics_ , 375:1339–1364, 2018\.
  * Thompson & Kramer (1994) Michael L Thompson and Mark A Kramer.  Modeling chemical processes using prior knowledge and neural networks.  _AIChE Journal_ , 40(8):1328–1340, 1994.
  * Toth et al. (2020) Peter Toth, Danilo Jimenez Rezende, Andrew Jaegle, Sébastien Racanière, Aleksandar Botev, and Irina Higgins.  Hamiltonian generative networks.  _International Conference on Learning Representations (ICLR)_ , 2020\.
  * Toubeau et al. (2018) Jean-François Toubeau, Jérémie Bottieau, François Vallée, and Zacharie De Grève.  Deep learning-based multivariate probabilistic forecasting for short-term scheduling in power markets.  _IEEE Transactions on Power Systems_ , 34(2):1203–1215, 2018.
  * Ummenhofer et al. (2020) Benjamin Ummenhofer, Lukas Prantl, Nils Thuerey, and Vladlen Koltun.  Lagrangian fluid simulation with continuous convolutions.  _International Conference on Learning Representations (ICLR)_ , 2020\.
  * Wang et al. (2019) Qi Wang, Feng Li, Yi Tang, and Yan Xu.  Integrating model-driven and data-driven methods for power system frequency stability assessment and control.  _IEEE Transactions on Power Systems_ , 34(6):4557–4568, 2019.
  * Wang et al. (2018) Yunbo Wang, Zhifeng Gao, Mingsheng Long, Jianmin Wang, and Philip S. Yu.  PredRNN++: Towards a resolution of the deep-in-time dilemma in spatiotemporal predictive learning.  In _International Conference on Machine Learning (ICML)_ , 2018.

##  Appendix A Reminder on proximinal and Chebyshev sets

We begin by giving a definition of proximinal and Chebyshev sets, taken from
(Fletcher & Moors, 2014):

######  Definition 1.

A proximinal set of a normed space $(E,\|\cdot\|)$ is a subset
$\mathcal{C}\subset E$ such that every $x\in E$ admits at least a nearest
point in $\mathcal{C}$.

######  Definition 2.

A Chebyshev set of a normed space $(E,\|\cdot\|)$ is a subset
$\mathcal{C}\subset E$ such that every $x\in E$ admits a unique nearest point
in $\mathcal{C}$.

Proximinality reduces to a compacity condition in finite dimensional spaces.
In general, it is a weaker one: Boundedly compact sets verify this property
for example.

In Euclidean spaces, Chebyshev sets are simply the closed convex subsets. The
question of knowing whether it is the case that all Chebyshev sets are closed
convex sets in infinite dimensional Hilbert spaces is still an open question.
In general, there exists examples of non-convex Chebyshev sets, a famous one
being presented in (Johnson, 1987) for a non-complete inner-product space.

Given the importance of this topic in approximation theory, finding necessary
conditions for a set to be Chebyshev and studying the properties of those sets
have been the subject of many efforts. Some of those properties are summarized
below:

  * •

The metric projection on a boundedly compact Chebyshev set is continuous.

  * •

If the norm is strict, every closed convex space, in particular any finite
dimensional subspace is Chebyshev.

  * •

In a Hilbert space, every closed convex set is Chebyshev.

##  Appendix B Proof of Propositions 1 and 2

We prove the following result which implies both propositions in the article:

######  Proposition 3.

The optimization problem:

$\underset{F_{p}\in{\mathcal{F}}_{p},F_{a}\in{\mathcal{F}}}{\min}~{}~{}~{}\left\|F_{a}\right\|~{}~{}~{}\mathrm{subject~{}to}~{}~{}~{}~{}\forall
X\in{\mathcal{D}},\forall
t,\frac{\mathop{}\\!\mathrm{d}X_{t}}{\mathop{}\\!\mathrm{d}t}=(F_{p}+F_{a})(X_{t})$
(5)

is equivalent a metric projection onto ${\mathcal{F}}_{p}$.

If ${\mathcal{F}}_{p}$ is proximinal, Eq. (5) admits a minimizing pair.

If ${\mathcal{F}}_{p}$ is Chebyshev, Eq. (5) admits a unique minimizing pair
which $F_{p}$ is the metric projection.

###### Proof.

The idea is to reconstruct the full functional from the trajectories of
${\mathcal{D}}$. By definition, ${\mathcal{A}}$ is the set of points reached
by trajectories in ${\mathcal{D}}$ so that:

${\mathcal{A}}=\\{x\in\mathbb{R}^{d}\ |\ \exists
X_{\cdot}\in{\mathcal{D}},\exists t,\ X_{t}=x\\}$

Then let us define a function $F^{\mathcal{D}}$ in the following way: For
$a\in{\mathcal{A}}$, we can find $X_{\cdot}\in{\mathcal{D}}$ and $t_{0}$ such
that $X_{t_{0}}=a$. Differentiating $X$ at $t_{0}$, which is possible by
definition of ${\mathcal{D}}$, we take:

$F^{\mathcal{D}}(a)=\left.\frac{\mathop{}\\!\mathrm{d}X_{t}}{\mathop{}\\!\mathrm{d}t}\right|_{t=t_{0}}$

For any $(F_{p},F_{a})$ satisfying the constraint in Eq. (5), we then have
that
$(F_{p}+F_{a})(a)=\nicefrac{{\mathop{}\\!\mathrm{d}X_{t}}}{{\mathop{}\\!\mathrm{d}t}}_{|t_{0}}=F^{\mathcal{D}}(a)$
for all $a\in{\mathcal{A}}$. Conversely, any pair such that
$(F_{p},F_{a})\in{\mathcal{F}}_{p}\times{\mathcal{F}}$ and
$F_{p}+F_{a}=F^{\mathcal{D}}$, verifies the constraint.

Thus we have the equivalence between Eq. (5) and the metric projection
formulated as: F_p∈F_p∥F^D\- F_p ∥

If ${\mathcal{F}}_{p}$ is proximinal, the projection problem admits a solution
which we denote $F^{\star}_{p}$. Taking
$F^{\star}_{a}=F^{\mathcal{D}}-F^{\star}_{p}$, we have that
$F^{\star}_{p}+F^{\star}_{a}=F^{\mathcal{D}}$ so that
$(F^{\star}_{p},F^{\star}_{a})$ verifies the constraint of Eq. (2). Moreover,
if there is $(F_{p},F_{a})$ satisfying the constraint of Eq. (2), we have that
$F_{p}+F_{a}=F^{\mathcal{D}}$ by what was shown above and
$\|F_{a}\|=\|F^{\mathcal{D}}-F_{p}\|\geq\|F^{\mathcal{D}}-F^{\star}_{p}\|$ by
definition of $F^{\star}_{p}$. This shows that $(F^{\star}_{p},F^{\star}_{a})$
is minimal.

Moreover, if ${\mathcal{F}}_{p}$ is a Chebyshev set, by uniqueness of the
projection, if $F_{p}\not=F^{\star}_{p}$ then $\|F_{a}\|>\|F^{\star}_{a}\|$.
Thus the minimal pair is unique.

∎

##  Appendix C Parameter estimation in incomplete physical models

Classically, when a set ${\mathcal{F}}_{p}\subset{\mathcal{F}}$ summarising
the most important properties of a system is available, this gives a
simplified model of the true dynamics and the adopted problem is then to fit
the trajectories using this model as well as possible, solving: F_p∈F_pE_X∼D
L(~X^X_0,X) ∀g∈I, ~X_0^g = g and ∀t, d~Xtgdt = F_p(~X_t^g) where $L$ is a
discrepancy measure between trajectories. Recall that $\widetilde{X}^{X_{0}}$
is the result trajectory of an ODE solver taking $X_{0}$ as initial condition.
In other words, we try to find a function $F_{p}$ which gives trajectories as
close as possible to the ones from the dataset. While estimation of the
function becomes easier, there is then a residual part which is left
unexplained and this can be a non negligible issue in at least two ways:

  * •

When $F\not\in{\mathcal{F}}_{p}$, the loss is strictly positive at the
minimum. This means that reducing the space of functions $\mathcal{F}_{p}$
makes us lose in terms of accuracy.444This is true in theory, although not
necessarily in practice when $F$ overfits a small dataset.

  * •

The obtained function $F_{p}$ might not even be the most meaningful function
from ${\mathcal{F}}_{p}$ as it would try to capture phenomena which are not
explainable with functions in ${\mathcal{F}}_{p}$, thus giving the wrong bias
to the calculated function. For example, if one is considering a dampened
periodic trajectory where only the period can be learned in
${\mathcal{F}}_{p}$ but not the dampening, the estimated period will account
for the dampening and will thus be biased.

This is confirmed in the paper in Section 4: the incomplete physical models
augmented with APHYNITY get different and experimentally better physical
identification results than the physical models alone.

Let us compare our approach with this one on the linearized damped pendulum to
show how estimates of physical parameters can differ. The equation is the
following:

$\frac{\mathop{}\\!\mathrm{d}^{2}\theta}{\mathop{}\\!\mathrm{d}t^{2}}+\omega_{0}^{2}\theta+\alpha\frac{\mathop{}\\!\mathrm{d}\theta}{\mathop{}\\!\mathrm{d}t}=0$

We take the same notations as in the article and parametrize the simplified
physical models as:

$F^{a}_{p}:X\mapsto(\frac{\mathop{}\\!\mathrm{d}\theta}{\mathop{}\\!\mathrm{d}t},-a\theta)$

where $a>0$ corresponds to $\omega_{0}^{2}$. The corresponding solution for an
initial state $X_{0}$, which we denote $X^{a}$, can then written explicitly
as:

$\theta^{a}_{t}=\theta_{0}\cos{\sqrt{a}t}$

Let us consider damped pendulum solutions $X$ written as:

$\theta_{t}=\theta_{0}e^{-t}\cos{t}$

which corresponds to:

$F:X\mapsto(\frac{\mathop{}\\!\mathrm{d}\theta}{\mathop{}\\!\mathrm{d}t},-2(\theta+\frac{\mathop{}\\!\mathrm{d}\theta}{\mathop{}\\!\mathrm{d}t}))$

It is then easy to see that the estimate of $a$ with the physical model alone
can be obtained by minimizing:

$\int_{0}^{T}|e^{-t}\cos{t}-\cos{\sqrt{a}t}|^{2}$

This expression depends on $T$ and thus, depending on the chosen time interval
and the way the integral is discretized will almost always give biased
estimates. In other words, the estimated value of $a$ will not give us the
desired solution $t\mapsto\cos{t}$.

On the other hand, for a given $a$, in the APHYNITY framework, the residual
must be equal to:

$F^{a}_{r}:X\mapsto(0,(a-2)\theta-2\frac{\mathop{}\\!\mathrm{d}\theta}{\mathop{}\\!\mathrm{d}t})$

in order to satisfy the fitting constraint. Here $a$ corresponds to
$1+\omega_{0}^{2}$ not to $\omega_{0}^{2}$ as in the simplified case.
Minimizing its norm, we obtain $a=2$ which gives us the desired solution:

$\theta_{t}=\theta_{0}e^{-t}\cos{t}$

with the right period.

##  Appendix D Discussion on supervision over derivatives

In order to find the appropriate decomposition $(F_{p},F_{a})$, we use a
trajectory-based error by solving: F_p∈F_p, F_a∈F∥F_a ∥ ∀g∈I, ~X_0^g = g and
∀t, d~Xtgdt = (F_p+F_a)(~X_t^g) ∀X∈D, L(X,~X^X_0) = 0

In the continuous setting where the data is available at all times $t$, this
problem is in fact equivalent to the following one: F_p∈F_pE_X∼D ∫∥dXtdt -
F_p(X_t) ∥ where the supervision is done directly over derivatives, obtained
through finite-difference schemes. This echoes the proof in Section B of the
Appendix where $F$ can be reconstructed from the continuous data.

However, in practice, data is only available at discrete times with a certain
time resolution. While Eq. (D) is indeed equivalent to Eq. (D) in the
continuous setting, in the practical discrete one, the way error propagates is
not anymore: For Eq. (D) it is controlled over integrated trajectories while
for Eq. (D) the supervision is over the approximate derivatives of the
trajectories from the dataset. We argue that the trajectory-based approach is
more flexible and more robust for the following reasons:

  * •

In Eq. (D), if $F_{a}$ is appropriately parameterized, it is possible to
perfectly fit the data trajectories at the sampled points.

  * •

The use of finite differences schemes to estimate $F$ as is done in Eq. (D)
necessarily induces a non-zero discretization error.

  * •

This discretization error is explosive in terms of divergence from the true
trajectories.

This last point is quite important, especially when time sampling is sparse
(even though we do observe this adverse effect empirically in our experiments
with relatively finely time-sampled trajectories). The following gives a
heuristical reasoning as to why this is the case. Let
$\widetilde{F}=F+\epsilon$ be the function estimated from the sampled points
with an error $\epsilon$ such that $\|\epsilon\|_{\infty}\leq\alpha$. Denoting
$\widetilde{X}$ the corresponding trajectory generated by $\widetilde{F}$, we
then have, for all $X\in{\mathcal{D}}$:

$\forall t,\
\frac{\mathop{}\\!\mathrm{d}(X-\widetilde{X})_{t}}{\mathop{}\\!\mathrm{d}t}=F(X_{t})-F(\widetilde{X}_{t})-\epsilon(\widetilde{X}_{t})$

Integrating over $[0,T]$ and using the triangular inequality as well as the
mean value inequality, supposing that $F$ has uniformly bounded spatial
derivatives:

$\forall t\in[0,T],\ \|(X-\widetilde{X})_{t}\|\leq\|\nabla
F\|_{\infty}\int_{0}^{t}\|X_{s}-\widetilde{X}_{s}\|+\alpha t$

which, using a variant of the Grönwall lemma, gives us the inequality:

$\forall t\in[0,T],\ \|X_{t}-\widetilde{X}_{t}\|\leq\frac{\alpha}{\|\nabla
F\|_{\infty}}(\exp(\|\nabla F\|_{\infty}t)-1)$

When $\alpha$ tends to $0$, we recover the true trajectories $X$. However, as
$\alpha$ is bounded away from $0$ by the available temporal resolution, this
inequality gives a rough estimate of the way $\widetilde{X}$ diverges from
them, and it can be an equality in many cases. This exponential behaviour
explains our choice of a trajectory-based optimization.

##  Appendix E Implementation details

We describe here the three use cases studied in the paper for validating
APHYNITY. All experiments are implemented with PyTorch (Paszke et al., 2019)
and the differentiable ODE solvers with the adjoint method implemented in
torchdiffeq.555https://github.com/rtqichen/torchdiffeq

###  E.1 Reaction-diffusion equations

The system is driven by a FitzHugh-Nagumo type PDE (Klaasen & Troy, 1984)

$\frac{\partial u}{\partial t}=a\Delta u+R_{u}(u,v;k),\frac{\partial
v}{\partial t}=b\Delta v+R_{v}(u,v)$

where $a$ and $b$ are respectively the diffusion coefficients of $u$ and $v$,
$\Delta$ is the Laplace operator. The local reaction terms are
$R_{u}(u,v;k)=u-u^{3}-k-v,R_{v}(u,v)=u-v$.

The state $X=(u,v)$ is defined over a compact rectangular domain
$\Omega=[-1,1]^{2}$ with periodic boundary conditions. $\Omega$ is spatially
discretized with a $32\times 32$ 2D uniform square mesh grid. The periodic
boundary condition is implemented with circular padding around the borders.
$\Delta$ is systematically estimated with a $3\times 3$ discrete Laplace
operator.

#### Dataset

Starting from a randomly sampled initial state
$X_{\text{init}}\in[0,1]^{2\times 32\times 32}$, we generate states by
integrating the true PDE with fixed $a$, $b$, and $k$ in a dataset ($a=1\times
10^{-3},b=5\times 10^{-3},k=5\times 10^{-3}$). We firstly simulate high time-
resolution ($\delta t_{\text{sim}}=0.001$) sequences with explicit finite
difference method. We then extract states every $\delta t_{\text{data}}=0.1$
to construct our low time-resolution datasets.

We set the time of random initial state to $t=-0.5$ and the time horizon to
$t=2.5$. 1920 sequences are generated, with 1600 for training/validation and
320 for test. We take the state at $t=0$ as $X_{0}$ and predict the sequence
until the horizon (equivalent to 25 time steps) in all reaction-diffusion
experiments. Note that the sub-sequence with $t<0$ are reserved for the
extensive experiments in Appendix G.1.

#### Neural network architectures

Our $F_{a}$ here is a 3-layer convolution network (ConvNet). The two input
channels are $(u,v)$ and two output ones are $(\frac{\partial u}{\partial
t},\frac{\partial v}{\partial t})$. The purely data-driven Neural ODE uses
such ConvNet as its $F$. The detailed architecture is provided in Table 2. The
estimated physical parameters $\theta_{p}$ in $F_{p}$ are simply a trainable
vector $(a,b)\in\mathbb{R}_{+}^{2}$ or $(a,b,k)\in\mathbb{R}_{+}^{3}$.

Table 2: ConvNet architecture in reaction-diffusion and wave equation experiments, used as data-driven derivative operator in APHYNITY and Neural ODE (Chen et al., 2018). Module | Specification
---|---
2D Conv. |  $3\times 3$ kernel, 2 input channels, 16 output channels, 1 pixel zero padding
2D Batch Norm. | No average tracking
ReLU activation | —
2D Conv. |  $3\times 3$ kernel, 16 input channels, 16 output channels, 1 pixel zero padding
2D Batch Norm. | No average tracking
ReLU activation | —
2D Conv. |  $3\times 3$ kernel, 16 input channels, 2 output channels, 1 pixel zero padding

#### Optimization hyperparameters

We choose to apply the same hyperparameters for all the reaction-diffusion
experiments: $Niter=1,\lambda_{0}=1,\tau_{1}=1\times 10^{-3},\tau_{2}=1\times
10^{3}$.

###  E.2 Wave equations

The damped wave equation is defined by

$\frac{\partial^{2}w}{\partial t^{2}}-c^{2}\Delta w+k\frac{\partial
w}{\partial t}=0$

where $c$ is the wave speed and $k$ is the damping coefficient. The state is
$X=(w,\frac{\partial w}{\partial t})$.

We consider a compact spatial domain $\Omega$ represented as a $64\times 64$
grid and discretize the Laplacian operator similarly. $\Delta$ is implemented
using a $5\times 5$ discrete Laplace operator in simulation whereas in the
experiment is a $3\times 3$ Laplace operator. Null Neumann boundary condition
are imposed for generation.

#### Dataset

$\delta t$ was set to $0.001$ to respect Courant number and provide stable
integration. The simulation was integrated using a 4th order finite difference
Runge-Kutta scheme for 300 steps from an initial Gaussian state, i.e for all
sequence at $t=0$, we have:

$w(x,y,t=0)=C\times\exp^{\frac{(x-x_{0})^{2}+(y-y_{0})^{2}}{\sigma^{2}}}$ (6)

The amplitude $C$ is fixed to $1$, and $(x_{0},y_{0})=(32,32)$ to make the
Gaussian curve centered for all sequences. However, $\sigma$ is different for
each sequence and uniformly sampled in $[10,100]$. The same $\delta t$ was
used for train and test. All initial conditions are Gaussian with varying
amplitudes. 250 sequences are generated, 200 are used for training while 50
are reserved as a test set. In the main paper setting, $c=330$ and $k=50$. As
with the reaction diffusion case, the algorithm takes as input a state
$X_{t_{0}}=(w,\frac{\mathop{}\\!\mathrm{d}w}{\mathop{}\\!\mathrm{d}t})(t_{0})$
and predicts all states from $t_{0}+\delta t$ up to $t_{0}+25\delta t$.

#### Neural network architectures

The neural network for $F_{a}$ is a 3-layer convolution neural network with
the same architecture as in Table 2. For $F_{p}$, the parameter(s) to be
estimated is either a scalar $c\in\mathbb{R}_{+}$ or a vector
$(c,k)\in\mathbb{R}_{+}^{2}$. Similarly, Neural ODE networks are build as
presented in Table 2.

#### Optimization hyperparameters

We use the same hyperparameters for the experiments:
${Niter=3,\lambda_{0}=1,\tau_{1}=1\times 10^{-4},\tau_{2}=1\times 10^{2}}$.

###  E.3 Damped pendulum

We consider the non-linear damped pendulum problem, governed by the ODE

$\frac{\mathop{}\\!\mathrm{d}^{2}\theta}{\mathop{}\\!\mathrm{d}t^{2}}+\omega_{0}^{2}\sin\theta+\alpha\frac{\mathop{}\\!\mathrm{d}\theta}{\mathop{}\\!\mathrm{d}t}=0$

where $\theta(t)$ is the angle, $\omega_{0}=\frac{2\pi}{T_{0}}$ is the proper
pulsation ($T_{0}$ being the period) and $\alpha$ is the damping coefficient.
With the state
$X=(\theta,\frac{\mathop{}\\!\mathrm{d}\theta}{\mathop{}\\!\mathrm{d}t})$, the
ODE can be written as
$\frac{\mathop{}\\!\mathrm{d}X_{t}}{\mathop{}\\!\mathrm{d}t}=F(X_{t})$ with
$F:X\mapsto(\frac{\mathop{}\\!\mathrm{d}\theta}{\mathop{}\\!\mathrm{d}t},-\omega_{0}^{2}\sin\theta-\alpha\frac{\mathop{}\\!\mathrm{d}\theta}{\mathop{}\\!\mathrm{d}t})$.

#### Dataset

For each train / validation / test split, we simulate a dataset with 25
trajectories of 40 timesteps (time interval $[0,20]$, timestep $\delta t=0.5$)
with fixed ODE coefficients $(T_{0}=12,\alpha=0.2)$ and varying initial
conditions. The simulation integrator is Dormand-Prince Runge-Kutta method of
order (4)5 (DOPRI5, Dormand & Prince, 1980). We also add a small amount of
white gaussian noise ($\sigma=0.01$) to the state. Note that our pendulum
dataset is much more challenging than the ideal frictionless pendulum
considered in Greydanus et al. (2019).

#### Neural network architectures

We detail in Table 3 the neural architectures used for the damped pendulum
experiments. All data-driven augmentations for approximating the mapping
$X_{t}\mapsto F(X_{t})$ are implemented by multi-layer perceptrons (MLP) with
3 layers of 200 neurons and ReLU activation functions (except at the last
layer: linear activation). The Hamiltonian (Greydanus et al., 2019; Toth et
al., 2020) is implemented by a MLP that takes the state $X_{t}$ and outputs a
scalar estimation of the Hamiltonian $\mathcal{H}$ of the system: the
derivative is then computed by an in-graph gradient of $\mathcal{H}$ with
respect to the input:
$F(X_{t})=\left(\frac{\partial\mathcal{H}}{\partial(\mathop{}\\!\mathrm{d}\theta/\mathop{}\\!\mathrm{d}t)},-\frac{\partial\mathcal{H}}{\mathop{}\\!\mathrm{d}\theta}\right)$.

Table 3: Neural network architectures for the damped pendulum experiments. n/a
corresponds to non-applicable cases.

Method Physical model Data-driven model Neural ODE n/a MLP(in=2, units=200,
layers=3, out=2) Hamiltonian MLP(in=2, units=200, layers=3, out=1) n/a
APHYNITY Hamiltonian MLP(in=2, units=200, layers=3, out=1) MLP(in=2,
units=200, layers=3, out=2) Param ODE ($\omega_{0}$) 1 trainable parameter
$\omega_{0}$ n/a APHYNITY Param ODE ($\omega_{0}$) 1 trainable parameter
$\omega_{0}$ MLP(in=2, units=200, layers=3, out=2) Param ODE
($\omega_{0},\alpha$) 2 trainable parameters $\omega_{0},\lambda$ n/a APHYNITY
Param ODE ($\omega_{0},\alpha$) 2 trainable parameters $\omega_{0},\lambda$
MLP(in=2, units=200, layers=3, out=2)

#### Optimization hyperparameters

The hyperparameters of the APHYNITY optimization algorithm
($Niter,\lambda_{0},\tau_{1},\tau_{2}$) were cross-validated on the validation
set and are shown in Table 4. All models were trained with a maximum number of
5000 steps with early stopping.

Table 4: Hyperparameters of the damped pendulum experiments.

Method Niter $\lambda_{0}$ $\tau_{1}$ $\tau_{2}$ APHYNITY Hamiltonian 5 1 1
0.1 APHYNITY ParamODE ($\omega_{0}$) 5 1 1 10 APHYNITY ParamODE
($\omega_{0},\lambda$) 5 1000 1 100

##  Appendix F Ablation study

We conduct ablation studies to show the effectiveness of APHYNITY’s adaptive
optimization and trajectory-based learning scheme.

###  F.1 Ablation to vanilla MB/ML cooperation

In Table 5, we consider the ablation case with the vanilla augmentation scheme
found in Le Guen & Thome (2020); Wang et al. (2019); Mehta et al. (2020),
which does not present any proper decomposition guarantee. We observe that the
APHYNITY cooperation scheme outperforms this vanilla scheme in all case, both
in terms of forecasting performances (e.g. log MSE= -0.35 vs. -3.97 for the
Hamiltonian in the pendulum case) and parameter identification (e.g. Err
Param=8.4% vs. 2.3 for Param PDE ($a,b$ for reaction-diffusion). It confirms
the crucial benefits of APHYNITY’s principled decomposition scheme.

###  F.2 Detailed ablation study

We conduct also two other ablations in Table 6:

  * •

derivative supervision: in which $F_{p}+F_{a}$ is trained with supervision
over approximated derivatives on ground truth trajectory, as performed in
Greydanus et al. (2019); Cranmer et al. (2020). More precisely, APHYNITY’s
$\mathcal{L}_{\text{traj}}$ is here replaced with
$\mathcal{L}_{\text{deriv}}=\|\frac{\mathop{}\\!\mathrm{d}X_{t}}{\mathop{}\\!\mathrm{d}t}-F(X_{t})\|$
as in Eq. (D), where
$\frac{\mathop{}\\!\mathrm{d}X_{t}}{\mathop{}\\!\mathrm{d}t}$ is approximated
by finite differences on $X_{t}$.

  * •

non-adaptive optim.: in which we train APHYNITY by minimizing $\|F_{a}\|$
without the adaptive optimization of $\lambda$ shown in Algorithm 1. This case
is equivalent to $\lambda=1,\tau_{2}=0$.

Table 5: Ablation study comparing APHYNITY to the vanilla augmentation scheme
(Wang et al., 2019; Mehta et al., 2020) for the reaction-diffusion equation,
wave equation and damped pendulum.

Dataset Method $\log$ MSE %Err Param. $\|F_{a}\|^{2}$ Reaction-diffusion
Param. PDE ($a,b$) with vanilla aug. -4.56$\pm$0.52 8.4 (7.5$\pm$1.4)e1
APHYNITY Param. PDE ($a,b$) -5.10$\pm$0.21 2.3 (6.7$\pm$0.4)e1 Param. PDE
($a,b,k$) with vanilla aug. -8.04$\pm$0.03 25.4 (1.5$\pm$0.2)e-2 APHYNITY
Param. PDE ($a,b,k$) -9.35$\pm$0.02 0.096 (1.5$\pm$0.4)e-6 True PDE with
vanilla aug. -8.12$\pm$0.05 n/a (6.1$\pm$2.3)e-4 APHYNITY True PDE
-9.17$\pm$0.02 n/a (1.4$\pm$0.8)e-7 Wave equation Param PDE ($c$) with vanilla
aug. -3.90 $\pm$ 0.27 0.51 88.66 APHYNITY Param PDE ($c$) -4.64$\pm$0.25 0.31
71.0 Param PDE ($c,k$) with vanilla aug. -5.96 $\pm$ 0.10 0.71 25.1 APHYNITY
Param PDE ($c,k$) -6.09$\pm$0.28 0.70 4.54 Damped pendulum Hamiltonian with
vanilla aug. -0.35$\pm$0.1 n/a 837$\pm$117 APHYNITY Hamiltonian -3.97$\pm$1.2
n/a 623$\pm$68 Param ODE ($\omega_{0}$) with vanilla aug. -7.02$\pm$1.7 4.5
148$\pm$49 APHYNITY Param ODE ($\omega_{0}$) -7.86$\pm$0.6 4.0 132$\pm$11
Param ODE ($\omega_{0},\alpha$) with vanilla aug. -7.60$\pm$0.6 4.65
35.5$\pm$6.2 APHYNITY Param ODE ($\omega_{0},\alpha$) -8.31$\pm$0.3 0.39
8.5$\pm$2.0 Augmented True ODE with vanilla aug. -8.40$\pm$0.2 n/a 3.4$\pm$0.8
APHYNITY True ODE -8.44$\pm$0.2 n/a 2.3$\pm$0.4

Table 6: Detailed ablation study on supervision and optimization for the
reaction-diffusion equation, wave equation and damped pendulum.

Dataset Method $\log$ MSE %Err Param. $\|F_{a}\|^{2}$ Reaction-diffusion
Augmented Param. PDE ($a,b$) derivative supervision -4.42$\pm$0.25 12.6
(6.8$\pm$0.6)e1 Augmented Param. PDE ($a,b$) non-adaptive optim.
-4.55$\pm$0.11 7.5 (7.6$\pm$1.0)e1 APHYNITY Param. PDE ($a,b$) -5.10$\pm$0.21
2.3 (6.7$\pm$0.4)e1 Augmented Param. PDE ($a,b,k$) derivative supervision
-4.90$\pm$0.06 11.7 (1.9$\pm$0.3)e-1 Augmented Param. PDE ($a,b,k$) non-
adaptive optim. -9.10$\pm$0.02 0.21 (5.5$\pm$2.9)e-7 APHYNITY Param. PDE
($a,b,k$) -9.35$\pm$0.02 0.096 (1.5$\pm$0.4)e-6 Augmented True PDE derivative
supervision -6.03$\pm$0.01 n/a (3.1$\pm$0.8)e-3 Augmented True PDE non-
adaptive optim. -9.01$\pm$0.01 n/a (1.5$\pm$0.8)e-6 APHYNITY True PDE
-9.17$\pm$0.02 n/a (1.4$\pm$0.8)e-7 Wave equation Augmented Param PDE ($c$)
derivative supervision -1.16$\pm$0.48 12.1 0.00024 Augmented Param PDE ($c$)
non-adaptive optim. -2.57$\pm$0.21 3.1 43.6 APHYNITY Param PDE ($c$)
-4.64$\pm$0.25 0.31 71.0 Augmented Param PDE ($c,k$) derivative supervision
-4.19$\pm$0.36 7.2 0.00012 Augmented Param PDE ($c,k$) non-adaptive optim.
-4.93$\pm$0.51 1.32 0.054 APHYNITY Param PDE ($c,k$) -6.09$\pm$0.28 0.70 4.54
Augmented True PDE derivative supervision -4.42 $\pm$ 0.33 n/a 6.02e-5
Augmented True PDE non-adaptive optim. -4.97$\pm$0.49 n/a 0.23 APHYNITY True
PDE -5.24$\pm$0.45 n/a 0.14 Damped pendulum Augmented Hamiltonian derivative
supervision -0.83$\pm$0.3 n/a 642$\pm$121 Augmented Hamiltonian non-adaptive
optim. -0.49$\pm$0.58 n/a 165$\pm$30 APHYNITY Hamiltonian -3.97$\pm$1.2 n/a
623$\pm$68 Augmented Param ODE ($\omega_{0}$) derivative supervision
-1.02$\pm$0.04 5.8 136$\pm$13 Augmented Param ODE ($\omega_{0}$) non-adaptive
optim. -4.30$\pm$1.3 4.4 90.4$\pm$27 APHYNITY Param ODE ($\omega_{0}$)
-7.86$\pm$0.6 4.0 132$\pm$11 Augmented Param ODE ($\omega_{0},\alpha$)
derivative supervision -2.61$\pm$0.2 5.0 3.2$\pm$1.7 Augmented Param ODE
($\omega_{0},\alpha$) non-adaptive optim. -7.69$\pm$1.3 1.65 4.8$\pm$7.7
APHYNITY Param ODE ($\omega_{0},\alpha$) -8.31$\pm$0.3 0.39 8.5$\pm$2.0
Augmented True ODE derivative supervision -2.14$\pm$0.3 n/a 4.1$\pm$0.6
Augmented True ODE non-adaptive optim. -8.34$\pm$0.4 n/a 1.4$\pm$0.3 APHYNITY
True ODE -8.44$\pm$0.2 n/a 2.3$\pm$0.4

We highlight the importance to use a principled adaptive optimization
algorithm (APHYNITY algorithm described in paper) compared to a non-adpative
optimization: for example in the reaction-diffusion case, log MSE= -4.55 vs.
-5.10 for Param PDE $(a,b)$. Finally, when the supervision occurs on the
derivative, both forecasting and parameter identification results are
systematically lower than with APHYNITY’s trajectory based approach: for
example, log MSE=-1.16 vs. -4.64 for Param PDE $(c)$ in the wave equation. It
confirms the good properties of the APHYNITY training scheme.

##  Appendix G Additional experiments

###  G.1 Reaction-diffusion systems with varying diffusion parameters

We conduct an extensive evaluation on a setting with varying diffusion
parameters for reaction-diffusion equations. The only varying parameters are
diffusion coefficients, i.e. individual $a$ and $b$ for each sequence. We
randomly sample $a\in[1\times 10^{-3},2\times 10^{-3}]$ and $b\in[3\times
10^{-3},7\times 10^{-3}]$. $k$ is still fixed to $5\times 10^{-3}$ across the
dataset.

In order to estimate $a$ and $b$ for each sequence, we use here a ConvNet
encoder $E$ to estimate parameters from 5 reserved frames ($t<0$). The
architecture of the encoder $E$ is similar to the one in Table 2 except that
$E$ takes 5 frames (10 channels) as input and $E$ outputs a vector of
estimated $(\tilde{a},\tilde{b})$ after applying a sigmoid activation scaled
by $1\times 10^{-2}$ (to avoid possible divergence). For the baseline Neural
ODE, we concatenate $a$ and $b$ to each sequence as two channels.

In Table G.1, we observe that combining data-driven and physical components
outperforms the pure data-driven one. When applying APHYNITY to Param PDE
($a,b$), the prediction precision is significantly improved ($\log$ MSE: -1.32
vs. -4.32) with $a$ and $b$ respectively reduced from 55.6% and 54.1% to 11.8%
and 18.7%. For complete physics cases, the parameter estimations are also
improved for Param PDE ($a,b,k$) by reducing over 60% of the error of $b$
(3.10 vs. 1.23) and 10% to 20% of the errors of $a$ and $k$ (resp. 1.55/0.59
vs. 1.29/0.39).

The extensive results reflect the same conclusion as shown in the main
article: APHYNITY improves the prediction precision and parameter estimation.
The same decreasing tendency of $\|F_{a}\|$ is also confirmed.

Table 7: Results of the dataset of reaction-diffusion with varying $(a,b)$.
$k=5\times 10^{-3}$ is shared across the dataset.

| Method |  $\log$ MSE | %Err $a$ | %Err $b$ | %Err $k$ | $\|F_{a}\|^{2}$
---|---|---|---|---|---|---
Data-driven | Neural ODE (Chen et al., 2018) | -3.61$\pm$0.07 | n/a | n/a | n/a | n/a
Incomplete physics | Param PDE ($a,b$) | -1.32$\pm$0.02 | 55.6 | 54.1 | n/a | n/a
APHYNITY Param PDE ($a,b$) | -4.32$\pm$0.32 | 11.8 | 18.7 | n/a | (4.3$\pm$0.6)e1
Complete physics | Param PDE ($a,b,k$) | -5.54$\pm$0.38 | 1.55 | 3.10 | 0.59 | n/a
APHYNITY Param PDE ($a,b,k$) | -5.72$\pm$0.25 | 1.29 | 1.23 | 0.39 | (5.9$\pm$4.3)e-1
2-7 | True PDE | -8.86$\pm$0.02 | n/a | n/a | n/a | n/a
| APHYNITY True PDE | -8.82$\pm$0.15 | n/a | n/a | n/a | (1.8$\pm$0.6)e-5

###  G.2 Additional results for the wave equation

We conduct an experiment where each sequence is generated with a different
wave celerity. This dataset is challenging because both $c$ and the initial
conditions vary across the sequences. For each simulated sequence, an initial
condition is sampled as described previously, along with a wave celerity $c$
also sampled uniformly in $[300,400]$. Finally our initial state is integrated
with the same Runge-Kutta scheme. $200$ of such sequences are generated for
training while $50$ are kept for testing.

For this experiment, we also use a ConvNet encoder to estimate the wave speed
$c$ from 5 consecutive reserved states $(w,\frac{\partial w}{\partial t})$.
The architecture of the encoder $E$ is the same as in Table 2 but with 10
input channels. Here also, $k$ is fixed for all sequences and $k=50$. The
hyper-parameters used in these experiments are the same than described in the
Section E.2.

The results when multiple wave speeds $c$ are in the dataset are consistent
with the one present when only one is considered. Indeed, while prediction
performances are slightly hindered, the parameter estimation remains
consistent for both $c$ and $k$. This extension provides elements attesting
for the robustness and adaptability of our method to more complex settings.
Finally the purely data-driven Neural-ODE fails to cope with the increasing
difficulty.

Table 8: Results for the damped wave equation when considering multiple $c$
sampled uniformly in $[300,400]$ in the dataset, $k$ is shared across all
sequences and $k=50$.

| Method |  $\log$ MSE | %Error $c$ | %Error $k$ | $\|F_{a}\|^{2}$
---|---|---|---|---|---
Data-driven | Neural ODE | 0.056$\pm$0.34 | n/a | n/a | n/a
Incomplete physics | Param PDE ($c$) | -1.32$\pm$0.27 | 23.9 | n/a | n/a
APHYNITY Param PDE ($c$) | -4.51$\pm$0.38 | 3.2 | n/a | 171
Complete physics | Param PDE ($c,k$) | -4.25$\pm$0.28 | 3.54 | 1.43 | n/a
APHYNITY Param PDE ($c,k$) | -4.84$\pm$0.57 | 2.41 | 0.064 | 3.64
2-6 | True PDE ($c,k$) | -4.51$\pm$0.29 | n/a | n/a | n/a
| APHYNITY True PDE ($c,k$) | -4.49$\pm$0.22 | n/a | n/a | 0.0005

###  G.3 Damped pendulum with varying parameters

To extend the experiments conducted in the paper (section 4) with fixed
parameters ($T_{0}=6,\alpha=0.2$) and varying initial conditions, we evaluate
APHYNITY on a much more challenging dataset where we vary both the parameters
($T_{0},\alpha$) and the initial conditions between trajectories.

We simulate 500/50/50 trajectories for the train/valid/test sets integrated
with DOPRI5. For each trajectory, the period $T_{0}$ (resp. the damping
coefficient $\alpha)$ are sampled uniformly in the range $[3,10]$ (resp.
$[0,0.5]$).

We train models that take the first 20 steps as input and predict the next 20
steps. To account for the varying ODE parameters between sequences, we use an
encoder that estimates the parameters based on the first 20 timesteps. In
practice, we use a recurrent encoder composed of 1 layer of 128 GRU units. The
output of the encoder is fed as additional input to the data-driven
augmentation models and to an MLP with final softplus activations to estimate
the physical parameters when necessary ($\omega_{0}\in\mathbb{R}_{+}$ for
Param ODE ($\omega_{0}$), $(\omega_{0},\alpha)\in\mathbb{R}_{+}^{2}$ for Param
ODE ($\omega_{0},\alpha$)).

In this varying ODE context, we also compare to the state-of-the-art
univariate time series forecasting method N-Beats (Oreshkin et al., 2020).

Results shown in Table G.3 are consistent with those presented in the paper.
Pure data-driven models Neural ODE (Chen et al., 2018) and N-Beats (Oreshkin
et al., 2020) fail to properly extrapolate the pendulum dynamics. Incomplete
physical models (Hamiltonian and ParamODE ($\omega_{0}$)) are even worse since
they do not account for friction. Augmenting them with APHYNITY significantly
and consistently improves forecasting results and parameter identification.

Table 9: Forecasting and identification results on the damped pendulum dataset
with different parameters for each sequence. log MSEs are computed over 20
predicted time-steps. For each level of incorporated physical knowledge,
equivalent best results according to a Student t-test are shown in bold. n/a
corresponds to non-applicable cases.