Yuan Yin
email: yuan.yin@sorbonne-universite.fr

Vincent Le Guen
email: vincent.le-guen@edf.fr

Jérémie Dona
email: jeremie.dona@sorbonne-universite.fr

Emmanuel De Bézenac
email: emmanuel.de_bezenac@sorbonne-universite.fr

Ibrahim Ayed
email: ibrahim.ayed@sorbonne-universite.fr

Nicolas Thome
email: nicolas.thome@cnam.fr

Patrick Gallinari
email: patrick.gallinari@sorbonne-universite.fr

Augmenting Physical Models with Deep Networks for Complex Dynamics Forecasting

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 nonnegligible 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. reactiondiffusion 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.

Introduction

Modeling and forecasting complex dynamical systems is a major challenge in domains such as environment and climate [START_REF] Rolnick | Tackling climate change with machine learning[END_REF], health science [START_REF] Choi | RETAIN: An interpretable predictive model for healthcare using reverse time attention mechanism[END_REF], and in many industrial applications [START_REF] Toubeau | Deep learning-based multivariate probabilistic forecasting for short-term scheduling in power markets[END_REF]. 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 [START_REF] Shi | Convolutional LSTM network: A machine learning approach for precipitation nowcasting[END_REF][START_REF] Wang | PredRNN++: Towards a resolution of the deep-in-time dilemma in spatiotemporal predictive learning[END_REF][START_REF] Oreshkin | N-BEATS: Neural basis expansion analysis for interpretable time series forecasting[END_REF][START_REF] Donà | Pde-driven spatiotemporal disentanglement[END_REF], and connections have been drawn between deep architectures and numerical ODE solvers, e.g. neural ODEs [START_REF] Chen | Neural ordinary differential equations[END_REF][START_REF] Ayed | Learning dynamical systems from partial observations[END_REF]. 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, [START_REF] Brunton | Discovering governing equations from data by sparse identification of nonlinear dynamical systems[END_REF] and [START_REF] Long | PDE-Net: Learning PDEs from data[END_REF] learn the explicit form of PDEs directly from data, [START_REF] Raissi | Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations[END_REF] and [START_REF] Sirignano | Dgm: A deep learning algorithm for solving partial differential equations[END_REF] use NNs as implicit methods for solving PDEs, Seo, Meng & Liu (2020) learn spatial differences with a graph network, [START_REF] Ummenhofer | Lagrangian fluid simulation with continuous convolutions[END_REF] introduce continuous convolutions for fluid simulations, [START_REF] De Bézenac | Deep learning for physical processes: Incorporating prior scientific knowledge[END_REF] learn the velocity field of an advection-diffusion system, [START_REF] Greydanus | Hamiltonian neural networks[END_REF] and [START_REF] Chen | Symplectic recurrent neural networks[END_REF] 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 [START_REF] Gentine | Could machine learning break the convection parameterization deadlock?[END_REF], computational constraints prescribing a fine grain modeling of the system [START_REF] Ayed | Ep-net: Learning cardiac electrophysiology models for physiology-based constraints in data-driven predictions[END_REF], unknown external factors, forces and sources which are present [START_REF] Large | Diurnal to decadal global forcing for ocean and sea-ice models: The data sets and flux climatologies[END_REF]. 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, Camps-Valls, Stevens, Jung, Denzler, Carvalhais & Prabhat 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 (reactiondiffusion, 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).

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 [START_REF] Pernot | A critical review of statistical calibration/prediction models handling data inconsistency and model inadequacy[END_REF]. Data assimilation techniques, e.g. the Kalman filter [START_REF] Kalman | Stationary wave solutions of a system of reaction-diffusion equations derived from the fitzhugh-nagumo equations[END_REF][START_REF] Becker | Recurrent kalman networks: Factorized inference in high-dimensional deep feature spaces[END_REF], 4D-var [START_REF] Courtier | A strategy for operational implementation of 4d-var, using an incremental approach[END_REF], 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 [START_REF] Chen | Disturbance observer based control for nonlinear systems[END_REF][START_REF] Li | Disturbance observer-based control: methods and applications[END_REF]). 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 [START_REF] Janner | When to trust your model: Model-based policy optimization[END_REF] or by model predictive control [START_REF] Nagabandi | Neural network dynamics for model-based deep reinforcement learning with model-free fine-tuning[END_REF]. 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: [START_REF] Psichogios | A hybrid neural network-first principles approach to process modeling[END_REF][START_REF] Thompson | Modeling chemical processes using prior knowledge and neural networks[END_REF][START_REF] Rico-Martinez | Continuous-time nonlinear signal processing: a neural network based approach for gray box identification[END_REF] 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 [START_REF] Long | PDE-Net: Learning PDEs from data[END_REF][START_REF] De Bézenac | Deep learning for physical processes: Incorporating prior scientific knowledge[END_REF], as discussed in introduction. Recently, some approaches have specifically targetted the MB/ML cooperation. HybridNet [START_REF] Long | Hybridnet: integrating model-based and data-driven learning to predict evolution of dynamical systems[END_REF] and PhICNet [START_REF] Saha | PhICNet: Physics-incorporated convolutional recurrent neural networks for modeling dynamical systems[END_REF] 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. [START_REF] Wang | Integrating model-driven and data-driven methods for power system frequency stability assessment and control[END_REF] and [START_REF] Mehta | Neural dynamical systems[END_REF] 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).

The APHYNITY Model

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

dX t dt = 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 ∈ R d for every t (pendulum equations in Section 4), or X t is a d-dimensional vector field over a spatial domain Ω ⊂ R k , with k ∈ {2, 3}, i.e. X t (x) ∈ R d for every (t, x) ∈ [0, T ] × Ω (reaction-diffusion and wave equations in Section 4). We suppose that we have access to a set of observed trajectories

D = {X • : [0, T ] → A | ∀t ∈ [0, T ], dXt /dt = F (X t )}
, where A is the set of X values (either R d or vector field). In our case, the unknown F has A as domain and we only assume that F ∈ F, with (F, ∥ • ∥) a normed vector space.

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 ∈ F p ⊂ F. The APHYNITY framework leverages the knowledge of F p while mitigating the approximations induced by this simplified model through the combination of physical and data-driven components. F being a vector space, we can write:

F = F p + F a
where F p ∈ F p encodes the incomplete physical knowledge and F a ∈ 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 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 ∈ 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:

min Fp∈Fp,Fa∈F ∥F a ∥ subject to ∀X ∈ D, ∀t, dX t dt = (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 ( 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 F p , and is formulated in the following proposition proven in Appendix B:

Proposition 1 (Existence of a minimizing pair). If F p is a proximinal set ‡, there exists a decomposition minimizing (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 (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 F p is a Chebyshev set ‡, (2) admits a unique minimizer. The F p in this minimizer pair is the metric projection of the unknown F onto 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 F = L 2 , 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.

Solving APHYNITY with deep neural networks

In the following, both terms of the decomposition are parametrized and are denoted as F θp p and F θa a . Solving APHYNITY then consists in estimating the parameters θ p and θ a . θ 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 F: we thus use deep neural networks, which have ‡ A 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.

shown promising performances for the approximation of differential equations [START_REF] Raissi | Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations[END_REF][START_REF] Ayed | Learning dynamical systems from partial observations[END_REF].

When learning the parameters of F θp p and F θa a , we have access to a finite dataset of trajectories discretized with a given temporal resolution ∆t: D train = {(X (i) k∆t ) 0≤k≤⌊ T /∆t⌋ } 1≤i≤N . Solving (2) requires estimating the state derivative dXt /dt appearing in the constraint term. One solution is to approximate this derivative using e.g. finite differences as in [START_REF] Brunton | Discovering governing equations from data by sparse identification of nonlinear dynamical systems[END_REF], [START_REF] Greydanus | Hamiltonian neural networks[END_REF][START_REF] Cranmer | Lagrangian neural networks[END_REF][START_REF] Cranmer | Lagrangian neural networks[END_REF]. 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 (2) using an integral trajectorybased approach: we compute X i k∆t,X 0 from an initial state X (i) 0 using the current F θp p +F θa a dynamics, then enforce the constraint X i k∆t,X 0 = X i k∆t . This leads to our final objective function on (θ p , θ a ):

min θp,θa F θa a subject to ∀i, ∀k, X (i) k∆t = X (i) k∆t (3) where X (i) k∆t is the approximate solution of the integral X (i) 0 + k∆t 0 (F θp p + F θa a )(X s ) ds obtained by a differentiable ODE solver.
In our setting, where we consider situations for which F θ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 θ p than when fitting the simplified physical model F θp p alone (Section 4). With only an incomplete knowledge on the physics, θ 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.

Adaptively constrained optimization

Algorithm 1: APHYNITY Initialization:

λ 0 ≥ 0, τ 1 > 0, τ 2 > 0; for epoch = 1 : N epochs do for iter in 1 : N iter do for batch in 1 : B do θ j+1 = θ j - τ 1 ∇ [λ j L traj (θ j ) + ∥F a ∥] λ j+1 = λ j + τ 2 L traj (θ j+1 )
The formulation in (3) involves constraints which are difficult to enforce exactly in practice. We considered a variant of the method of multipliers [START_REF] Bertsekas | Constrained Optimization and Lagrange Multiplier Methods (Optimization and Neural Computation Series[END_REF] which uses a sequence of Lagrangian relaxations L λ j (θ p , θ a ):

L λ j (θ p , θ a ) = ∥F θa a ∥ + λ j • L traj (θ p , θ a ) (4) where L traj (θ p , θ a ) = N i=1 T /∆t h=1 ∥X (i) h∆t - X (i) h∆t ∥.
This method needs an increasing sequence (λ j ) j such that the successive minima of L λ j converge to a solution (at least a local one) of the constrained problem (3). We select (λ j ) j by using an iterative strategy: starting from a value λ 0 , we iterate, minimizing L λ j by gradient descent §, then update λ j with: λ j+1 = λ j + τ 2 L traj (θ j+1 ), where τ 2 is a chosen hyper-parameter and θ = (θ p , θ 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.

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.

Experimental setting

We describe the three families of equations studied in the experiments. In all experiments, F = L 2 (A) where A is the set of all admissible states for each problem, and the L 2 norm is computed on D train by: ∥F ∥ 2 ≈ i,k ∥F (X (i) k∆t )∥ 2 . All considered sets of physical functionals F p are closed and convex in 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)∥. Implementation details and architectures are given in Appendix E.

Reaction-diffusion equations

We consider a 2D FitzHugh-Nagumo type model [START_REF] Kalman | Stationary wave solutions of a system of reaction-diffusion equations derived from the fitzhugh-nagumo equations[END_REF]. The system is driven by the PDE ∂u ∂t = a∆u + R u (u, v; k), ∂v ∂t = b∆v + R v (u, v) where a and b are respectively the diffusion coefficients of u and v, § Convergence to a local minimum isn't necessary, a few steps are often sufficient for a successful optimization. ∥ This 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.

∆ 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 Ω with periodic boundary conditions. The considered physical models are: • Param PDE (a, b), with unknown (a, b) diffusion terms and without reaction terms:

F p = {F a,b p : (u, v) → (a∆u, b∆v) | a ≥ a min > 0, b ≥ b min > 0};
• Param PDE (a, b, k), the full PDE with unknown parameters:

F p = {F a,b,k p : (u, v) → (a∆u + R u (u, v; k), b∆v + R v (u, v) | a ≥ a min > 0, b ≥ b min > 0, k ≥ k min > 0}.
Damped wave equations We investigate the damped-wave PDE:

∂ 2 w ∂t 2 -c 2 ∆w + k ∂w ∂t = 0
where k is the damping coefficient. The state is X = (w, ∂w ∂t ) and we consider a compact spatial domain Ω 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:

F p = {F c p : (u, v) → (v, c 2 ∆u) | c ∈ [ϵ, +∞) with ϵ > 0}; • Param PDE (c, k): F p = {F c,k p : (u, v) → (v, c 2 ∆u -kv) | c, k ∈ [ϵ, +∞) with ϵ > 0}.
Damped pendulum The evolution follows the ODE d 2 θ /dt 2 + ω 2 0 sin θ + α dθ /dt = 0, where θ(t) is the angle, ω 0 the proper pulsation (T 0 the period) and α the damping coefficient. With state X = (θ, dθ /dt), the ODE is F ω 0 ,α p : X → ( dθ /dt, -ω 2 0 sin θ -α dθ /dt). Our physical models are: • Hamiltonian [START_REF] Greydanus | Hamiltonian neural networks[END_REF]), a conservative approximation, with

F p = {F H p : (u, v) → (∂ y H(u, v), -∂ x H(u, v)) | H ∈ H 1 (R 2 )}, H 1 (R 2
) is the first order Sobolev space. • Param ODE (ω 0 ), the frictionless pendulum:

F p = {F ω 0 ,α=0 p | ω 0 ∈ [ϵ, +∞) with ϵ > 0} • Param ODE (ω 0 , α), the full pendulum equation: F p = {F ω 0 ,α p | ω 0 , α ∈ [ϵ, +∞) with ϵ > 0}.
Baselines As purely data-driven baselines, we use Neural ODE [START_REF] Chen | Neural ordinary differential equations[END_REF] for the three problems and PredRNN++ [START_REF] Wang | PredRNN++: Towards a resolution of the deep-in-time dilemma in spatiotemporal predictive learning[END_REF], 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 [START_REF] Wang | Integrating model-driven and data-driven methods for power system frequency stability assessment and control[END_REF] and [START_REF] Mehta | Neural dynamical systems[END_REF]. 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 [START_REF] Greydanus | Hamiltonian neural networks[END_REF][START_REF] Toth | Hamiltonian generative networks[END_REF] and to the the deep Galerkin method (DGM, [START_REF] Sirignano | Dgm: A deep learning algorithm for solving partial differential equations[END_REF]). See additional details in Appendix E.

Results

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 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 × 10 -3 , b = 5 × 10 -3 , k = 5 × 10 -3 , for (b) c = 330, k = 50 and for (c) T 0 = 6, α = 0.2 as true parameters. log MSEs are computed respectively over 25, 25, and 1. The Deep Galerkin method for the pendulum in complete physics DGM (ω 0 , α), being constrained by the equation, outperforms Neural ODE but is far inferior to APHYNITY models. In the incomplete physics case, DGM (ω 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 (ω 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 1. 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 1. 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 [START_REF] Wang | Integrating model-driven and data-driven methods for power system frequency stability assessment and control[END_REF] and [START_REF] Mehta | Neural dynamical systems[END_REF]. Results shown in Table 1 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 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 1 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 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).

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 dw /dt as the forecast horizon increases: it misses crucial details for the second component dw /dt 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).

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.

Then let us define a function F D in the following way: For a ∈ A, we can find X • ∈ D and t 0 such that X t 0 = a. Differentiating X at t 0 , which is possible by definition of D, we take:

F D (a) = dX t dt t=t 0
For any (F p , F a ) satisfying the constraint in (B.1), we then have that (F p + F a )(a) = dXt /dt |t 0 = F D (a) for all a ∈ A. Conversely, any pair such that (F p , F a ) ∈ F p × F and F p + F a = F D , verifies the constraint.

Thus we have the equivalence between (B.1) and the metric projection formulated as:

min Fp∈Fp F D -F p (B.2)
If F p is proximinal, the projection problem admits a solution which we denote

F ⋆ p . Taking F ⋆ a = F D -F ⋆ p , we have that F ⋆ p + F ⋆ a = F D so that (F ⋆ p , F ⋆ a )
verifies the constraint of (2). Moreover, if there is (F p , F a ) satisfying the constraint of (2), we have that F p + F a = F D by what was shown above and

∥F a ∥ = ∥F D -F p ∥ ≥ ∥F D -F ⋆ p ∥ by definition of F ⋆ p . This shows that (F ⋆ p , F ⋆ a ) is minimal. Moreover, if F p is a Chebyshev set, by uniqueness of the projection, if F p ̸ = F ⋆ p then ∥F a ∥ > ∥F ⋆ a ∥.
Thus the minimal pair is unique.

Appendix C. Parameter estimation in incomplete physical models

Classically, when a set F p ⊂ F summarizing 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:

min Fp∈Fp E X∼D L( X X 0 , X) subject to ∀g ∈ I, X g 0 = g and ∀t, d X g t dt = F p ( X g t ) (C.1)
where L is a discrepancy measure between trajectories. Recall that 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 ̸ ∈ F p , the loss is strictly positive at the minimum. This means that reducing the space of functions F p makes us lose in terms of accuracy. ¶

• The obtained function F p might not even be the most meaningful function from F p as it would try to capture phenomena which are not explainable with functions in F p , thus giving the wrong bias to the calculated function. For example, if one is ¶ is true in theory, although not necessarily in practice when F overfits a small dataset.

considering a dampened periodic trajectory where only the period can be learned in 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:

d 2 θ dt 2 + ω 2 0 θ + α dθ dt = 0
We take the same notations as in the article and parametrize the simplified physical models as:

F a p : X → ( dθ dt , -aθ)
where a > 0 corresponds to ω 2 0 . The corresponding solution for an initial state X 0 , which we denote X a , can then written explicitly as:

θ a t = θ 0 cos √ at
Let us consider damped pendulum solutions X written as:

θ t = θ 0 e -t cos t
which corresponds to:

F : X → ( dθ dt , -2(θ + dθ dt ))
It is then easy to see that the estimate of a with the physical model alone can be obtained by minimizing:

T 0 |e -t cos t -cos √ at| 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 → cos t.

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

F a r : X → (0, (a -2)θ -2 dθ dt )
in order to satisfy the fitting constraint. Here a corresponds to 1 + ω 2 0 not to ω 2 0 as in the simplified case. Minimizing its norm, we obtain a = 2 which gives us the desired solution:

θ t = θ 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:

min Fp∈Fp,Fa∈F ∥F a ∥ subject to ∀g ∈ I, X g 0 = g and ∀t, d X g t dt = (F p + F a )( X g t ) ∀X ∈ D, L(X, X X 0 ) = 0 (D.1)
In the continuous setting where the data is available at all times t, this problem is in fact equivalent to the following one:

min Fp∈Fp E X∼D dX t dt -F p (X t ) (D.2)
where the supervision is done directly over derivatives, obtained through finite-difference schemes. This echoes the proof in Appendix B 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 (D.2) is indeed equivalent to (D.1) in the continuous setting, in the practical discrete one, the way error propagates is not anymore: For (D.1) it is controlled over integrated trajectories while for (D.2) 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 (D.1), 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 (D.2) 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 F = F + ϵ be the function estimated from the sampled points with an error ϵ such that ∥ϵ∥ ∞ ≤ α. Denoting X the corresponding trajectory generated by F , we then have, for all X ∈ D:

∀t, d(X -X) t dt = F (X t ) -F ( X t ) -ϵ( 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:

∀t ∈ [0, T ], ∥(X -X) t ∥ ≤ ∥∇F ∥ ∞ t 0 ∥X s -X s ∥ + αt
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 ( ∂u ∂t , ∂v ∂t ). The purely datadriven Neural ODE uses such ConvNet as its F . The detailed architecture is provided in Table E1. The estimated physical parameters θ p in F p are simply a trainable vector (a, b) ∈ R 2 + or (a, b, k) ∈ R 3 + . Optimization hyperparameters We choose to apply the same hyperparameters for all the reaction-diffusion experiments: N iter = 1, λ 0 = 1, τ 1 = 1 × 10 -3 , τ 2 = 1 × 10 3 .

Appendix E.2. Wave equations

The damped wave equation is defined by

∂ 2 w ∂t 2 -c 2 ∆w + k ∂w ∂t = 0
where c is the wave speed and k is the damping coefficient. The state is X = (w, ∂w ∂t ). We consider a compact spatial domain Ω represented as a 64 × 64 grid and discretize the Laplacian operator similarly. ∆ is implemented using a 5 × 5 discrete Laplace operator in simulation whereas in the experiment is a 3 × 3 Laplace operator. Null Neumann boundary condition are imposed for generation.

Dataset δ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 × exp (x-x 0 ) 2 +(y-y 0 ) 2 σ 2 (E.1)
The amplitude C is fixed to 1, and (x 0 , y 0 ) = (32, 32) to make the Gaussian curve centered for all sequences. However, σ is different for each sequence and uniformly sampled in [10,100]. The same δ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, dw dt )(t 0 ) and predicts all states from t 0 + δt up to t 0 + 25δt.

Neural network architectures

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

Optimization hyperparameters We use the same hyperparameters for the experiments:

N iter = 3, λ 0 = 1, τ 1 = 1 × 10 -4 , τ 2 = 1 × 10 2 . Appendix E.3. Damped pendulum
We consider the non-linear damped pendulum problem, governed by the ODE

d 2 θ dt 2 + ω 2 0 sin θ + α dθ dt = 0
where θ(t) is the angle, ω 0 = 2π T 0 is the proper pulsation (T 0 being the period) and α is the damping coefficient. With the state X = (θ, dθ dt ), the ODE can be written as

dXt dt = F (X t ) with F : X → ( dθ dt , -ω 2 0 sin θ -α dθ dt ).
Dataset For each train / validation / test split, we simulate a dataset with 25 trajectories of 40 timesteps (time interval [0, 20], timestep δt = 0.5) with fixed ODE coefficients (T 0 = 12, α = 0.2) and varying initial conditions. The simulation integrator is Dormand-Prince Runge-Kutta method of order (4)5 (DOPRI5, [START_REF] Dormand | A family of embedded runge-kutta formulae[END_REF]). We also add a small amount of white gaussian noise (σ = 0.01) to the state. Note that our pendulum dataset is much more challenging than the ideal frictionless pendulum considered in [START_REF] Greydanus | Hamiltonian neural networks[END_REF].

Neural network architectures

We detail in Table E2 the neural architectures used for the damped pendulum experiments. All data-driven augmentations for approximating the mapping X t → 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 [START_REF] Greydanus | Hamiltonian neural networks[END_REF][START_REF] Toth | Hamiltonian generative networks[END_REF]) is implemented by a MLP that takes the state X t and outputs a scalar estimation of the Hamiltonian H of the system: the derivative is then computed by an in-graph gradient of H with respect to the input:

F (X t ) = ∂H ∂(dθ/ dt) , -∂H dθ .
Table F1. Ablation study comparing APHYNITY to the vanilla augmentation scheme [START_REF] Wang | Integrating model-driven and data-driven methods for power system frequency stability assessment and control[END_REF][START_REF] Mehta | Neural dynamical systems[END_REF] for the reaction-diffusion equation, wave equation and damped pendulum. 1 × 10 -2 (to avoid possible divergence). For the baseline Neural ODE, we concatenate a and b to each sequence as two channels.

In Table G1, 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 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. to more complex settings. Finally the purely data-driven Neural-ODE fails to cope with the increasing difficulty. To extend the experiments conducted in the paper (section 4) with fixed parameters (T 0 = 6, α = 0.2) and varying initial conditions, we evaluate APHYNITY on a much more challenging dataset where we vary both the parameters (T 0 , α) 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 α) 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 (ω 0 ∈ R + for Param ODE (ω 0 ), (ω 0 , α) ∈ R 2 + for Param ODE (ω 0 , α)). In this varying ODE context, we also compare to the state-of-the-art univariate time series forecasting method N-Beats [START_REF] Oreshkin | N-BEATS: Neural basis expansion analysis for interpretable time series forecasting[END_REF].

Results shown in Table G3 are consistent with those presented in the paper. Pure data-driven models Neural ODE [START_REF] Chen | Neural ordinary differential equations[END_REF]) and N-Beats [START_REF] Oreshkin | N-BEATS: Neural basis expansion analysis for interpretable time series forecasting[END_REF] fail to properly extrapolate the pendulum dynamics. Incomplete physical models (Hamiltonian and ParamODE (ω 0 )) are even worse since they do not account for friction. Augmenting them with APHYNITY significantly and consistently improves forecasting results and parameter identification.

  Figure 1. Predicted dynamics for the damped pendulum vs. ground truth (GT) trajectories d 2 θ /dt 2 + ω 2 0 sin θ + α dθ /dt = 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).

Figure 2 .Figure 3 .

 23 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).

Figure 4 .

 4 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

  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
		Data-	Neural ODE	-3.76±0.02	n/a	n/a
		driven	PredRNN++	-4.60±0.01	n/a	n/a
	(a) Reaction-	Incomplete physics	Param PDE (a, b) APHYNITY Param PDE (a, b)	-1.26±0.02 -5.10±0.21	67.6 2.3	n/a 67
	diffusion
			Param PDE (a, b, k)	-9.34±0.20	0.17	n/a
		Complete	APHYNITY Param PDE (a, b, k) -9.35±0.02	0.096	1.5e-6
		physics	True PDE	-8.81±0.05	n/a	n/a
			APHYNITY True PDE	-9.17±0.02	n/a	1.4e-7
		Data-driven	Neural ODE	-2.51±0.29	n/a	n/a
		Incomplete	Param PDE (c)	0.51±0.07	10.4	n/a
	(b)	physics	APHYNITY Param PDE (c)	-4.64±0.25	0.31	71.
	Wave
	equa-tion	Complete	Param PDE (c, k) APHYNITY Param PDE (c, k)	-4.68±0.55 -6.09±0.28	1.38 0.70	n/a 4.54
		physics	True PDE	-4.66±0.30	n/a	n/a
			APHYNITY True PDE	-5.24±0.45	n/a	0.14
		Data-driven	Neural ODE	-2.84±0.70	n/a	n/a
			Hamiltonian	-0.35±0.10	n/a	n/a
			APHYNITY Hamiltonian	-3.97±1.20	n/a	623
	(c) Damped	Incomplete physics	Param ODE (ω 0 )	-0.14±0.10	13.2	n/a
	pendu-lum		Deep Galerkin Method (ω 0 )	-3.10±0.40	22.1	n/a
			APHYNITY Param ODE (ω 0 )	-7.86±0.60	4.0	132
			Param ODE (ω 0 , α)	-8.28±0.40	0.45	n/a
			Deep Galerkin Method (ω 0 , α)	-3.14±0.40	7.1	n/a
		Complete physics	APHYNITY Param ODE (ω 0 , α) -8.31±0.30	0.39	8.5
			True ODE	-8.58±0.20	n/a	n/a
			APHYNITY True ODE	-8.44±0.20	n/a	2.3
	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

Table E1 .

 E1 ConvNet architecture in reaction-diffusion and wave equation experiments, used as data-driven derivative operator in APHYNITY and Neural ODE[START_REF] Chen | Neural ordinary differential equations[END_REF]).

	Module	Specification
	2D Conv.	3 × 3 kernel, 2 input channels, 16 output channels, 1 pixel zero padding
	2D Batch Norm. No average tracking
	ReLU activation -
	2D Conv.	3 × 3 kernel, 16 input channels, 16 output channels, 1 pixel zero padding
	2D Batch Norm. No average tracking
	ReLU activation -
	2D Conv.	3 × 3 kernel, 16 input channels, 2 output channels, 1 pixel zero padding

Table F2 .

 F2 Detailed ablation study on supervision and optimization for the reactiondiffusion equation, wave equation and damped pendulum.

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

Table G1 .

 G1 Results of the dataset of reaction-diffusion with varying (a, b). k = 5 × 10 -3 is shared across the dataset.Methodlog MSE %Err a %Err b %Err k ∥F a ∥ 2

	Data-driven	Neural ODE (Chen et al. 2018)	-3.61±0.07	n/a	n/a	n/a	n/a
	Incomplete Param PDE (a, b)	-1.32±0.02	55.6	54.1	n/a	n/a
	physics	APHYNITY Param PDE (a, b)	-4.32±0.32 11.8	18.7	n/a	(4.3±0.6)e1
		Param PDE (a, b, k)	-5.54±0.38 1.55	3.10	0.59	n/a
	Complete	APHYNITY Param PDE (a, b, k) -5.72±0.25 1.29	1.23	0.39 (5.9±4.3)e-1
	physics	True PDE	-8.86±0.02 n/a	n/a	n/a	n/a
		APHYNITY True PDE	-8.82±0.15 n/a	n/a	n/a (1.8±0.6)e-5

Table G2 .

 G2 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±0.34	n/a	n/a	n/a
	Incomplete	Param PDE (c)	-1.32±0.27	23.9	n/a	n/a
	physics	APHYNITY Param PDE (c)	-4.51±0.38	3.2	n/a	171
		Param PDE (c, k)	-4.25±0.28	3.54	1.43	n/a
	Complete	APHYNITY Param PDE (c, k) -4.84±0.57	2.41	0.064	3.64
	physics	True PDE (c, k)	-4.51±0.29	n/a	n/a	n/a
		APHYNITY True PDE (c, k)	-4.49±0.22	n/a	n/a	0.0005
	Appendix G.3. Damped pendulum with varying parameters

Acknowledgements

Part of this work has been supported by project DL4CLIM, ANR-19-CHIA-0018-01.

//github.com/

Appendix A. Reminder on proximinal and Chebyshev sets

We begin by giving a definition of proximinal and Chebyshev sets, taken from [START_REF] Fletcher | Chebyshev sets[END_REF]: Definition 1. A proximinal set of a normed space (E, ∥ • ∥) is a subset C ⊂ E such that every x ∈ E admits at least a nearest point in C. Definition 2. A Chebyshev set of a normed space (E, ∥ • ∥) is a subset C ⊂ E such that every x ∈ E admits a unique nearest point in 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 [START_REF] Johnson | A nonconvex set which has the unique nearest point property[END_REF] 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:

is equivalent a metric projection onto F p . If F p is proximinal, (B.1) admits a minimizing pair. If F p is Chebyshev, (B.1) 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 D. By definition, A is the set of points reached by trajectories in D so that:

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

When α tends to 0, we recover the true trajectories X. However, as α is bounded away from 0 by the available temporal resolution, this inequality gives a rough estimate of the way 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 [START_REF] Paszke | Pytorch: An imperative style, highperformance deep learning library[END_REF] and the differentiable ODE solvers with the adjoint method implemented in torchdiffeq. +

Appendix E.1. Reaction-diffusion equations

The system is driven by a FitzHugh-Nagumo type PDE [START_REF] Kalman | Stationary wave solutions of a system of reaction-diffusion equations derived from the fitzhugh-nagumo equations[END_REF])

where a and b are respectively the diffusion coefficients of u and v, ∆ is the Laplace operator. The local reaction terms are

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

Dataset Starting from a randomly sampled initial state X init ∈ [0, 1] 2×32×32 , we generate states by integrating the true PDE with fixed a, b, and k in a dataset (a = 1 × 10 -3 , b = 5 × 10 -3 , k = 5 × 10 -3 ). We firstly simulate high time-resolution (δt sim = 0.001) sequences with explicit finite difference method. We then extract states every δt 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.

+ https://github.com/rtqichen/torchdiffeq Optimization hyperparameters The hyperparameters of the APHYNITY optimization algorithm (N iter, λ 0 , τ 1 , τ 2 ) were cross-validated on the validation set and are shown in Table E3. All models were trained with a maximum number of 5000 steps with early stopping.

Table E3. Hyperparameters of the damped pendulum experiments. We conduct ablation studies to show the effectiveness of APHYNITY's adaptive optimization and trajectory-based learning scheme.

Method

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

In Table F1, we consider the ablation case with the vanilla augmentation scheme found in (Le Guen & Thome 2020, [START_REF] Wang | Integrating model-driven and data-driven methods for power system frequency stability assessment and control[END_REF][START_REF] Mehta | Neural dynamical systems[END_REF]), 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.

Appendix F.2. Detailed ablation study

We conduct also two other ablations in

2), where dXt dt 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 λ shown in Algorithm 1. This case is equivalent to λ = 1, τ 2 = 0.

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.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 ∈ [1 × 10 -3 , 2 × 10 -3 ] and b ∈ [3 × 10 -3 , 7 × 10 -3 ]. k is still fixed to 5 × 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 E1 except that E takes 5 frames (10 channels) as input and E outputs a vector of estimated (ã, b) after applying a sigmoid activation scaled by Appendix 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, ∂w ∂t ). The architecture of the encoder E is the same as in Table E1 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 Appendix 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


Stéphanie Engelen

Aleksandra Alicja Wróblewska

Kevin De Bruycker

Resat Aksakal

Vincent Ladmiral

Sylvain Caillol

Filip E Du Prez

Sustainable Design of Vanillin-Based Vitrimers Using Vinylogous Urethane Chemistry

Research on bio-based covalent adaptable networks is popular nowadays in the search for an optimal implementation of thermoset materials and composites in a circular context. Herein, a vanillin derivative is integrated into vitrimers with promising material properties in which the vinylogous urethane associative chemistry has been used as dynamic covalent chemistry platform. The vanillin derivative, 2-methoxyhydroquinone, is epoxidised and aminated by aqueous ammonia, with the formation of a bi-functional aromatic β-hydroxy-amine. The straightforward synthesis protocol is high yielding and up-scalable, without the need of any chromatographic purification step. The presented rigid, catalyst-free vitrimers have a high renewable carbon content (up to 86 %), glass transition temperatures up to 80 °C, show very fast reprocessing and consequently a swift recyclability with relaxation times in the range of seconds by virtue of the applied β-hydroxyamine functionality. This research thus provides a sustainable approach to vanillin-based vitrimers and fits in the growing interest for the design of recyclable crosslinked polymer materials.

Introduction

The pursuit for more sustainable alternatives for commodity materials and petrochemistry-based molecules is nowadays one of the key research themes because of the ever-increasing environmental concerns, social awareness, as well as the persistently rising prices of fossil feedstocks related to their decreasing availability. 1,2 Likewise, overall greenness implies the prevention of waste and thus indirectly the recyclability of materials, a well-known issue with the class of thermoset materials. 3,4 These crosslinked materials are mechanically robust and chemically stable, yet lack the recyclability, reprocessability and malleability of thermoplastics. [START_REF] Billmeyer | Textbook of Polymer Science[END_REF] The introduction of exchangeable covalent bonds that lead to network topological rearrangements upon a trigger such as temperature or light, enables the recycling of thermoset-like materials, resulting in a sustainable alternative to commodity materials. Covalent adaptable networks (CANs) have attracted extensive interest in the past years. [START_REF] Zheng | [END_REF]7 Emerging from the network chemistry, the crosslink exchange mechanism can follow either an associative or a dissociative path, resulting in different material properties. 8 Whereas associative CANs exhibit a fixed crosslinking density, dissociative CANs loose network integrity upon a trigger. Vitrimers, introduced by the pioneering work of Leibler et al.

9

, are thermally triggered associative CANs that experience an Arrhenius-like viscosity decrease, rather than a sudden viscosity decrease, as seen for dissociative CANs upon a thermal trigger. 10 Among bio-based CANs, many types of renewable building blocks, ranging from terpenes to oils, 11,12 fatty acids, 13 fructose, 14 lipoic acid 15 , ferulic acid, 16 or vanillin [17][18][19][20] have been incorporated into dynamic materials. Very recently, the field of bio-based CANmaterials has been reviewed. 21 With lignin being the second most abundant biopolymer, after cellulose, and vanillin the only bio-based aromatic building block so far produced on an industrial scale from lignin, vanillin has attracted an immense interest. [22][23][24] The aromatic structure of vanillin is one of its attractive features, in contrast to vegetable oils that often do not result in materials with sufficiently high glass transition temperatures (T g ).

25

Besides its use in thermoset resins or thermoplastics, 26 vanillin has also paved its way to CANs. Liu et al., for instance, transformed vanillin into a trifunctional epoxy that could dynamically react with an anhydride, via transesterification in the presence of a catalyst. 27 Yet, so far, mainly the dynamic imine-based covalent bonds have been explored for vanillincontaining materials. The imine moiety is typically formed by a condensation reaction between vanillin, its aldehyde functionality and an amine. For example, Memon et al. combined diglycidyl ether of ethylene glycol and diglycidyl ether of bisphenol A resins with vanillin-derived imine hardeners, resulting in dynamic epoxy materials with a closed loop recycling. 20 In parallel, Zhao et al. combined epoxidised soybean oil with a dynamic vanillin-derived Schiff base. 28 In both cases, vanillin was combined with a di-amine, resulting in a phenolic dynamic hardener, which is able to cure the epoxy compounds.

Next to polyimine-based vitrimers, vinylogous urethane (VU)-based vitrimers are at least equally interesting. VU moieties are easily obtained by combining primary amines with any alcohol-derived acetoacetate. Moreover, by varying different parameters, ultra-fast processing can be targeted, thus making thermal recycling more convenient. Denissen et al. chemically controlled the viscoelastic properties of VU by studying the influence of different catalysts.

Later on, Spiesschaert et al. showed that the introduction of vinylogous urethanes into an epoxy matrix resulted in an accelerating effect of the transamination, as a result of the presence of hydroxyl groups after reaction of an epoxy with a nucleophilic amine. 30 In addition, the presence of these epoxy units resulted in better network properties, as a result of a partial permanent crosslinking. 30 More recently, Taplan et al.

were able to continuously process low-viscosity VU-based vitrimers, by altering the stoichiometry of amines and by the addition of an organic acid catalyst p-toluenesulphonic acid.

31

Conventionally, the synthesis of suitable bio-based building blocks often requires tedious chromatographic purification steps. The resulting process mass intensity (PMI) 32 , a parameter used to assess the sustainability of reactions conjointly introduced in the CHEM21 toolkit, 33,[START_REF]Chemical Manufacturing Methods for the 21st century Pharmaceutical Industries[END_REF] is often high as large amounts of solvents and silica are consumed. Keeping the green chemistry principles in mind 3 and aiming for very fast processable, bio-based and high T g vitrimers, the work reported here will focus on the sustainable design of vanillin-based building blocks, containing a β-hydroxy-amine functionality and their incorporation in dynamic vinylogous urethane-based networks with high bio-based content. In this context, emphasis was put on 2-methoxyhydroquinone (MHQ), a vanillin derivative that can easily be epoxidised with epichlorohydrin and ring-opened with aqueous ammonia.

35,36

The first part of the study focused on a significant optimisation of the synthetic protocols for the MHQ-derivatives, starting from previously reported procedures, 24,[36][37][38] with the aim to obtain relatively pure monomers without the need of tedious purification techniques such as column chromatography. In a second part, the corresponding bis-amines will be reacted, without further purification, with a terpene-based bisacetoacetate, 39 resulting in a fast processable vinylogous urethane curing agent, for which the reactivity is ascribed to the increased nucleophilicity of the amine by the adjacent hydroxyl group. Figure 1 illustrates the bio-mass origin of the bis-amine, referred to as dihydroxyamino propane of methoxyhydroquinone (DHAP-MHQ), and the 1,4bis(hydroxymethyl)cyclohexane bis-acetoacetate (CDM-AA) building blocks as well as their incorporation into a VU-based polymer network. It will be described that the latter step will be realised either by means of a tri-functional amine crosslinker tris (2-aminoethyl)amine (TREN), added in different ratios, or by using the diglycidyl ether of methoxyhydroquinone (DGE-MHQ) as crosslinker. The chemical nature of the used compounds was chosen with the aim to design vitrimer materials with elevated T g 's, helping to restrain molecular motions of polymeric fragments at typical usage temperature.

Results and discussion

Optimised synthesis of DGE-MHQ and DHAP-MHQ

Because of the known beneficial influence of a protic environment on the transamination exchange rate in vinylogous urethanes, as investigated in previous studies 30,31 , the starting point of this research was the straightforward design and use of a bio-based β-hydroxy-amine. By ringopening an epoxy group with aqueous ammonia, this protic functional handle can be incorporated within the amine building block, rather than in the matrix. 40,41 However, amination of epoxides with aqueous ammonia often results in oligomeric products, as the resulting primary amine is a more reactive nucleophile than ammonia itself.

Moreover, epoxidation can result in oligomer formation by polyaddition of a diphenol. 25 At the start of this research, we aimed for the significant suppression of these side reactions in order to establish well-defined crosslinked materials based on biobased β-hydroxy-amines. Methoxyhydroquinone can be synthesised by a straightforward and environmentally-friendly Dakin oxidation starting from renewable vanillin and sodium percarbonate. 43 In contrast to other renewable vanillin-based diols, such as vanillyl alcohol, MHQ possesses two hydroxyl groups with similar reactivity. Furthermore, these hydroxyl groups are placed in para position with respect to each other, which prevents the formation of benzodioxane derivatives by oglycidylation, a side reaction observed e.g. with gallic acid.

38

In a first step, epichlorohydrin is used as a reactive solvent to epoxidise MHQ (Figure 2b, step 1). Epichlorohydrin owes its 2c, is formed by ring-opening of epichlorohydrin with a phenolic group in the presence of a phase transfer catalyst, here benzyltriethylammonium chloride. The Cl-intermediate itself can be easily transferred back to the epoxy group by further reacting with 5M NaOH (Figure 2b, step 1). Consequently, when less or almost no oligomers are formed, a chromatographic work-up can be avoided, which is beneficial for the upscaling of the monomer synthesis. The corresponding PMI, which represents the total mass of materials used to produce a specified product, including reactants, reagents, solvents and other materials used in the reaction or the work-up, was calculated to be 34 g/g ( Equation 1, ESI †), whereas a PMI involving the use of silica column chromatography would have resulted in a 2.5 times higher value. The prevention of tedious work-ups is thus in line with the green chemistry principles.

In a next step (Figure 2b, step 2), the amination of DGE-MHQ is carried out for the synthesis of the desired bisfunctional β-hydroxy-amine, DHAP-MHQ. Aqueous ammonia is the reactive solvent of choice, because of its low cost, availability and its bio-based nature. Other amination reactions, using e.g. a Gabriel synthesis, a Curtius reaction, a Hofmann rearrangement or direct catalytic conversion of alcohols with ammonia, 44,[START_REF] Basf | US 20120232293A[END_REF] have been proven to be lowyielding, while requiring harsh conditions and generating a lot of by-products. In previous studies, amination either took place in a microwave, in an autoclave or in a reactor. 36,41,[START_REF] Mora | [END_REF] We opted for the use of a reactor because of its upscaling potential.

The amination of the bis-epoxide can result in oligomeric products by nucleophilic attack of the formed primary amine onto unreacted epoxide moieties (Figure 2d). In order to suppress the formation of those oligomers, it is important to work in a diluted environment. Indeed, it was observed that the volume ratio of ammonia to solvent has a significant impact on the oligomer formation and thus purity of the obtained compound (Table 1). The purest product is found by using a large excess of ammonia (Entry 3 and 4, Table 1 and Fig. S2, ESI †). Also a higher volume fraction of ammonia in water to ethanol is desired from the viewpoint of working in a diluted environment.

Model compound study

With the intention to verify whether a β-hydroxy-amine functionality truly has an accelerating effect, compared to a non-hydroxy amine in the transamination, a low molecular weight (MW) model compound study was conducted. The latter allowed to assess the kinetics of the conjugate transamination in the vinylogous urethane chemistry. To this purpose, N-octyl vinylogous urethane (1) and amino-3phenoxypropan-2-ol (2) were prepared by condensation of ethyl acetoacetate with octylamine and, amination of the glycidyl ether of phenol, respectively. 3a). Two distinct proton signals at 1.48 and 1.53 ppm, corresponding to compounds (1) and (3), respectively, could be easily followed during the reaction (Fig. S3, ESI †). The decrease of the fraction of octyl VU (1) as a function of time and temperature is given in Figure 3b. By using an excess of compound (2), the system was shifted kinetically towards a pseudo-first-order reaction, whereby linearity during the first time points was maintained for 80, 100 and 120 °C, evolving in a shift towards chemical equilibrium at higher conversions. In the previously published model compound study by Denissen et al.

47

, when APP (2) was replaced by benzylamine in the pseudo-first-order reaction, a slower transamination was observed in the same reaction conditions. While the transamination of β-hydroxy-amine APP reaches 40 % at 120 °C after 60 minutes (Figure 3b), only 50 % of benzylamine was transaminated. Thus, as expected, the relatively close proximity of the OH-group to the dynamic centre, results in an increased nucleophilicity of the primary amine and thus slightly enhances the transamination reaction rate.

7

Synthesis and properties of the bio-based dynamic materials

The gained insights of the built-in protic handle on the dynamic exchange in the model compound study helped us to design vinylogous urethane vitrimers with an increased dynamic character. Similar to APP used in the model compound study, DHAP-MHQ was reacted with a readily available bifunctional acetoacetate (CDM-AA) and a crosslinker to prepare the network. The bulk condensation reaction between the CDM-AA, DHAP-MHQ and TREN crosslinker resulted in a range of β-hydroxy-amine vinylogous urethane networks with varying DHAP-MHQ to TREN ratio, further referred to as HA-VU X, with X indicating the molar equivalents of DHAP-MHQ (see Table 2). By altering this ratio, materials with varying mechanical properties, T g -values and exchange dynamics can be obtained. In addition to the HA-VU networks, a fully bio-based epoxy network (HA-VU-EP) was prepared with a renewable carbon content of 86 %, by replacing TREN with DGE-MHQ, which is the synthesised bi-functional epoxy. It should be noted that in this latter case, the primary amine can react twice with the epoxy unit to form the network.

All materials were prepared at 90 °C, homogenised by mixing in a Speedmixer, further cured at 100 °C and post-cured at 120 °C. The prepared dynamic materials had 5 mol % excess of primary amines with respect to the formed vinylogous urethane moieties, so that a sufficient amount of amines could enable the fast network rearrangement. For all networks, the acetoacetate was consumed, as observed by ATR-FTIR with the disappearance of the vibration bands of the acetoacetate at 1720 and 1700 cm -1 and the appearance of the vinylogous urethane ones at 1640 and 1600 cm -1 (Fig. S4, ESI †).

The material properties can be tailored by using different stoichiometric ratios. A more bio-based and faster relaxing material can be obtained by decreasing the amount of trifunctional crosslinker, whereas a higher thermal degradation temperature and storage modulus is observed when the crosslinking density, thus the amount of TREN crosslinker is higher (Table 2).

The thermal properties of the materials are displayed in Table 2 and Fig. S5-6, ESI †. DSC-analysis showed T g 's ranging between 50 and 85 °C, linked to the stoichiometry of the material. The most densely crosslinked HA-VU 0.15 exhibits the highest T g of 85 °C, yet is the least bio-based in the HA-VU series. Although a decreasing trend in T g was expected for HA-VU 0.20, HA-VU 0.25 and HA-VU 0.30 related to the decrease in crosslinking density, HA-VU 0.30 showed a slightly higher T g . This might be attributed to the higher DHAP-MHQ content of this material, compared to the other ones, resulting in other interactions such as π-π stacking. A T g of 50 °C is obtained for the epoxy-based network (HA-VU-EP) with 50 % reversible and 50 % permanent crosslinks. This lower T g value might be associated with homopolymerisations of epoxy molecules at elevated temperatures and the absence of the trifunctional crosslinker TREN. Additionally, relatively broad transitions in the DSC thermograms were observed (Fig. S5,ESI †), which might be attributed to hydrogen bonding interactions, related to the aforementioned β-hydroxy-amine group.

The thermal stability of the materials, on the contrary, is only slightly affected by the amount of DHAP-MHQ. Yet, a slower decrease in onset degradation can be visualised for higher crosslinking densities as expected (Figure S6, ESI †). HA-VU 0.15 has a higher onset degradation with a mass loss of 5 % at 254 °C compared to HA-VU 0.20 (240 °C). For the HA-VU-EP, a different degradation profile was observed. Isothermal analysis of HA-VU 0.25 at 150 °C for 1h (Fig. S7,ESI †) demonstrates the stability under the aforementioned synthesis conditions, envisioning reprocessability (vide infra).

Table 2. Overview of the compositions and different properties of the vinylogous urethane-based materials including their renewable carbon content (RCC), relaxation time, swelling and solubility values, as well as mechanical and thermal properties.

a Characteristic relaxation time obtained from stress relaxation experiments at 160 °C. b Swelling ratio and soluble fraction in THF for 24h at r.t. obtained from four sample measurements. c Apparent plateau storage modulus G' obtained from frequency sweep measurements at 160 °C and 1 rad s -1 . d Crosslinking density obtained using equation 3 ESI †. e Glass transition temperature obtained from second heating (rate 10°C/min). f TGA onset-degradation temperatures for 5% weight-loss. g Char yield derived from remaining mass % at 600 °C.

Next to the thermal degradation, different char yields in connection to their stoichiometry are obtained for the materials (Table 2 and Figure S6, ESI †). For example, for HA-VU-EP with the highest amount of DHAP-MHQ, , a char yield of 25 % is present, which is relevant for applications targeting flame retardancy.

48

The resistance to creep is a recurring problem in the design of dynamic materials. In order to minimize this phenomenon, one strategy is to restrain molecular motion of polymeric fragments at service temperature, e.g. by increasing the T g , and thereby suppressing macroscopic material flow. 9 When the usage temperature is below T g (e.g. 30 °C), the material is showing limited creep as seen by the resulting strain being constant as a function of time (Figure 4). However, the material starts to exhibit more creep behaviour around and above its T g (e.g. 55 °C for HA-VU 0.25).

Swelling experiments of the HA-VU(-EP) vitrimers have been conducted in THF to assess the resistance of the network materials to solvent and the corresponding ratios are reported in Table 2. The highest soluble and swelling fraction in THF are obtained for weaker network structures, yet no sample integrity was lost. Besides a solubility and swelling test in THF, HA-VU 0.25 has also been submerged in aqueous ammonia (25 %) for one week to verify the stability and resistance in an aqueous, yet basic solvent. A similar soluble fraction of 4 wt.% is observed in comparison to the solubility test in THF, yet a lower swelling ratio of 170 wt.% is found. Overall a good stability could be observed, despite the sensitivity of vinylogous urethane bonds to hydrolysis. Moreover, comparable soluble fractions to previously reported vitrimers are obtained. 47 The viscoelastic material properties were analysed by stress-relaxation and frequency sweep experiments and are reported in Table 2 and Figure 5. For the stress relaxation experiments, the relaxation modulus was followed as a function of time between 160 and 130 °C, after applying a constant strain of 1% in the linear viscoelastic region. For all materials, with the exception of HA-VU-EP, which is the partially permanent crosslinked network, full stress relaxation was revealed (Figure 5

Recyclability

Recyclability is a key aspect for vitrimeric materials and CANs. Herein, reprocessability was assessed by breaking the sample into pieces, followed by remoulding for 15 min at 150 °C. It is observed that the sample properties, amongst which the T g and degradation temperatures as well as FTIR bands, were merely unaltered after the 4 th recycling step (Figure 6 In this work, an optimised synthesis of a vanillin-based bifunctional compound with β-hydroxy-amine functionalities was shown, allowing to obtain and upscale relatively pure building blocks without the need of tedious purification techniques such as column chromatography. Next to this biobased bifunctional aromatic amine compound, the applied acetoacetate can be readily obtained from a bio-sourced diol. The synthesised monomers were combined in bulk, with the formation of a bio-based vinylogous urethane network, with a renewable carbon content up to 86 %. Overall, greenness also implies the prevention of waste and thus indirectly, the recyclability of materials, a constitutional issue for thermoset materials that can be overcome by the formation of CANs. We reported fast relaxing vitrimers with good thermal properties and recyclability, without impairing in solvent resistance or mechanical properties. The networks were designed with higher T g 's, to restrain molecular motions of polymeric fragments at room temperature and thus to prevent creep.

Figure 1 .

 1 Figure 1. Schematic overview of the vinylogous urethane dynamic network synthesis. Diglycidyl ether of methoxyhydroquinone (DGE-MHQ) can be obtained from lignocellulosic bio-mass and can be transformed into a bis-functional β-hydroxy amine, dihydroxyamino propane of methoxyhydroquinone (DHAP-MHQ), by reaction with aqueous ammonia. 1,4-Bis(hydroxymethyl)cyclohexane (CDM) can be obtained from terpene bio-mass and can be easily transformed into 1,4bis(hydroxymethyl)cyclohexane bis-acetoacetate (CDM-AA) by reaction with tert-butyl acetoacetate. By combining DHAP-MHQ and CDM-AA with either TREN or DGE-MHQ, fast dynamic materials have been obtained.

29 After mixing compound ( 1 )

 291 with five equivalents of compound (2), deuterated benzene was added to the reaction mixture and the kinetics of formation of N-3-phenoxypropan-

Figure 2 .

 2 Figure 2. Synthesis of DHAP-MHQ from MHQ. a. Vanillin can be obtained from lignocellulosic biomass. b. Synthesis of DHAP-MHQ from MHQ. c. Competitive dimerisation upon reaction with epichlorohydrin. d. Possible undesired dimerisation during the amination process.

2-ol VU ( 3 )

 3 were conducted via 1 H-NMR at different temperatures being 80, 100 and 120 °C (Figure

Figure 3 .

 3 Figure 3. a) Reaction scheme showing the catalyst-free exchange of low-MW model compounds for the kinetic study. b) Disappearance of octyl VU compound (1) as a function of time at different temperatures.

  and Fig. S8-11, ESI †). The most densely crosslinked network HA-VU 0.15 shows the slowest relaxation time in the HA-VU series. Those relaxation times in the HA-VU

Figure 4 .

 4 Figure 4. Creep experiment of HA-VU 0.25 with an applied stress of 2000Pa for 5000 seconds at 30 °C and 50 °C.

ConclusionsFigure 6 .

 6 Figure 6. a) Normalised stress relaxation measurement and b) frequency sweep measurement of HA-VU 0.25.

Figure 5 .

 5 Figure 5. a) Recycling up to 4 times by compression moulding at 150°C for 15 minutes, starting from broken sample pieces. b) FT-IR of the recycled materials.

Table 1 .

 1 Different molar equivalents and solvent combinations tested in a pressure tube/reactor for 1h30 at 85 °C, resulting in different DHAP-MHQ purities. For all entries, 1 eq. of DGE-MHQ is used.

	Entry	Eq. NH3.H2O (25 %)	NH3.H2O (25 %)/EtOH (v/v)	Purity (%) a
	1	50	6/5	70
	2	50	3/4	25
	3	50	5/1	97
	4	60	5/1	97

a The product purity was evidenced using LC-MS (λ = 214 nm).

Acknowledgements

Bernhard De Meyer, Dieter Buyst, Jan Goeman are thanked for the technical support while Prof. Dr. Johan M. Winne, Dr. Christian Taplan, Dr. Yann Spiesschaert, Dr. Nezha Badi, Dr. Tapas Debsharma, Filip Van Lijsebetten and Ir. Laurent Van Durme are thanked for fruitful discussions. S.E. acknowledges the Research Foundation-Flanders (FWO) for her PhD fellowship (Grant 1SD4821N). This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme 101021081 (ERC-AdG-2020, CMaC-project).

Notes and references

Conflicts of interest

There are no conflicts to declare.


Romain David

Dharma Akmon

Nichola Burton

Sarah C Davidson

Steven Mceachern

Kristan Kang

Aleksandra Michalewicz

Priyanka Pillai

Kristal Spreadborough

Frankie Stevens

RDA Sensitive Data IG: A shared understanding of Sensitive Data: an update on current progress and looking to future work

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

RDA Sensitive Data IG: A shared understanding of Sensitive Data: an update on current progress and looking to future work.

More information and contact details are available at https://www.rd-alliance.org/groups/sensitive-data-interest-group The RDA Sensitive Data Interest Group aims to promote the FAIR principles and reproducible research, while drawing attention to the unique challenges and opportunities associated with sensitive data and exploring strategies to mitigate risks.

Sensitive data level: how, where and why?

This poster is presented in tandem with the Sensitive Data IG session at RDA Plenary 18 (which shares the same name). In previous meetings, the Sensitive Data IG has scoped the interest and needs of the RDA community. At this meeting, we focus on specific challenges and opportunities associated with working with Sensitive Data. This poster explores one of these:

• how do we understand sensitive data definitions in different regions and disciplines, and

• how might we develop a shared language around sensitive data. For example, developing an understanding of different community agreed vocabularies and how these relate to each other.

The first step will be to identify together common terminology on sensitive data usable across disciplinaries.

This poster complements the RDA Sensitive Data Interest Group session.

It is developed by the Sensitive Data IG co-chairs and is an opportunity for us to convey our ideas and seek input from the community in different fora.

You can join us!

The future work for sensitive data terminology building consist of three actions.

First, case studies of different types of sensitive data must be identified. These will include, for example, case studies from the humanities, medicine, biodiversity, and military cases.

Second, different classification systems for the level of sensitive data must be shown and compared. This should explore, for example, what constitutes a high level of sensitive data in one classification system compared to another.

Third, examples of how sensitive data is managed across different regions must be detailed. For example, what protocols govern the use of sensitive data in a national context compared to the European context?

By presenting these examples, we aim to:

1. promote wider discussion on how sensitive data is currently defined.

draw out how context can impact whether data is considered sensitive over time

In doing so, we take steps towards developing a common language for discussing sensitive data classifications and levels across contexts/disciplines/localities.

Methodology/poll for identification and comparison of classification systems for the level of sensitive data Methodology/poll to describe how sensitive data are managed across regions

• In your region / place / country (to be defined by participants), what are the relevant governance policies or protocols for sensitive data? • Within the region, does governance differ by discipline or industry?

• By what process are definitions of sensitive data re-assessed or updated?

General objectives of the Sensitive Data Interest Group

The IG objectives will be informed by the community. We take the following objectives as a starting point:

• Develop a shared understanding and refined definition of sensitive data.

• Define various levels of sensitivity for data and develop tools to assess this sensitivity.

• Develop an understanding of how sensitivity relates to openness.

• Identify different consent models.

• Identify types of sensitive data holdings and resources across various domains. • Identify existing data definitions and standards for different types of sensitive data. • Identify challenges in collecting, using, archiving and sharing sensitive data.

• Engage with key stakeholders working in the area of sensitive data management/analytics. • Identify existing solutions for sensitive data collection, analysis, storage and dissemination. • Identify differences in how sensitive data is managed between groups and regions.

Take home message

The Interest Group is focusing on developing an understanding of the terminology around Sensitive Data.

We seek a diversity of perspectives across disciplines and regions.

We warmly welcome new and existing members to join this conversation at RDA18 and beyond.

  POSTER LEAD: Romain David (European Research Infrastructure on Highly Pathogenic Agents [ERINHA], 0000-0003-4073-7456) CO-AUTHORS (in alphabetical order): Dharma Akmon (University of Michigan, 0000-0002-1359-0586); Nichola Burton (Australian Research Data Commons, 0000-0003-4470-4846); Sarah C. Davidson (Max Planck Institute of Animal Behavior, The Ohio State University, 0000-0002-2766-9201); Steven McEachern (Australian Data Archive, 0000-0001-7848-4912); Kristan Kang (Australian Research Data Commons,0000-0002-2057-1033); Aleksandra Michalewicz (University of Melbourne, 0000-0002-7328-2470); Priyanka Pillai (University of Melbourne, 0000-0002-3768-8895); Kristal Spreadborough (University of Melbourne, 0000-0002-7022-3213); Frankie Stevens (Australia's Academic and Research Network, 0000-0002-2556-473X). DOI: 10.5281/zenodo.5663084 Presented at Research Data Alliance 18th Plenary -remotely, November 2021.

•

  In your discipline (to be defined by participants), what are the relevant terms to describe sensitivity? • In your region / place / country (to be defined by participants), what are the official classifications for sensitive data? Methodology/poll for case studies to identify types of sensitive data • Six open ended questions during the session calling for sensitive data examples concerning military topics / health data / biodiversityenvironment risk / personal data / industrial data / other • Sharing specific case studies, including horror / success stories. • Open call for helping to build a shared terminology around sensitive data, drawing on your examples, experiences, or case studies.