Year: 2021

Abstract:
Modeling dynamical systems combining prior physical knowledge and machine learning (ML) is promising in scientific problems when the underlying processes are not fully understood, e.g. when the dynamics is partially known. A common practice to identify the respective parameters of the physical and ML components is to formulate the problem as supervised learning on observed trajectories. However, this formulation leads to an infinite number of possible decompositions. To solve this ill-posedness, we reformulate the learning problem by introducing an upper bound on the prediction error of a physical-statistical model. This allows us to control the contribution of both the physical and statistical components to the overall prediction. This framework generalizes several existing hybrid schemes proposed in the literature. We provide theoretical guarantees on the wellposedness of our formulation along with a proof of convergence in a simple affine setting. For more complex dynamics, we validate our framework experimentally.

heading: 1 INTRODUCTION
text: Dynamical systems prediction and identification find crucial applications ranging from medicine and the study of tumors (Hanahan & Weinberg, 2011; Lu & Fei, 2014) to oceanic and climate forecasting (Oreskes et al., 1994; Caers, 2011). The modeling of such systems traditionally rely on ordinary or partial differential equations (ODE/PDE) (Madec, 2008; Marti et al., 2010), and their resolution via numerical solvers and data assimilation (Ghil & Malanotte-Rizzoli, 1991). In real world applications, two main pitfalls occur: first the dynamics may only be partially known and thus do not fully represent the studied phenomena (Rio & Santoleri, 2018); second, the system state may only be partially observed as in ocean models (Gaultier et al., 2013). Machine learning (ML) has become a complementary approach to traditional physics based models (denoted MB for model based) (Reichstein et al., 2019; Dueben & Bauer, 2018). Both offer advantages: whereas MB approaches generalize and extrapolate better, ML high expressivity approaches benefit from the ongoing explosive growth of available data such as satellite observations, with reduced costs compared to classical data assimilation.
In that perspective, recent lines of work tackle the learning of hybrid models relying on prior physical knowledge and machine learning (Yin et al., 2021; Mehta et al., 2020). Efficiently learning such decompositions actually means solving two different tasks: system identification, i.e. estimating the parameters of the physical model, and prediction, i.e. recovering the trajectories associated to the dynamics. Both are essential for hybrid MB/ML models of dynamical systems. Whereas prediction aims at robust extrapolation, identification accounts for physical interpretability of the MB/ML model. While solving both problems using model-based formulation admits well-known numerical solutions, for example using the adjoint method (Le Dimet & Talagrand, 1986; Courtier et al., 1994), the combination of physical models and deep learning is still an open area of research. In this context, ML applications mainly focus on the prediction task, at the expense of the system identification: Willard et al. (2020) underlines the lack of generalizability of black-box ML models and their inability to produce physically sound results. Indeed, Ayed et al. (2020) show that without any prior knowledge, the recovered estimates of a dynamical system states are not physically plausible despite accurate predictions. Moreover, as noted by Yin et al. (2021), learning a linear MB/ML decompo-
sition with the sole supervision on the system trajectories is ill-posed and admits an infinite number of decompositions. Such observations highlight the need to incorporate physically motivated constraints in the learning of hybrid models, e.g. through regularization penalties, and several works propose additional constraints to guide the model towards physical solutions (Jia et al., 2019; Yin et al., 2021; Linial et al., 2021). Finally, to complete prior dynamical knowledge with a data-driven component and ensure interpretability of the decomposition, we work out a principled framework that generalizes previous attempts in the regularization of hybrid models. Our contributions are :
• In section 3.1, we introduce a novel way to recover well-posedness and interpretability in the learning of hybrid MB/ML models through the control of an upper bound. We further extend our framework to incorporate auxiliary data when available to handle complex realworld problems.
• In section 3.2, we propose a novel alternate-optimization algorithm to learn hybrid models. • In section 3.3, we provide an analysis of the convergence of the proposed algorithm on a
simplified case and experimentally evidence the soundness of our approach on more complex settings of increasing difficulty including challenging real world problems (section 4).

heading: 2 BACKGROUND AND PROBLEM SETUP
text: We consider a dynamical system with state at time t denoted Zt = Z(t). Zt might be fully or only partially observed: we write Zt = (Xt, Yt), where Xt is the observed component and Yt the unobserved one. The evolution of Z is governed by a differential equation with dynamics :
dZt dt = d dt ( Xt Yt ) = ( fX(Zt) fY (Zt) ) (1)
The objective is to predict trajectories of X , i.e. to model the evolution of the observable part following dXtdt = fX(Zt). For simplicity, we omit the index X in fX and write f(.) ∆ = fX(.).
Dynamical Hypothesis We assume partial knowledge of the dynamics of the observed Xt:
dXt dt = f(Zt) = fk(Zt) + fu(Zt) (2)
where fk ∈ Hk is a known operator with unknown parameters θ∗, and fu ∈ Hu is the unknown residual dynamics. Hk and Hu denote function spaces. Note that the additive hypothesis in eq. (2) is not restrictive (see discussion in appendix B).
Learning Problem Our objective is to approximate f with a function h learned from the observed data. According to eq. (2), we assume h = hk +hu. hk ∈ Hk belongs to the same hypothesis space as fk, i.e. it has the same parametric form. Its parameters are denoted θk. Note that hk(., θ∗) = fk. hu ∈ Hu is represented by a free form functional with parameters θu, e.g. a neural network. The learning problem is then to estimate from data the parameters of hk so that they match the true physical ones and those of hu so as to approximate at best the unknown dynamics f . In this regard, an intuitive training objective is to minimize a distance d between h = hk + hu and f :
d(h, f) = EZ∼pZ‖h(Z)− f(Z)‖2, (3)
where pZ is the distribution of the state values Z that accounts for varying initial states. Each Z then defines a training sample. Minimizing eq. (3) with h = hk + hu enables to predict accurate trajectories but may have an infinite number of solutions. For instance, hu may bypass the physical hypothesis hk. Thus, interpretability is not guaranteed. We develop in the following section a method to overcome this ill-posedness.

heading: 3 METHOD
text: In hybrid modeling, two criteria are essentials: 1. identifiability, i.e. the estimated parameters of hk should correspond to the true physical ones; 2. prediction power, i.e. the statistical component hu should complete hk so that h = hk + hu performs accurate prediction over the system states. To
control the contribution of each term hk and hu, we work upper bounds out of eq. (3) (section 3.1). We then propose to minimize d(h, f) while constraining the upper bounds, which provide us with a well-posed learning framework (section 3.2). Besides, we show that several previous works that introduced constrained optimization to solve related problems are specific cases of our formulation (Yin et al., 2021; Jia et al., 2019; Linial et al., 2021). Finally, we introduce an alternate optimization algorithm which convergence is shown in section 3.3 for a linear approximation of f .

heading: 3.1 STRUCTURAL CONSTRAINTS FOR DYNAMICAL SYSTEMS
text: To ensure identifiability, we derive regularizations on hk and hu flowing from the control of an upper bound of d(h, f). In particular, to minimize d(hk, fk) would enable us to accurately interpret hk as the true fk, and thus hu as the residual dynamics fu. However, since we do not access the parameters of fk, computing d(hk, fk) is not tractable. We then consider two possible situations. In the first one, the only available information on the physical system is the parametric form of fk (or equivalently of hk), training thus only relies on observed trajectories (section 3.1.1). In the second one, we consider that auxiliary information about fk is available and will be used to minimize the distance between hk and fk (section 3.1.2). While the first setting is the more general, the physical prior it relies on is often insufficient to effectively handle real world situations. The second setting makes use of more informative priors and better corresponds to real cases as will be shown in the experimental section (see section 4.2).

heading: 3.1.1 CONTROLLING THE ML COMPONENT AND THE MB HYPOTHESIS
text: We propose a general approach to constrain the learning of hybrid models when one solely access the functional form of hk. In this case, to make hk accountable in our observed phenomena, a solution is to minimize d(hk, f). Following the triangle inequality we link up both errors d(h, f) and d(hk, f) (computations available in appendix C.1):
d(h, f) ≤ d(h, hk) + d(hk, f) = d(hu, 0) + d(hk, f) (4)
We want the physical-statistical model h = hk + hu to provide high quality forecasts. Minimizing the sole upper bound does not ensure such aim, as hu is only penalized through d(hu, 0) and is not optimized to contribute to predictions. We thus propose to minimize d(h, f) while controlling both d(hu, 0) and d(hk, f). Such a control of the upper bound of eq. (4) amounts to balancing the contribution of the ML and the MB components. This will be formally introduced in section 3.2.
Link to the Literature The least action principle on the ML component i.e. constraining d(hu, 0) is invoked for a geometric argument in (Yin et al., 2021), and appears as a co-product of the introduction of d(hk, f) in eq. (4). Optimizing d(hk, f) to match the physical model with observations is investigated in (Forssell & Lindskog, 1997).
The general approach of eq. (4) allows us to perform prediction (via h) and system identification (via hk) on simple problems (see section 4.1). The learning of real-world complex dynamics, via data-driven hybrid models, often fails at yielding a physically sound estimation, as illustrated in section 4.2. This suggests that learning complex dynamics requires additional information. In many real-world cases, auxiliary information is available in the form of measurements providing complementary information on fk. Indeed, a common issue in physics is to infer an unobserved variable of interest (in our case fk parameters θ?) from indirect or noisy measurements that we refer to as proxy data. For instance, one can access a physical quantity but only at a coarse resolution, as in (Um et al., 2020; Belbute-Peres et al., 2020) and in the real world example detailed in section 4.2. We show in the next subsection how to incorporate such an information in order to approximate d(hk, fk).

heading: 3.1.2 MATCHING THE PHYSICAL HYPOTHESES: INTRODUCING AUXILIARY DATA
text: We here assume one accesses a proxy of fk, denoted f pr k ∈ Hk. Our goal is to adapt our framework to incorporate such auxiliary information, bringing the regularization induced by fprk within the scope of the control of an upper bound. This enables us to extend our proposition towards the solving of real world physical problems, still largely unexplored by the ML community. We have:
d(h, f) ≤ d(h, hk) + d(hk, fprk ) + Γ = d(hu, 0) + d(hk, f pr k ) + Γ (5)
where Γ is a constant of the problem that cannot be optimized (see appendix C.2). In that context, we can benefit from auxiliary information providing us with coarse estimates of θ?, denoted θpr, such that fprk = hk(., θ
pr) ≈ fk. To use the available θpr to guide our estimation towards the true parameters θ? of fk, a simple solution is to directly enforce the minimization of d(hk, f pr k ) in the parameter space by minimizing ‖θk − θpr‖2, where θk are the parameters of hk. Indeed, because fk and f pr k have identical parametric forms (as both belong to the same functional space Hk), minimizing ‖θk − θpr‖2 will bring hk closer to fprk and thus to fk. As above, we propose to minimize d(h, f) while controlling both d(hu, 0) and d(hk, f pr k ), as described in section 3.2.
Link to the Literature In (Linial et al., 2021) fprk stands for true observations used to constrain a learned latent space, minimizing d(hk, f pr k ). Jia et al. (2019) uses synthetic data as f pr k to pre-train their model which amounts to the control an upper bound, see appendix C.3. Finally, this setting finds an extension, when the model fprk is a learned model, for example trained using eq. (4), leading to a self-supervision approach described in appendix C.4.

heading: 3.2 LEARNING ALGORITHM AND OPTIMIZATION PROBLEM
text: From the upper bounds, we first recover the well-posedness of the optimization and derive a theoretical learning scheme (section 3.2.1). We then discuss its practical implementation (section 3.2.2).

heading: 3.2.1 WELL-POSEDNESS AND ALTERNATE OPTIMIZATION ALGORITHM
text: Recovering Well-Posedness We reformulate the ill-posed learning of minhk,hu∈Hk×Hu d(h, f), by instead optimizing d(h, f) while constraining the upper bounds. Let us define Sk and Su as
Sk = {hk ∈ Hk | `(hk) ≤ µk } Su = {hu ∈ Hu | d(hu, 0) ≤ µu } (6)
where µk, µu are two positive scalars and `(hk) = d(hk, f) in the case of section 3.1.1 and `(hk) = d(hk, f pr k ) in the case of section 3.1.2. Our proposition then amounts to optimizing d(h, f) over the Minkowski-sum Sk + Su = {h = hk + hu | hk ∈ Sk, hu ∈ Su } :
min h∈Sk+Su d(h, f), (7)
This constrained optimization setting enables us to recover the well-posedness of the optimization problem under the relative compactness of the family of functionHk (proof in appendix D.3). Proposition 1 (Well-posedness). Under the relative compactness of Sk, eq. (7) finds a solution h that writes as h = hk + hu ∈ Sk + Su. Moreover, this solution is unique.
Alternate Optimization Algorithm As the terms in both upper bounds of eqs. (4) and (5) specifically address either hk or hu, we isolate losses relative to hk and hu and alternate projections of hk on Sk and hu on Su, as described in Algorithm 1. Said otherwise, we learn h by alternately optimizing hk (hu being fixed) and hu (hk being fixed). In practice, we rely on a dual formulation (see section 3.2.2 and the SGD version of Algorithm 1 in Appendix F).
Algorithm 1 Alternate estimation: General Setting Result: Converged hk and hu Set h0u = 0, h 0 k = minhk∈Hk d(hk, f), tol ∈ R+
while d(h, f) > tol do hn+1k = arg min
hk∈Sk d(hk + h
n u, f); h n+1 u = arg min hu∈Su d(hn+1k + hu, f) (8)
n← n+ 1 end
The convergence of the alternate projections is well studied for the intersection of convex sets or smooth manifolds (von Neumann, 1950; Lewis & Malick, 2008) and has been extended in our setting of Minkowski-sum of convex sets (Lange et al., 2019). Because d as defined in eq. (3) is convex, Su and Sk are convex sets as soon asHk andHu are convex (Appendix A). Thus, if d(., f) is strongly convex, eq. (8) finds one and only one solution (Boyd et al., 2004). However, neither the
convexity of Hu nor of Hk is practically ensured. Nonetheless, we recover the well-posedness of eq. (7) and show the convergence of Algorithm 1 in the simplified case where h is an affine function of Xt (see section 3.3). For complex PDE where convexity may not hold, we validate our approach experimentally and we evidence in section 4 that this formulation enables us to recover both an interpretable decomposition h = hk +hu and improved prediction and identification performances.

heading: 3.2.2 PRACTICAL OPTIMIZATION
text: Equation (6) involves the choice of µk and µu. In practice, we implement the projection algorithm by descending gradients on the parameters of hk and hu, with respect to the following losses: Lk(hk) = λhd(h, f) + λhk`(hk) Lu(hu) = λhd(h, f) + λhud(hu, 0) (9) where λh, λhk , λhu are positive real values, dynamically increased/decreased during training. Indeed, d(hu, 0) can be interpreted as a stability loss, preventing the neural networks to trump the physical component. On the other hand, d(hk, f) can be interpreted has an initialization loss yield a first estimate of θk explaining the dynamics. Yet, f being unknown: d(h, f) is not tractable. To estimate d(h, f), we rely on the trajectories associated to the dynamics. We minimize the distance between the ODE flows φh and φf defined by h and f , dφ(φh, φf ), over all initial conditions X0:
dφ(φh, φf ) = EX0 ∫ t t0 ‖φh(τ,X0)− φf (τ,X0)‖2dτ (10)
We have: dφ(φh, φf ) = 0 ⇔ d(h, f) = 0. Definitions of flows for ODE and in depth consideration on these distances are available in appendix A. The gradients of dφ(φh, φf ) with respect to the parameters of hk or hu can be either estimated analytically using the adjoint method (Chen et al., 2018) or using explicit solvers, e.g. Rk45, and computing the gradients thanks to the backpropagation, see (Onken & Ruthotto, 2020). To compute eq. (10), we rely on a temporal sampling of X: our datasets are composed of n sequences of observations of length N , Xi = (Xit0 , . . . , X i t0+N∆t
), where each sequence Xi follows eq. (2) and corresponds to one initial condition Xit0 . We then sample the space of initial conditions Xit0 to compute a Monte-Carlo approximation to dφ(φh, φf ). Let ODESolve be the function integrating any arbitrary initial state xt0 up to time t with dynamics h, so that xt = ODESolve(xt0 , h, t). The estimate of dφ(φh, φf ) then writes as:
dφ(φh, φf ) ≈ 1
n n∑ i=1 N∑ j=1 ∥∥∥ODEsolve(Xit0 , h, tj)−Xitj∥∥∥ 2
Note that the way to compute ODEsolve differs across the experiments (see section 4).

heading: 3.3 THEORETICAL ANALYSIS FOR A LINEAR APPROXIMATION
text: We investigate the validity of our proposition when approximating an unknown derivative with an affine function (interpretable first guess approximators). We here consider hk as a linear function. We do not assume any information on f , thus relieving this section from the need of an accurate prior knowledge fk. In this context, we show the convergence of the learning scheme introduced in Algorithm 1 with ` = d(hk, f), hence demonstrating the validity of our framework in this simplified setting. For more complex cases, for which theoretical analysis cannot be conducted, our framework is validated experimentally in section 4. All proofs of this section are conducted using the distance dφ. Let Xs be the unique solution to the initial value problem:
dXt dt = f(Xt) with Xt=0 = X0 (11)
With hk(X) = AX and hu(X) = DA, the affine approximation of f writes as: dXt dt = AXt +DA with Xt=0 = X0 (12) where A ∈ Mp,p(R), DA ∈ Rp. We write XD the solution to eq. (12) and XA the solution to eq. (12) when DA = 0. The alternate projection algorithm with the distance dφ writes as:
Â = arg min A ∫ t t0 ∥∥Xs(τ)−XD(τ)∥∥ 2 dτ + λA ∫ t t0 ∥∥Xs(τ)−XA(τ)∥∥ 2 dτ (13)
D̂A = arg min DA ∫ t t0 ∥∥Xs(τ)−XD(τ)∥∥ 2 dτ + λD‖DA‖2 (14)
where λD, λA > 0. As the optimization of eq. (13) is not convex on A, the solution existence and uniqueness is not ensured. The well-posedness w.r.t A can be recovered by instead considering a simple discretization scheme, e.g. Xt+1 ≈ (AXt + DA)∆t + Xt and solving the associated least square regression, which well-posedness is guaranteed, see details in appendix D.2. Such strategy is common practice in system identification. Generic theoretical considerations on existence and uniqueness of solutions to eqs. (13) and (14) are hard to retrieve. Nonetheless, if A is an invertible matrix, we prove in appendix D.4:
Proposition 2 (Existence and Uniqueness). If Â is invertible, There exists a unique DA, hence a unique XD, solving eq. (14).
Finally, formulating Algorithm 1 as a least square problem in an affine setting (see appendix D.5), we prove the convergence of the alternate projection algorithm (appendix D.6) : Proposition 3. For λD and λA sufficiently high, the algorithm that alternates between the estimation of A and the estimation of DA following eqs. (13) and (14) converges.

heading: 4 EXPERIMENTS
text: We validate Algorithm 1 on datasets of increasing difficulty (see appendix E), where the system state is either fully or partially observed (resp. section 4.1 and section 4.2). We no longer rely on an affine prior and explicit hk and hu for each dataset. Performances are evaluated via standard metrics: MSE (lower is better) and relative Mean Absolute Error (rMAE, lower is better). We assess the relevance of our proposition based on eqs. (4) and (5), against NeuralODE (Chen et al., 2018) and state of the art Aphynity (Yin et al., 2021) and ablation studies. We denote Ours eq. (4) (resp. Ours eq. (5)) the results when ` = d(hk, f) i.e eq. (4), (resp. ` = d(hk, f pr k ) i.e. eq. (5)) When d(hk, f) (resp. d(hu, 0)) is not considered in the optimization, we refer to the results as d(h, f)+d(hu, 0) (resp. d(h, f)+d(hk, f)). When h is trained by only minimizing the discrepancy between actual and predicted trajectories the results are denoted «Only d(h, f)». We report between brackets the standard deviation of the metrics over 5 runs and refer to Appendices F and G for training information and additional results. Code and datasets will be provided upon acceptance.

heading: 4.1 FULLY OBSERVABLE DYNAMICS
text: To illustrate the learning scheme induced by eq. (4), we focus on fully observed low dimensional dynamics: a simple example emerging from Newtonian mechanics and a population dynamics model.
Damped Pendulum (DPL) Now a standard benchmark for hybrid models, we consider the motion of a pendulum of length L damped due to viscous friction (Greydanus et al., 2019; Yin et al., 2021). Newtonian mechanics provide an ODE describing the evolution of the angle x of the pendulum:
ẍ− g/L sin(x) + kẋ = 0 (15) We suppose access to observations of the system state Z = (x, ẋ). We consider as physical motion hypothesis hk(x, θk) = θk sin(x). The true pulsation θ∗ = g/L of the pendulum has to be estimated with θk. The viscous friction term kẋ remains to be estimated by hu.
Population Dynamics (LV) Lotka-Volterra ODE system models a prey/predator population dynamics describing the growth of the preys (x) without predators (y), and the extinction of predators without preys, the non linear terms expressing the encounters between both species:
ẋ = αx− βxy, and ẏ = −γy + δxy (16) We observe the system state Z = (x, y) and set as prior knowledge: hk(x, y) = (θ1kx,−θ2ky). θ? = (α, γ) has to be estimated by θk = (θ1k, θ 2 k). hu accounts for the non linear terms (βxy, δxy).
Experimental Setting For both DPL and LV experiments, we consider the following setting: we sample the space of initial conditions building 100/50/50 trajectories for the train, validation and test sets. The sequences share the same parameters; respectively ( gL , k), for DPL, and (α, β, γ, δ) for LV. The parameter θk is set to a neuron (of dimension 1 in the pendulum and 2 for LV) and hu is a 2-layer MLP. Further experimental details are available in appendices E.1, E.2 and F.
Table 1: Experimental Results for PDL and LV data. The presented metric for parameter evaluation is the rMAE reported in %. Pred. columns report the prediction log MSE on trajectories on test set.
Model PDL LVrMAE(θk, θ?) Pred. logMSE rMAE(θk, θ?) Pred. logMSE
Ours eq. (4) 1.56 (0.009) -13.7 (0.84) 7.80 (0.011) -9.28 (0.75) Only d(h, f) 9.35 (0.04) -13.3 (0.65) 24.5 (0.017) -9.21 (0.91) d(h, f) + d(hk, f) 1.82 (0.01) -13.4 (0.56) 7.91 (0.02) -9.01 (0.99) d(h, f) + d(hu, 0) 11.1 (0.03) -12.9 (0.29) 9.80 (0.098) -9.45 (0.55) Aphynity 6.15 (0.009) -12.2 (0.13) 21.1 (0.016) -9.89 (0.53)
Identification and Prediction Results Table 1 shows that despite accurate trajectory forecasting, the unconstrained setting «Only d(h, f)» fails at estimating the models parameters, showing the need for regularization for identification. Constraining the norm of the ML component can be insufficient: for LV data, both Aphynity and d(h, f) + d(hu, 0) do not accurately estimate the model parameters. However, the control of d(hk, f), following eq. (4), significantly improves the parameter identification for both datasets. Indeed, in the PDL case, hk and f are (pseudo)-periodic of the same period, hence the gain in the performances. Finally, our proposition based on eq. (4) is able to identify the parameters of DPL and LV equation with a precision of respectively 1.56% and 7.8% beating all consid-
ered baselines. Regarding prediction performances, in under-constrained settings ( «Only d(h, f)» in Table 1), hu learns to corrects the inaccurate hk. Table 1 and figs. 4 and 5 (appendix G.1) show that our proposition provides more consistent prediction performances. These experiments confirm that the constraints on hk and hu arising from the control of the upper bound of eq. (4) increase interpretability and maintain prediction performances.
Throwback to the Affine Case We verify the convergence proved in section 3.3 using the damped pendulum (eq. (15)) linearized in the small oscillations regime (see appendix E.1). Making an affine hypothesis following eq. (12), we apply our alternate projection algorithm and optimize A and DA alternately using SGD. Figure 1 shows that we are able to accurately estimate A and D using our proposition, recovering both the oscillation pulsation and the damping coefficient.

heading: 4.2 HIGH DIMENSIONAL DYNAMICS
text: We now address the learning of transport equations, describing a wide range of physical phenomena such as chemical concentration, fluid dynamics or material properties. We evaluate the learning setting induced by eq. (4) and (5) on two physical datasets depicting the evolution of the temperature T advected by a time-dependent velocity field U and subject to forcing S, following:
∂T ∂t +∇.(TU) = S(U) (17)
The system state Z = (T,U, S) is partially observed, we only access T . Every quantities, observed or to estimate, are regularly sampled on a spatiotemporal grid: at each timestep t, the time varying velocity field Ut writes as Ut = (ut, vt) and ut, vt, Tt and the forcing term St are all of size 64×64.
Experimental Setting We consider as physical prior the advection i.e hk(T, θk) = −∇.(Tθk). Thus, θk is time-dependent, as we learn it to approximate θ? = U . We identify the velocity field θk from observations of T , learning a mapping between T and U parameterized by a neural network Gψ , so that θk = Gψ(Tt−l, ..., Tt) ≈ Ut, which is common practice in oceanography (Béréziat & Herlin, 2015). Gψ is optimized following eq. (9). S remains to be learned by hu. hk implements a differentiable semi-Lagrangian scheme (Jaderberg et al., 2015) (see appendix E.3) and hu is a ResNet. Gψ is a UNet. Training details and a schema of our model are to be found in appendix F.
Table 2: Results for Adv+S and Natl data. We report the MSE (× 100) on the predicted observations T , the velocity fields U and the source term S over 6 time steps on test set.
Models Adv+S Natl
T U S T U S
Ours eq. (4) 0.74 (0.05) 1.99 (0.13) 0.17 (0.01) 8.27 (0.06) 11.72 (0.07) 6.01 (0.08) Ours eq. (5) – – – 6.86 (0.12) 6.81 (0.07) 4.35 (0.11) Aphynity 0.85 (0.35) 3.07 (0.74) 0.18 (0.05) 8.18 (0.16) 11.75 (0.49) 6.02 (0.02) NODE 1.35 (0.02) – – 8.83 (0.98) – –
Figure 2: Best viewed in color. Estimations of S, T and U = (u, v) on Adv+S. Prediction ranges from 1 to 20 half-days.
Synthetic Advection and Source (Adv+S) To test the applicability of the learning setting induced by eq. (4) on partially observed settings, we first study a synthetic setting (denoted Adv+S) of eq. (17) by generating velocity fields U , simulated following (Boffetta et al., 2001) and adding a source term S inspired by (Frankignoul, 1985). The simulation details are given in appendix E.3.
Real Ocean Dynamics (Natl) We consider a complex dataset emulating real world observations of the North ATLantic ocean (denoted Natl) (Ajayi et al., 2019). Modeling the evolution of T in Natl is particularly challenging as its dynamics is chaotic and highly non-linear. This simulation is representative of the complexity encountered in real world data. The principled approach corresponding to eq. (4) is insufficient here and one must resort to additional physical information. To illustrate the extension developed in section 4.2,
we make use of available auxiliary data: measurements from satellite observations provide a coarse estimate of the surface ocean currents velocity fields (see appendix E). Then, the learning goal is to refine the approximated velocity fields to fit the ocean dynamics. Therefore we proceed as described in eq. (5) and enforce d(hk, f pr k ) by supervising Gψ with the proxy data (see appendix E.3).
Identification and Prediction Results Table 2 indicates that for Adv+S dataset, we estimate accurately the unobserved velocity fields. Qualitatively, Figure 2 shows that controlling our proposed upper bound eq. (4) facilitates the recovery of truthful velocity fields U along with an accurate prediction of T . For the highly complex Natl, Table 2 shows that the introduction of auxiliary data following the formulation in eq. (5) significantly helps identification, as the dynamics is too complex to be able to recover physically interpretable velocity fields using the bound of eq. (4). Regarding prediction performances on the Adv+S data, Table 2 shows that thanks to our truthful estimates of U , our model provides more precise prediction than NODE and Aphynity. For real world data, thanks to the proxy data our model recovers better velocity fields terms while providing a better estimate for T . Besides, adding prior knowledge in the prediction systems improves prediction performances: appendix G shows that NODE minimizes d(h, f) by predicting average and blurred frames. This shows the need for regularization when learning on structured physical data.
Ablation Study We present in Table 3 an ablation study on the Adv+S dataset evidencing the influence of our learning choices on the resolution of both identification and prediction tasks (see appendix G for detailed results). “Joint” rows of Table 3 indicate that the learning of hu and hk is done simultaneously. As shown in Table 3, the sole optimization of d(h, f) fails at estimating physically sounded U . This evidences the ill-posedness in such unconstrained optimization. Table 3 indicates that all introduced regularizations improve the recovery of U w.r.t. the «Only d(h, f)» baseline, while adding d(hu, 0) significantly improves both prediction performances and velocity fields estimation. We highlight that the alternate optimization performs better compared to optimiz-
ing jointly all parameters of hk and hu. Notably, our proposition to optimize hk and hu alternately beats all baselines on both T prediction and U identification (Table 3, Joint rows). Finally, jointly trained models fail at estimating U in Table 3, forcing hu to capture the whole dynamics.

heading: 5 RELATED WORK
text: Grey-box or hybrid modeling, combining ODE/PDE and data based models, has received an increasing focus in the machine learning community (Rico-Martinez et al., 1994; Thompson & Kramer, 1994; Raissi et al., 2020b). Hybrid approaches allow for alleviated computational costs for fluid simulation (Tompson et al., 2017; De Avila Belbute-Peres et al., 2020; Wandel et al., 2021), and show better prediction performances through data specific constraints that preserve physics (Raissi et al., 2020a; Jia et al., 2019). They offer increased interpretability via constraints on convolutional filters (Long et al., 2018; 2019) or on learned residual (Geneva & Zabaras, 2020). Physical knowledge, introduced through ODE/PDE regularization (Psichogios & Ungar, 1992; Bongard & Lipson, 2007; de Bézenac et al., 2018) or Hamiltonian priors (Greydanus et al., 2019; Lee et al., 2021), increases generalization power w.r.t pure ML approaches. Closer to our work, (Mehta et al., 2020; San & Maulik, 2018; Young et al., 2017; Saha et al., 2020) study the learning of a physical model augmented with a statistical component. Yin et al. (2021) tackle the same task, ensuring the uniqueness in the decomposition by constraining the norm of the ML component. We generalize latter approaches and address the well-posedness in the learning of hybrid ML/MB models through additional regularization on the estimated parameters of the physical part. Indeed, to describe natural phenomena relying on hybrid MB/ML models, one major task lies in the estimation of the MB part parameters. This can be done using neural networks (Raissi et al., 2019; Mehta et al., 2020). However, identification tasks being intrinsically ill-posed (Sabatier, 2000), imposing prior knowledge or regularization is necessary to ensure sound estimations (Stewart & Ermon, 2017). Yet, using only prediction as supervision, the recovered parameters are not physically interpretable (de Bézenac et al., 2018; Ayed et al., 2020). To ensure uniqueness of the estimation solution, Ardizzone et al. (2018) use invertible neural networks. Linial et al. (2021); Tait & Damoulas (2020); Saemundsson et al. (2020) combine variational encoding (Kingma & Welling, 2013) and a PDE model, sampling the space of initial conditions and parameters to solve both identification and prediction. However, such methods only deal with low-dimensional dynamics. Besides low dimensional systems, we also show the soundness of our approach on complex high dimensional and partially observed dynamics.

heading: 6 DISCUSSION
text: We propose in this work an algorithm to learn hybrid MB/ML models. For interpretability purposes, we impose constraints flowing from an upper bound of the prediction error and derive a learning algorithm in a general setting. We prove its well posedness and its convergence in a linear approximation setting. Besides theoretical considerations, we empirically evidence the soundness of our approach thanks to ablation studies and comparison with recent baselines on several low and high dimensional datasets. This work can see several extensions, e.g. considering non uniform 3-D grid for climate models, further considerations on the investigated upper bounds, or less restrictive decomposition hypothesis.

heading: A DISTANCE

heading: A.1 DISTANCE BETWEEN DYNAMICS
text: We here give the definition of the distance d. Let u and v be two functions of L2(Rp,Rp). We consider the distance:
d(u, v) = EX∼pX‖u(X)− v(X)‖2 (18)
Naturally, eq. (18) verifies the triangle inequality, the symmetry and the positiveness. Moreover, in this case, for all functions f , d(., f) is convex. Indeed, for u, v two functions, and λ ∈ [0, 1]:
d(λu+ (1− λ)v, f) = EX∼pX‖λu(X) + (1− λ)v(X)− f(X)‖2 = EX∼pX‖λu(X)− λf(X)− (1− λ)f(X) + (1− λ)v(X)‖2 ≤ λEX∼pX‖u(X)− f(X)‖2 + (1− λ)EX∼pX‖v(X)− f(X)‖2
Hence the convexity of d(., f). This consideration suffices to ensure the convexity of Sk and Su defined in section 3.

heading: A.2 DISTANCE BETWEEN FLOWS
text: Consider the ODE with X(t), X0 ∈ Rp: dX(t)
dt = f(X(t)), X(t = 0) = X0 (19)
Equation (19) admits a unique solution as soon as f is Lipschitz. We note X? this solution. Then, we can defined the flow φf of such ODE as :
[0,T ]×Rp → Rp
t, X0 → φf (t,X0) = X?(t) (20)
With the definition of eq. (20), we can define the distance between two flows of ODE as:
dφ(φu, φf ) = EX0∼pX0 ∫ τ t0 ‖φu(t,X0)− φf (t,X0)‖ dt (21)
dφ is positive and symmetric. Let φu, φv be two flows, we have the triangle inequality:
dφ(φu, φf ) = EX0∼pX0 ∫ τ t0 ‖φu(t,X0)− φf (t,X0)‖ dt
= EX0∼pX0 ∫ τ t0 ‖φu(t,X0)− φv(t,X0) + φv(t,X0) + φf (t,X0)‖ dt
≤ EX0∼pX0 ∫ τ t0 ‖φu(t,X0)− φv(t,X0)‖+ ‖φv(t,X0) + φf (t,X0)‖ dt ≤ dφ(φv, φv) + dφ(φv, φf )
Let φf be fixed, we also have the convexity of dφ(., φf ) with respect to the first argument. Indeed for λ ∈ [0, 1]:
dφ(λφu + (1− λ)φv, f) = EX0∼pX0 ∫ τ t0 ‖λφu(t,X0) + (1− λ)φv − φf (t,X0)‖ dt
= EX0∼pX0 ∫ τ t0 ‖λφu(t,X0) + (1− λ)φv − λφf (t,X0)− (1− λ)φf (t,X0)‖ dt ≤ λ dφ(φu, φv) + (1− λ)dφ(φv, φf )
However, in this case the convexity is not ensured with respect to u and vThis is the reason why for theoretical investigations, we consider the distance d instead of dφ.
Nonetheless, dφ(φu, φf ) = 0 =⇒ φu = φf =⇒ u = f .

heading: B REMARK ON ADDITIVE DECOMPOSITION
text: First, note that in the case of a metric space the decomposition as defined in eq. (2) always exists.
We now detail an intuition for the well-posedness of such decomposition.
Let Hk be a closed convex subset of functions of an Hilbert space (E,<,>), and f the function we want to approximate with partial knowledge (represented by the space of hypothesisHk). Then, thanks to Hilbert projection lemma, we have the uniqueness of the minimizer of ming∈Hk ‖f − g‖, i.e it exists one unique hk ∈ Hk such that: ∀g ∈ Hk, ‖f − hk‖ ≤ ‖f − g‖. Finally, the additive decomposition hypothesis presents a remarkable advantage in the case of a vector space. Indeed, if Hk is a (closed) vector space, let H⊥k be its supplementary in E, then we have the uniqueness in the decomposition: f = fHk + fH⊥k , where fH⊥k ∈ H ⊥ k and fHk ∈ Hk.
The existence and uniqueness flowing directly from the additive decomposition hypothesis, this can explain why such assumption is common when bridging ML and MB hypothesis.

heading: C UPPER BOUNDS

heading: C.1 DERIVATION OF EQUATION (4)
text: The first upper bound is a simple use of the triangle inequality:
d(h, f) = d(h, f) + d(hk, f)− d(hk, f) ≤ d(hk, f) + |d(h, f)− d(hk, f)| ≤ d(hk, f) + d(h, hk)

heading: C.2 DERIVATION OF EQUATION (5)
text: To derive the second upper bound, we assume that fprk comes from an overall dynamics f pr obeying the additive decomposition hypothesis of eq. (2) so that fpr and fprk verifies: f pr = fprk +f pr u . First, with computations similar to eq. (4), we have:
d(h, f) ≤ d(h, fpr) + d(fpr, f) (22)
Then:
d(h, fpr) = d(h, fpr) + d(hk, f pr k )− d(hk, f pr k )
≤ d(hk, fprk ) + |d(h, f pr)− d(hk, fprk )| ≤ d(hk, fprk ) + |d(h, f pr)− d(h, fprk )− d(h, f pr k ) + d(hk, f pr k )| ≤ d(hk, fprk ) + |d(h, f pr)− d(h, fprk )|+ |d(hk, f pr k )− d(h, hk)| ≤ d(hk, fprk ) + d(f pr, fprk ) + d(h, hk) (23)
Combining Equations (22) and (23), we retrieve eq. (5):
d(h, f) ≤ d(hk, fprk ) + d(h, hk) + d(f pr, fprk ) + d(f pr, f) (24)
and we have: Γ = d(fpr, fprk ) + d(f pr, f). Γ is a constant of the problem that cannot be optimized.
C.3 UPPER-BOUND USING AUXILIARY DYNAMICS fpr
Let fpr be the dynamics of model data, we can link up the error made by h on true data (following dynamics f ) and the error made by h on model data (with dynamics fpr) via:
d(h, f) ≤ d(h, fpr) + d(fpr, f) (25)
Thus a pre-training on auxiliary data of dynamics fpr amounts to control the term d(h, fpr) in the upper-bound of eq. (25).

heading: C.4 SELF-SUPERVISION
text: Let h = hk + hu be the function to learn and Gψ the recognition network providing an estimate θ̂ik of the parameters from an initial sequence (Xit0 , . . . , X i t0+k∆t
). This learning setting corresponds to how velocity fields are learned from consecutive measurements of the tracer fields T in section 4.2.
To compute d(hk, f pr k ) in the case where f pr = h?, where h? = h?k+h ? u is a learned model, we rely on the computed θk associated to h?k (thanks for example to the algorithm of section 3.2 associated to eq. (4)) to generate a synthetic dataset with achievable supervision in the space of the parameters θk.
From a real initial sequence (Xit0 , . . . , X i t0+k∆t ), we can estimate the unknown parameter θik associated to sequence i with the recognition network G ? ψ learned with h
?, i.e θik = G ? ψ ( Xit0 , . . . , X i t0+k∆t ) . Then, integrating from the initial condition Xit0 , we generate a trajectory of known parameters θik with dynamics h ? denoted by: X̃i = ( X̃it0 , . . . , X̃ i tn ) . Sampling
the space of initial conditions, we obtain a synthetic dataset: ( (X̃1, θ1k), . . . , (X̃ m, θmk ) )
enabling us to perform self-supervision for Gψ . Let θ̂ik be the parameters estimated by Gψ from the simulated( X̃it , . . . , X̃ i t+k∆t ) , we make the following approximation:
d(hk, fk) ≈ 1
m m∑ i=1 ∥∥∥θ̂ik − θik∥∥∥ 2
(26)

heading: D PROOFS
text: D.1 NOTE ON THE CONVEXITY OF Sk AND Su

heading: Convexity of Sk
text: Proof. Let u, v ∈ Sk: d(tu+ (1− t)v, f) = ‖tu+ (1− t)v − f‖ = ‖tu− tf + (1− t)v − (1− t)f‖
≤ tµ1 + (1− t)µ1 = µ1 Hence the convexity of Sk.

heading: Convexity of Su
text: Proof. Let t ∈ [0, 1] and u, v ∈ Su. d(hk, hk + tu+ (1− t)v) = d(0, tu+ (1− t)v)
≤ td(u, 0) + (1− t)d(v, 0) ≤ µ2
Hence the convexity of Su.

heading: D.2 ODE IDENTIFICATION
text: Consider the following set: SA = {X(t) ∈ C1([0, T ],Rp) such that: ∃A ∈ Mp,p(R), X ′ = AX}, where T > 0.
SA is not a convex set. Consider u and v in SA, and consider Au and Av so that u′(t) = Auu(t) and v′(t) = Avv(t). For λ ∈]0, 1[: we have:
[λu+ (1− λ)v]′ = λu′ + (1− λ)v′
= λAuu+ (1− λ)Avv In general the last term is not equal to Aλu+(1−λ)v(λu + (1 − λ)v), for some matrix Aλu+(1−λ)v . Thus SA is not a convex set. However, discretizing the trajectories and employing a simple integration scheme leads to considering the following cost function:
L(A) = ∑ t ‖ ( Xs(t+ 1)− (A∆t+ Id)XA(t) ) ‖22, (27)
As a least square regression problem, L(A) is convex with respect to A. A least square regression setting can also be recovered using more complex integration schemes, or several time steps integration.

heading: D.3 PROOF FOR WELL-POSEDNESS OF EQUATION (7)
text: We set ourselves in the Hilbert space of squared integrable functions with the canonical scalar product ( L2(Rp,Rp), <,> ) . For further consideration on such functional space we refer to (Droniou, 2001).
We assume thatHk hence Sk is convex and a relatively compact family of functions.
Convexity of Sk + Su Let S = Sk + Su = {f |∃fk ∈ Sk, fu ∈ Su, f = fk + fu}. Let f, g ∈ S and λ ∈]0, 1[:
λf + (1− λ)g = λfk + (1− λ)gk + λfu + (1− λ)gu ∈ Sk + Su
Hence the convexity of S.
CLOSENESS OF Su We show that Su is a closed set. Indeed, Su = g−1([0, µu]), where g(u) = ‖u‖, Because g is 1-Lipschitz (using the triangle inequality), g is continuous. Therefore Su is closed set as the inverse image of a closed set by continuous function.
Sequential Limit We now show that S is a closed set thanks to the sequential characterisation: let fn a converging sequence of elements of S and denote f its limit. We prove that fn converges in S. Because ∀n, fn ∈ S, we have: fn = fnk + fnu , where fnu ∈ Su and fnk ∈ Sk. Thanks to the relative compactness of Sk, we can extract a converging sub-sequence, of indexes nj , from fnk so that f nj k → fk ∈ Sk.
Because fn → f , the sub-sequence fnj converges: fnj → f . By definition, fnju is a sequence of Su and we also have that: f nj u = fnj − fnjk . Because the right member of the equation converges (as a sum of converging functions), the left member of the equation converges i.e. fnju converges.
Since Su is a closed set f nj u converges in Su. We write fu its limit. Therefore, f nj u = fnj − fnjk → f − fk = fu ∈ Su. Hence, f = fu + fk with fu ∈ Su and fk ∈ Sk. Therefore S is a closed set. Finally, we can apply Hilbert projection lemma on the closed convex set S and retrieve the uniqueness of the minimizer of eq. (7).
Remark The relative compactness of a family of functions is a common assumption in functional analysis. For example, in the study of differential equation Cauchy-Peano theorem provides the existence to the solution of an ODE under the assumption of relative compactness.
Also, Ascoli theorem provides the relative compactness of a family of function F under the hypothesis of the equi-continuity of F and the relative compactness of the image space A(x) = {f(x)|f ∈ F}.

heading: D.4 PROOF OF PROPOSITION 2
text: We now set ourselves in the Hilbert space ( L2([0, T ],Rp), <,> ) of squared integrable functions,
where <,> is the canonical scalar product of L2([0, T ],Rp).
Proof. Let A be a given invertible matrix. We consider the following space SD = {X ∈ C1([0, T ],Rp) such that: ∃D ∈ Rp, X ′ = AX + D and X(t = 0) = X0}, where T > 0. We show that SD is a closed convex set.
Convexity Indeed, let λ ∈]0, 1[ and u, v ∈ SD. λu+ (1− λ)v is differentiable and: [λu+ (1− λ)v)]′ = λu′ + (1− λ)v′ = A(λu+ (1− λ)v) +D,
Where D = λDu + (1− λ)Dv . Hence λu+ (1− λ)v ∈ SD.
Closeness via Affine-Space To prove the closeness of SD, we prove that it is an affine space of finite dimension.
Let g the application that to any vector D ∈ Rd associate the solution XD. Let D0 ∈ RD, we show that gD0 : D → g(D0 +D)− g(D0) is a linear application. Naturally, for gD0(0Rp) = 0L2 . Then for D 6= 0Rp we have:
gD0(D) = e At(X0 +A −1(D0 +D))−A−1(D0 +D)− eAt(X0 +A−1(D0) +A−1D0 = eAtA−1D
Therefore gD0 is a linear function and g is an affine function.
Moreover, g is an injection. Indeed, if two functions are equals, then they have at most one inverse image by g thanks to Cauchy-Lipschitz theorem. Therefore it defines a bijection of Rd in g(Rd). Since, SD = g(Rd), SD is an affine space of dimension p and g is continuous in particular for the canonical norm induced on L2([0, T ],Rp). Therefore SD is an affine space of finite dimension and is a closed set.
Finding a Unique Minimizer We conclude by applying Hilbert projection lemma: our problem of minimizing ∫ T 0 ∥∥Xs(τ)−XD(τ)∥∥, amounts to an orthogonal projection problem. Because SD is a closed convex set, we have existence and uniqueness of such projection. Therefore, it exists a unique function XD ∈ SD and a unique vector D minimizing its distance to the function Xs.

heading: D.5 ALGORITHM IN LINEAR SETTING
text: We detail in Algorithm 2 the alternate projection algorithm in a linear setting. We denote Y = (Xit0+∆t, X i t0+n∆t ) and X = (Xit0 , X i t0+(n−1).∆t). For readability purposes we set ∆t = 1.
Algorithm 2 Alternate estimation: Linear Setting

heading: Result: A ∈Mp,p(R), D ∈ Rp
text: k = 0, D0 = 0, A−10 = 0 A 0 0 = minA‖Y −XA‖
while ‖Dk −Dk−1‖ > and ‖Ak −Ak−1‖ > do Dk+1 = minD ‖Y −XAk −D‖22 + λ‖D‖22 Ak+1 = minA ‖Y −XA−Dk+1‖22 + γ‖Y −XA‖22 k ← k + 1
end

heading: D.6 PROOF TO PROPOSITION 3
text: Naturally, one could estimate jointly D and A using least square regression. However, the idea is to verify the convergence of such alternate algorithm in a simple case. We conduct the proof for the first dimension of Y to lighten notations, meaning that we are regressing the first dimension of Y against the X .
A similar reasoning for the other dimension completes the proof.
Proof. We first give the analytical solution for D. Let An be fixed.
Estimation of D Consider:
LD = ‖Y −XAn −D‖22 + λ‖D‖ 2 2 (28)
where D = (d, . . . , d) ∈ RQ. For Q samples, we find d so that ∂L∂d = 0:
∂L ∂d = 0⇔ −2 ∗ Q∑ i=1 (yi −XiAn − d) + 2λd = 0
⇔ Qd+ λd = Q∑ i=1 (yi −XiAn)
⇔ d(Q+ λ) = Q∑ i=1 (yi −XiAn)
⇔ d = Y −XA 1 + λ/Q
where Y −XA = 1Q ∑Q i=1(yi −XiAn).
Estimation of A Let D be fixed and consider:
LA = ‖Y −XA−D‖22 + γ‖Y −XA‖ 2 2 (29)
Similarly, we aim to cancel the first derivative of LA with respect to all parameters of A = (a1, .., ap):
∂LA ∂aj = 0⇔− 2 ∗ Q∑ i=1 xi,j(yi − a0xi,0 + · · ·+ apxi,p − d)
− 2γ ∗ Q∑ i=1 xi,j(yi − a0xi,0 + · · ·+ apxi,p) = 0
⇔− 2Xt(Y −XA−D)− 2γXt(Y −XA) = 0 ⇔(1 + γ)XtXA−Xt(Y −D)− γXtY = 0 ⇔(1 + γ)XtXA = Xt ( γY + (Y −D)
) ⇔A = B −1Xt
1 + γ
( (1 + γ)Y −D ) (30)
whereB = XtX . Equation (30) indicates that as soon aD converges, An converges. Thus, we now prove the convergence of (Dn). Then, for n > 1 consider:∥∥Dn+1 −Dn∥∥ = 1
1 + λ/Q ∥∥∥Y −XAn − Y −XAn−1∥∥∥ = 1
1 + λ/Q ∥∥∥X(An −An−1)∥∥∥ = 1
(1 + λ/Q)(1 + γ) ∥∥∥XB−1Xt([(1 + γ)Y −Dn]− [(1 + γ)Y −Dn−1)]∥∥∥ = 1
(1 + λ/Q)(1 + γ) ∥∥∥XB−1Xt[Dn−1 −Dn]∥∥∥ ≤ K
(1 + λ/Q)(1 + γ) ∥∥Dn−1 −Dn∥∥ where K = ‖XB−1Xt‖.
Therefore, for λ, γ, sufficiently large, K(1+λ/Q)(1+γ) < 1. ‖Dn −Dn−1‖ converges as a positive decreasing sequence. Finally, the sequence of (Dn) converge and so the sequence of (An).
In conclusion, the proposed algorithm converges.

heading: E DATASETS
text: In this section, we provide exhaustive simulation details for the damped pendulum, Lotka-Volterra, and both geophysical datasets.

heading: E.1 DAMPED-PENDULUM
text: For the damped pendulum data, eq. (15) is integrated with ∆t = 0.2s using a Runge-Kutta 4- 5 scheme from t = 0 up to t = 10s. Both the pulsation ω0 and the damping coefficient k are fixed across the dataset. We generate 100/50/50 sequences respectively for train, validation and test sampling over the initial conditions so that (θ, θ̇) ∼ U ( [−π/2, π/2]× [−0.1, 0.1] ) .
Small Oscillations To linearize the pendulum, we consider the small oscillations regime and take the initial conditions so that : (θ, θ̇) ∼ U ( [−π/6, π/6] × [−0.1, 0.1] ) . In that case eq. (15) writes as: d
dt ( θ̇ θ ) = ( −λ gL 1 0 )( θ̇ θ ) (31)
and following notations of section 3.3, we have: DA = 0 and A = ( −λ gL 1 0 )

heading: E.2 LOTKA-VOLTERRA
text: For Lotka-Volterra data, eq. (16) is integrated with ∆t = 0.05 using a Runge-Kutta 4-5 scheme from t = 0 up to t = 20. All parameters α, β, γ, δ are set to 1 across the dataset. We generate 100/50/50 sequences respectively for train, validation and test sampling over the initial prey and predators populations so that (x, y) ∼ U ( [0, 2]2 ) .
Practical Issues and Adaptation Assuming that α and γ have positive values makes the following problem arises: the trajectories defined by hk for the prey are unbounded, whereas the trajectories defined by eq. (16) are. Minimizing d(hk, f) over long term horizon will lead in an underestimation of α to match the bounded behaviour of true data. Therefore, we enforce d(hk, f) on the prey component as soon as the number of predator is small. In practice, we set this threshold to 0.15.

heading: E.3 GEOPHYSICAL DATASETS
text: We present in this section introductory tools for the understanding of the fluid dynamics data presented in section 4.2. We first introduce the physical modeling of ocean dynamics. Then, we outline the Adv+S dataset simulation which draws from ocean modeling. Finally, we introduce the Natl dataset and the proxy data used in the experiments.
Introduction To Ocean Modeling The increase in ocean observations thanks to satellites and floats enabled a great development in Earth modeling over the last decades. The ocean circulation, that is the current velocity fields dynamics, are now realistically modeled in tri-dimensional structured models such as NEMO (Madec, 2008).
Such models rely on in-depth physical knowledge of the studied system and its representation through partial differential equations. Integrated over depth, the equations associated to the transport of the Sea Surface Temperature T by a time-varying horizontal velocity field U can be written as:
∂T ∂t = −∇.(TU) +DT + FT (32)
∂U
∂t = −(U.∇)U + f ∧ U − g′∇h+DU + FU (33)
where f is the Coriolis parameter, h the depth of the surface layer obtained from sea surface height (SSH) observations, g′ the reduced gravity which takes the stratification in density of the ocean into account such that g′ ≈ g.10−3. In a two-dimensional setting, ∇(TU) refers to the advection of a
scalar quantity T by a velocity field U = (u, v) and writes as : ∇(TU) = ∂T∂x u+ ∂T ∂y v. The mixing terms, referred to as DT/U and the forcings FT/U , are not known.
In the context of the presented work, the physical state Zt = (Tt, Ut), fX and fY from eq. (1) can be interpreted as follows: fX represents the dynamics of the observed T , i.e. fX(T ) = −∇.(TU)+ DT + FT in eq. (32). fY represents the dynamics of U in eq. (33), i.e. fY (U, h) = −(U.∇)U + f ∧ U − g′∇h+DU + FU . Whereas T is observed by satellites, U is not known. However, the Sea Surface Height (SSH) could be used to compute coarse estimates of U . Indeed, under hypothesis such as stationarity (∂U∂t = 0), incompressibility ((U.∇)U = 0)), forcings can be omited. In this case, eq. (33) can be rewritten into f ∧ U = −g′∇h (34) When projected onto x and y axis, eq. (34) becomes
−fv = −g′ ∂h ∂x , fu = −g′ ∂h ∂y , (35)
Note that eq. (34) and eq. (35) do not hold at fine scales as the stationarity and incompressibility assumptions only hold at large scale. In this case, the SSH h can be regarded as a stream function.
Both datasets considered in the paper follow the same equations approximating the tracer equation (eq. (17)) inspired by eq. (32):
∂T ∂t = −∇.(TU) + S (36)
We study the equations 32 and 33 in an incremental approach. In the following parts, we describe how T , U and S are computed in both datasets Adv+S and Natl.
E.3.1 ADV+S
We first investigate a dataset generated following simplifying assumptions (Adv+S). We don’t rely on true U and S, we instead build them so that they correspond to our hypothesis.
Building a Velocity Field U Under stationarity and incompressibility hypothesis, U can be approximated from a stream function H. Note that, in this dataset, H is not equal to the SSH h, it is simulated following (Boffetta et al., 2001):
H(x, y, t) = − tanh[ y −B(t)× cos kx√ 1 + k2B(t)2 × sin2kx ] + cy, (37)
As introduced precedently (see eq. (34)), eq. (33) can be simplified and we compute U = (u, v) so that it follows:
u = −∂H ∂y , v = ∂H ∂x
(38)
Note that B varies periodically with time according to B = B0 + cos(ωt + φ). We compute 10 different velocity fields sampling random parameters B0, k, c, ω, , φ.
Building a Source Term S In eq. (32), the diffusion term DT is omitted. We generate the source term S so that it represents the forcing term FT in eq. (36). To illustrate heat exchanges, we draw from Frankignoul (1985). This source term is a non linear transformation of U = (u, v) multiplied by the difference between the ocean temperature and a reference temperature:
S(U, T ) = we × (T − Te) where we = { 0 if ∂H∂t < 10 −4
1 otherwise. where Te is the sequence mean image (computed without source).
Dataset Generation Using computed U and S, we integrate eq. (36) with ∆t = 8640s over 30 days, using a Semi-Lagrangian scheme (see explanations below). We generate 800/100/200 sequences respectively for train, validation and test sampling over the initial conditions, which are images of size 64 × 64 sampled from Natl dataset. Finally, for integration, we impose East-West periodic conditions, implying that what comes out the left part re-enters at the right, and reciprocally. We also impose velocity to be null on both top and bottom parts of the image.
Semi-Lagrangian Integration Unlike Eulerian scheme, relying on time discretization of the derivative, the semi Lagrangian scheme relies on the constancy of the solution of a PDE along a characteristic curve. Consider a solution to the advection equation, i.e. eq. (36) with S = 0. The method of characteristics consists in exhibiting curves (x(s), t(s)) along which the derivative of the solution T is simple, i.e ∂T∂s (x(s), t(s)) = 0. For a 1D constant advection scheme, computations lead to:
dt ds = 1 =⇒ s = t
dx ds = U =⇒ x = x0 + Ut
giving therefore, T (x, t) = T0(x − Ut), linking the value of the solution at all time to its initial condition. Therefore from a single observation at t0, it suffices to estimate the original departure points x0 − Ut to infer the prediction at t. However, when U is not constant in time, the method remains doable, not along characteristic lines defined by : (x0 + Ut), but along characteristic curves which are given by:
dt ds = 1 =⇒ s = t
dx ds = U(x, t) (39)
A great deal in the semi-Lagrangian literature involves solving correctly eq. (39). We use the conventional mid-point integration rule and the semi-Lagrangian is implemented using Pytorch function gridsample, following in (Jaderberg et al., 2015). Further developments can be found for example in (Diamantakis, 2014).
E.3.2 NATL
This second dataset depicts the actual ocean circulation, i.e. we consider both eq. (32) and eq. (33). In this case, no assumptions are made on U and S represents both diffusion terms DT and forcing terms FT . We access daily data over a year of ocean surface temperature of the North Atlantic observations model resulting from (Ajayi et al., 2019) 1. The dataset covers a 2300km× 2560km zone at 1.5km resolution, in the North Atlantic Ocean.
In this real-life dataset, sea surface height (SSH) partial derivative with respect to x and y serves as proxies to the (unobserved) velocity fields U . Indeed, recall that simplifying hypotheses led us to eq. (35).
We divide the Natl zone into 270 patches of size 64 × 64. For each region, we extract sea surface temperatures, velocity fields, source terms and height variables. We sample 200/20/50 sequences of 1 year, for respectively train, validation and test. In this case, ∆t = 86400s (1 day).

heading: F TRAINING DETAILS
text: All experiments were conducted on NVIDIA TITAN X GPU using Pytorch (Paszke et al., 2019).
Hyper-Parameters Interpretation From eq. (4), two independent terms appear justifying an alternate projections approach.
First, we highlight that strictly minimizing d(hk, f) biases our estimation of hk. However, it may yield a good estimation of hk provided that fk contributes significantly to the prediction of f . Hence, we interpret this loss as an initialization loss. Thus, in most applications, we progressively decrease its magnitude along training as detailed in appendices F.1 to F.3.
On the other hand, d(hu, 0) aims at constraining the free form function hu to make its action as small as possible. We interpret this loss as a stability penalty.
1Details available at : https://meom-group.github.io/swot-natl60/access-data.html
Finally, aiming to recover exact trajectories of observations, we proceed as suggested in (Yin et al., 2021) progressively increasing the hyper-parameters associated to d(h, f).
The practical implementation is summarized in the following algorithm:
Algorithm 3 Alternate estimation: Practical Setting Initialization: θ0u = 0, θ 0 k = minhk∈Hk d(hk, f), λh, λhk , λhu for epoch = 1 : Nepochs do for batch = 1 : Bk do
θn+1k = θ n k − τ1∇θk [λhd(h, f) + λhk`(hk)]
end for batch = Bk : Bu do
θn+1u = θ n u − τ1∇θu [λhd(h, f) + λhud(hu, 0)]
end λh = τ2λh; λhk = 1 τ2 λhk ; λhu = 1 τ2 λhu
end

heading: F.1 DAMPED PENDULUM
text: Architecture Details The physical parameters to be learned is a scalar of dimension 1, and hu is a 1-hidden layer MLP with 200-hidden neurons with leaky-relu activation.
Optimization For this dataset we use RMSProp optimizer with learning rate 0.0004 for 100 epochs with batch size 128. We supervise the trajectories up to t = ∆t × 50, i.e we enforce dφ over (t0 + ∆t, .., t0 + 50∆t). Overall the number of optimization subsequences for training is 17000. We alternate projection on Sk and Su by descending the gradient 10-batches on hk then 10-batches on hu.
Hyperparameters We initialize λhk = 0.1 and decrease it geometrically down to λhk = 0.001. We initialize λh = 0.1 and increase it geometrically up to λh = 100. λhu is fixed through training at 0.1.
The hyper-parameters were chosen by randomly exploring the hyper-parameters space by sampling them so that λ ∼ U(1, 0.1, . . . , 10−4). We select the ones with the lowest prediction errors, i.e with lowest dφ(h, f).
For the ablation study of Table 1, we set to 0 the hyper-parameters associated to the non-considered loss.
The training time for this dataset is 1 hour.

heading: F.2 LOTKA-VOLTERRA
text: Architecture Details The physical parameters to be learned is a vector of dimension 2 accounting for (α, β) in eq. (16), and hu is a 1-hidden layer MLP with 200-hidden neurons with leaky-relu activation.
Optimization We use Adam optimizer with learning rate 0.0005 for 200 epochs with batch size 128. Overall the number of sequences for training is 15000. We supervise the trajectories up to t = ∆t× 25, i.e we enforce dφ over (t0 + ∆t, .., t0 + 25∆t). We alternate projection on Sk and Su by descending the gradient 10-batches on hk then 10-batches on hu.
Hyperparameters We initialize λhk = 0.1 and decrease it geometrically down to λhk = 0.001. We initialize λh = 0.001 and increase it geometrically up to λh = 1. λhu is fixed through training at 0.001.
The hyper-parameters were chosen by randomly exploring the hyper-parameters space by sampling them so that λ ∼ U(1, 0.1, . . . , 10−4). We select the ones with the lowest prediction errors (i.e lowest d(h, f)).
For the ablation study of Table 1, we set to 0 the hyper-parameters associated to the non-considered loss.
The training time for this dataset is 2 hours.

heading: F.3 ADV+S
text: Architectures Details The physical parameters to be estimated are the velocity fieldsU , of dimension (2, 64, 64). As U varies over time, we follow data assimilation principles to map a sequence of 4 consecutive measurements of the tracer field T to the associated velocity field (Gaultier et al., 2013). To do so, we parameterize a recognition network Gψ by U-net with at most 512 latent channels also following the implementation of (Isola et al., 2017), taking as input a sequence of 4 time steps of T : (Tt0 , .., Tt0+3∆t). The residual dynamics hu is learned by a convolutional ResNet, with 1 residual block taking as entry the same sequence of T . We implement hk via a semi-lagrangian scheme, taking as input Tt and the estimated Ut to predict Tt+1.
Optimization We use Adam optimizer with learning rate 0.0001 for 30 epochs with batch size 32. We supervise the trajectories up to t = ∆t× 6, i.e we enforce dφ on (Tt0+∆t, ..., Tt0+6∆t). Overall the number of sequences for training is 36800. We alternate projection on Sk and Su by descending the gradient 4-batches on hk then 6-batches on hu.
Hyperparameters, setting of eq. (4) We initialize λhk = 0.1 and decrease it geometrically down to λhk = 0.00001. We initialize λh = 0.01 and increase it geometrically every epoch up to λh,f = 1000. λhu is fixed through training at 0.1. We select the hyperparameters with the lowest prediction errors (i.e lowest d(h, f)). For the ablation study of Table 1, we set to 0 the hyper-parameters associated to the non-considered loss.
The training time for this dataset is 8 hours.

heading: F.4 NATL
text: Architecture Details The architectures in this setting are identical to the ones described in appendix F.3.
Optimization We use Adam optimizer with learning rate 0.00001 for 50 epochs with batch size 32. Overall the number of sequences for training is 67000. We enforce dφ over 6 time-steps, i.e we supervise the predictions on timesteps: (t0 + ∆t, .., t0 + 6∆t). We use dropout in both Gψ and hu.
Hyperparameters, setting of eq. (4) For this setting, λh geometrically increases from 0.01 up to 100. We initialize λhk = 0.1 and decrease it geometrically down to λhk = 0.00001. λhu is fixed through training at 0.1. We alternate projection on Sk and Su by descending the gradient 10-batches on both hk and hu.
The selected model is the one with lowest prediction errors on validation set (i.e lowest d(h, f)), sampling uniformly the hyperparameters: λ ∼ U(1, 0.1, . . . , 10−4).
Hyperparameters, setting of eq. (5) Because the dynamics of Natl is highly non linear and chaotic, we follow (Jia et al., 2019) and first warm-up the parameters recognition network Gψ on the velocity fields proxies for 10 epochs. For this setting, λh geometrically increase from 0.01 up to 1. λhk is set equal to λh. λhu is fixed through training at 0.01.
After warm-up, we alternate projection on Sk and Su by descending the gradient 100-batches on hk and 300 on hu. In this setting of eq. (5), the selected model is the one with lowest d(h, f) + d(hk, f pr k ) error, sampling uniformly the hyperparameters: λ ∼ U(1, 0.1, . . . , 10−4).
The training time for this dataset is 12 hours.
Baselines We train NODE (Chen et al., 2018) and Aphynity (Yin et al., 2021) on both the Adv+S and Natl dataset. For the training of Aphinity, we set the learning rate at 0.0001 and train on 30 epochs. We initialize λh = 0.01 and increase it geometrically every epoch up to λh = 100. λhu is fixed through training at 0.1. For the training of NODE, we set the learning rate at 0.00004 and train on 50 epochs. To perform prediction, we first encode the 4-consecutive measurements of T (as a 3 × 64 × 64 state) then learn to integrate this state in time thanks to a network h. h is a 3-layer convolutional networks, with 64 hidden channels. It is integrated using RK4 scheme available from https://github.com/rtqichen/torchdiffeq.

heading: G ADDITIONAL RESULTS AND SAMPLES

heading: G.1 RESULTS FOR PENDULUM AND LOTKA-VOLTERRA DATASETS
text: We provide respectively in figs. 4 and 5 phase diagrams for the damped pendulum and Lotka-Volterra experiments. Both graphs in the phase space indicate that the trajectories and their nature are well handled by the learned decomposition, providing a periodic phase space for Lotka-Volterra (fig. 5), and a converging spiral for the damped pendulum (fig. 4).

heading: G.2 RESULTS FOR ADV+S AND NATL
text: In this section, we provide additionial results on both Adv+S and Natl datasets. A thorough ablation study (table 4) gives results with constant hyperparameters λh and λhk (row Vanilla Optim), which validates our hyper-parameters interpretation. Indeed, the results are better when respectively increasing and decreasing λh and λhk . Besides, the row Ours eq. (5) refers to a training performed as introduced in appendix C.4 with fpr = h? trained on eq. (4). Figure 7 shows predictions up to 4 days on the Adv+S data. Finally, figs. 9 and 10 provide results on Natl dataset associated to training relying on both eq. (4) and eq. (5) and with NODE (Chen et al., 2018).

heading: Neural-ODE

Review summaries:
Authors propose a general framework for learning and ensuring identifiability with hybrid models. They derive strong theoretical guarantees as well as a proof of convergence in a simple affine setting. They validate their framework in several simulated experiments and the Natl dataset, and draw comparisons with state-of-the-art hybrid model Aphynity as well as Neural ODEs. Experiments confirm that the proposed framework increase interpretability and maintain high prediction performances

The paper proposes a method for hybrid model-based/machine-leanring learning, wherein a model is decomposed into an interpretable parametric prior and a residual (typically represented by a neural net). The paper demonstrates that prediction error minimization does not accurately identify the parametric component and proposes an alternating optimization method that augments prediction error loss with component-specific losses. The paper also introduces a variant where auxiliary data can be used to define the loss.  Theoretical guarantees are provided where auxiliary losses are given as upper-bound constraints and exact minimization is performed. Also, uniqueness and convergence guarantees are provided for the linear approximation case using a more practical algorithm.  Experiments show the efficacy of the method in linear and non-linear settings, in both identification and prediction.

This paper studies the identification and prediction of dynamical systems from observed data by combining statistical approaches with physical priors. The main idea is to decompose the approximation of unknown dynamics into two parts: a parametric physics-based model $h_k$ (used to incorporate prior knowledge, or to identify parameters), and a non-parametric statistical model $h_u$ (used to fit the unknown, or unobserved degrees of freedom).   Since the statistical part has a priori no restrictions, the overall problem is ill-posed in the sense that the statistical part can fit all of the dynamics. To ensure well-posedness, the authors propose to convert the joint learning problem to learning $h_k$ and $h_u$ subject to constraints: $h_k$ is given an allowed error tolerance and $h_u$ is given an allowed deviation from 0. The authors show that this modification implies well-posedness, and subsequently design an alternative minimization algorithms to solve for $h_k$ and $h_u$ iteratively. It is proved that this algorithm converges in the affine setting (via sequential projection to find intersection of convex sets) and it is demonstrated experimentally in low and high dimensional systems that this method works well in practice outside of the convex setting.

This paper investigate a novel method to recover well-posedness and interpretability for the problem of learning a physical and machine-learning hybrid model, by controlling the ML component and the MB hypothesis. The authors also study the case where auxiliary data are introduced. Experiments show the performance of the proposed method in both synthetic data and real-world ocean dynamics.