Projected Neural Differential Equations for Learning Constrained Dynamics

## Abstract

Abstract Neural differential equations offer a powerful approach for learning dynamics from data.
However, they do not impose known constraints that should be obeyed by the learned model.
It is well-known that enforcing constraints in surrogate models can enhance their generalizability and numerical stability.
In this paper, we introduce projected neural differential equations (PNDEs), a new method for constraining neural differential equations based on projection of the learned vector field to the tangent space of the constraint manifold.
In tests on several challenging examples, including chaotic dynamical systems and state-of-the-art power grid models, PNDEs outperform existing methods while requiring fewer hyperparameters.
The proposed approach demonstrates significant potential for enhancing the modeling of constrained dynamical systems, particularly in complex domains where accuracy and reliability are essential.

## 1 Introduction

Figure: Figure 1: Schematic of projected neural differential equations (PNDEs). The vector field of the unconstrained NDE (red arrow) is projected to the tangent space $T_{u}\mathcal{M}$ of the constraint manifold $\mathcal{M}$. For initial conditions $u_{0}\in\mathcal{M}$, solutions (orange) of the projected vector field (green arrow) remain on the manifold and thereby satisfy the constraints.

Numerical simulation of dynamical systems plays a pivotal role in science and engineering. Across fields, simulations extend the reach of scientific inquiry far beyond the limitations of direct observation and experimentation.
However, as the scale and complexity of the systems being studied increases, the computational cost of traditional numerical methods can become prohibitive.
Deep learning surrogates, enabled by the emergence of efficient architectures and specialized accelerator hardware, offer promising alternatives.
By circumventing the computational limitations of conventional approaches, these approaches promise to further advance the use of simulations for scientific inquiry [Kochkov et al., , Jumper et al., , Pathak et al., , Kovachki et al., , Merchant et al., , Azizzadenesheli et al., ].

In power grid modeling, for example, simulations allow grid operators and researchers to examine the stability of the complex electrical networks that underpin modern life.
However, the transition to renewable energy sources presents unprecedented challenges.
As wind and solar generation increasingly contribute to the energy mix, grid operators must contend with distributed and highly variable production, alongside the loss of system inertia that was previously provided by conventional generators [Breithaupt et al., , Milano et al., ].
On the demand side, driven by electric vehicle adoption and the rapid growth of artificial intelligence, forecasts of peak demand have increased rapidly, further threatening grid stability [North American Electric Reliability Corporation
(2023), ].
Unfortunately, existing methods for modeling grid dynamics are computationally infeasible at the scale and accuracy required for realistically modeling future scenarios dominated by renewable energy sources [Matevosyan et al., ].
This problem is further exacerbated by legal mandates requiring grid operators to perform large-scale probabilistic simulations [European Network of Transmission System Operators for Electricity
(2018), ].
By leveraging large amounts of data to learn and predict grid dynamics, machine learning can enable grid operators to efficiently forecast system states, mitigate faults, and maintain stability in future power grids.

Among existing machine learning approaches, neural ordinary differential equations (NODEs) offer a promising approach for data-driven modeling of dynamical systems, combining the expressiveness of neural networks with the interpretability of differential equations [Chen et al., ].
Building upon NODEs, universal differential equations (UDEs) offer a hybrid framework that integrates neural networks with differentiable physical models [Rackauckas et al., ].
We refer to these approaches collectively as *neural differential equations* (NDEs).
NDEs have demonstrated significant potential across various scientific disciplines, including weather and climate modeling [Kochkov et al., , Yu et al., , Gelbrecht et al., ], robotics [Richards et al., ], epidemiology [Dandekar et al., ], bioengineering [Narayanan et al., ], and chaotic systems [Schötz et al., ].
Despite excellent performance within the training data, NDEs often generalize poorly, for example, to new initial conditions or future times.
Due to the iterative nature in which neural networks are applied within NDEs, inaccuracies and uncertainties are quickly amplified and can push the trajectory far outside the training distribution, where performance degrades.

Inductive biases play an important role in avoiding such issues.
In scientific machine learning, for example, inductive biases can take the form of physical constraints, conservation laws, or other domain-specific knowledge that can enhance model performance and interpretability [Karniadakis et al., ].
Physical constraints, in particular, are often incorporated as additional loss terms, for example, in physics-informed neural networks [Raissi et al., , , ] and physics-informed operator learning [Li et al., , Wang et al., , Goswami et al., , White et al., ]. However, a variety of issues can arise with these soft constraints.
For instance, the resulting composite loss can be expensive to evaluate, due to the derivatives in the physics loss, and
challenging to minimize.
Furthermore, a number of numerical issues are known to arise, such as worse conditioning, conflicting loss terms, inaccuracies in the approximated derivatives, and the existence of spurious
solutions satisfying the constraints on the collocation points of the physics loss [Leiteritz and Pflüger, , Krishnapriyan et al., , Huang and Agarwal, , Wang et al., , , Rohrhofer et al., ].
Yet perhaps the greatest drawback of soft constraints is that they provide no guarantees that the physical constraints will be satisfied at inference time; at best, they can only hope to guide the model weights towards physically plausible solutions during training.
For certain safety-critical applications, where the ability to provably and robustly obey known physical constraints is crucial, this is inadequate.
These critical drawbacks motivate stricter enforcement of the constraints, both during training and at inference time, in the form of *hard* constraints.

##### Contribution.

In this paper, we introduce *projected neural differential equations* (PNDEs), a method for incorporating hard constraints in NDEs based on the projection of the vectors of the NDE’s vector field onto tangent spaces of the constraint manifold (see Figure 1).
Since the constraint manifold is defined by arbitrary algebraic constraints on the state of the system, our method can enforce diverse constraints, such as conservation laws, holonomic constraints, and even external forcings, all without requiring the system to be measured in specific (e.g., canonical) coordinates.
Our method enforces constraints more accurately than existing methods, requires fewer hyperparameters to tune, and enhances the numerical stability and generalizability of the systems we study.

##### Outline.

We motivate and derive the proposed method in Section 2, before surveying the existing literature and comparing to related works in Section 3.
Finally, in Section 4, we demonstrate empirically on several challenging examples, including highly chaotic systems and elaborate power grid models, that our approach satisfies the constraints more accurately than state-of-the-art baselines, all while requiring fewer hyperparameters and permitting larger timesteps of the numerical solver.

## 2 Projected Neural Differential Equations

### 2.1 Background

We first review some differential geometry definitions that are relevant to the proposed approach and refer the reader to [Lang, , Lee, , Boumal, ] for a more detailed treatment.
The set of tangent vectors to a manifold $\mathcal{M}$ at a point $u\in\mathcal{M}$ is a vector space called the *tangent space* to $\mathcal{M}$ at $u$ and is denoted by $\text{T}_{u}\mathcal{M}$. The disjoint union of all the tangent spaces to $\mathcal{M}$ forms the *tangent bundle* $\text{T}\mathcal{M}=\{(u,v)\ |\ u\in\mathcal{M},\ v\in\text{T}_{u}\mathcal{M}\}$ of $\mathcal{M}$. A *vector field* on a manifold $\mathcal{M}$ is a continuous map $X:\mathcal{M}\rightarrow\text{T}\mathcal{M}$ with the property that $X(u)\in\text{T}_{u}\mathcal{M}$ for all $u\in\mathcal{M}$. The *integral curve* at $u$ of a vector field $X$ on $\mathcal{M}$ is the smooth curve $\gamma$ on $\mathcal{M}$ such that $\gamma(0)=u$ and $\gamma^{\prime}(t)=X(\gamma(t))$ for all $t$ in the interval of interest. The integral curve $\gamma$ is simply the solution of a first-order ODE with right-hand side $X_{\gamma(t)}$ for $t$ in the interval of interest, so theorems for the existence, uniqueness, and smoothness of ODE solutions translate directly to integral curves of vector fields on manifolds.

### 2.2 Problem Setting

##### Neural Differential Equations.

Assume we have an ambient $n$-dimensional Euclidean phase space $\mathcal{E}$, and consider the neural differential equation

$$ $\dot{u}=f_{\theta}(u(t),t)\qquad u(0)=u_{0},$ (1) $$

where $u:{\mathbb{R}}\to\mathcal{E}$, $u_{0}\in\mathcal{E}$ is an initial condition, $\dot{u}=du/dt$, and $f_{\theta}$ is parameterized by a neural network with parameters $\theta\in{\mathbb{R}}^{d}$.

Given observations $\{t_{i},u(t_{i})\}_{i=1}^{N}$, the parameters $\theta$ can be optimized using stochastic gradient descent by integrating a trajectory $\hat{u}(t)=\mathrm{ODESolve}(u_{0},f_{\theta},t)$ and taking gradients of the loss $\mathcal{L}(u,\hat{u})$ with respect to $\theta$.
Gradients of the $\mathrm{ODESolve}$ operation can be obtained using adjoint sensitivity analysis (sometimes called *optimize-then-discretize*) or via direct automatic differentiation of the solver operations (*discretize-then-optimize*) [Chen et al., , Kidger, ].

Since any non-autonomous system can be equivalently represented as an autonomous system by adding time as an additional coordinate, for ease of notation we will drop the time-dependence in Equation 1 and consider only autonomous systems from now on, without loss of generality.

##### Constraints.

Suppose that the state $u(t)$ of the system must be constrained to the set

$$ $\mathcal{M}=\{u\in\mathcal{E}\,;\,g(u)=0\},\qquad\text{where}\quad g\colon \mathcal{E}\to{\mathbb{R}}^{m},\ \ g:u\mapsto\big{(}g_{1}(u),\ldots,g_{m}(u) \big{)},$ (2) $$

given $m<n$ independent explicit algebraic constraints $g_{i}(u)=0$. Assuming that the Jacobian $Dg(u)$ exists and has full rank for all $u\in\mathcal{M}$, the Regular Level Set Theorem [Lee, , #5.14] implies that $\mathcal{M}$ is a smooth $(n-m)$-dimensional embedded submanifold of ${\mathbb{R}}^{n}$. To satisfy the constraints $g(u)=0$, we seek solutions of Equation 1 that are constrained to the embedded submanifold $\mathcal{M}$, which we refer to as the *constraint manifold*. In other words, we seek integral curves on $\mathcal{M}$.

### 2.3 Proposed Approach

##### Proposed Approach.

We want to learn $f_{\theta}:\mathcal{M}\to\mathrm{T}\mathcal{M}$ as a smooth vector field on the embedded submanifold $\mathcal{M}$ with local defining function $g$. We realize this by first parameterizing a smooth extension $\bar{f}_{\theta}\colon\mathcal{E}\to\mathcal{E}$ of $f_{\theta}$ in Euclidean space using a neural network, and then projecting $\bar{f}_{\theta}(u)$ using a suitable projection operator $\mathrm{Proj}_{u}\colon\mathcal{E}\to\mathrm{T_{u}}\mathcal{M}$ to obtain the desired tangent vector,

$$ $f_{\theta}(u)=\mathrm{Proj}_{u}\bigl{(}\bar{f}_{\theta}(u)\bigr{)}\in\mathrm{T _{u}}\mathcal{M}.$ (3) $$

Since $\mathcal{M}$ is a smooth embedded submanifold of $\mathcal{E}$, the restriction $f_{\theta}=\bar{f}_{\theta}|_{\mathcal{M}}$ is also smooth.
Note that the projection acts fibrewise (i.e. acting "point-by-point"); for each $u\in\mathcal{M}$, it projects the vector $\bar{f}_{\theta}(u)$ at $u$ onto $\mathrm{T_{u}}\mathcal{M}$, as opposed to projecting directly the full vector field $\bar{f}_{\theta}$ onto $\mathrm{T}\mathcal{M}$.

##### Construction of the projection operator.

We now derive a suitable projection operator, following the approach of Boumal [].
Given an inner product $\langle\cdot,\cdot\rangle$ and induced norm $\lVert\cdot\rVert$ on $\mathcal{E}$, the differential $\mathrm{D}g(u)$ and its adjoint $\mathrm{D}g(u)^{*}$ define linear maps

$$ $\mathrm{D}g(u)[v]=\bigl{(}\langle\nabla g_{1}(u),v\rangle,\ldots,\langle\nabla g _{m}(u),v\rangle\bigr{)},\qquad\mathrm{D}g(u)^{*}[\alpha]=\sum_{i=1}^{m}\alpha _{i}\nabla g_{i}(u).$ (4) $$

Denoting the orthogonal projection from $\mathcal{E}$ to $\mathrm{T}_{u}\mathcal{M}$ by $\mathrm{Proj}_{u}:\mathcal{E}\to\mathrm{T}_{u}\mathcal{M}$, we can uniquely decompose any vector $v\in\mathcal{E}$ into its components parallel and perpendicular to $\mathrm{T}_{u}\mathcal{M}$,

$$ $v=\mathrm{Proj}_{u}(v)+\mathrm{D}g(u)^{*}[\alpha].$ (5) $$

The coefficients $\alpha$ are obtained as the solution of the least-squares problem

$$ $\alpha=\operatorname*{arg\,min}_{\alpha\in{\mathbb{R}}^{m}}\lVert v-\mathrm{D} g(u)^{*}[\alpha]\rVert^{2}=\bigl{(}\mathrm{D}g(u)^{*}\bigr{)}^{\dagger}[v],$ (6) $$

where $(\cdot)^{\dagger}$ denotes the Moore-Penrose pseudoinverse.
Substituting Equation 6 into Equation 5, we obtain a formula for the orthogonal projection from $\mathcal{E}$ to $\mathrm{T}_{u}\mathcal{M}$,

$$ $\boxed{\mathrm{Proj}_{u}(v)=v-\mathrm{D}g(u)^{*}\left[\bigl{(}\mathrm{D}g(u)^{ *}\bigr{)}^{\dagger}[v]\right].}$ (7) $$

##### Projected Neural Differential Equation.

Projecting the NDE in Equation 1 using the projection operator defined in Equation 7 gives the following *Projected Neural Differential Equation* (PNDE),

$$ $\boxed{\dot{u}=\mathrm{Proj}_{u}(f_{\theta}(u))\in\text{T}_{u}\mathcal{M},}$ (8) $$

which is a differential equation on the constraint manifold $\mathcal{M}$.

###### Proposition 1 .

Solutions to the PNDE 8 with $u_{0}\in\mathcal{M}$ satisfy $g(u(t))=0$ for all $t\geq 0$.

###### Proof.

We note that

$$ $\frac{dg(u(t))}{dt}=Dg(u)[\dot{u}]=Dg(u)\left[\mathrm{Proj}_{u}(f_{\theta}(u)) \right]=0$ (9) $$

since $\mathrm{Proj}_{u}(f_{\theta}(u))\in\text{T}_{u}\mathcal{M}=\ker\mathrm{D}g(u)$. As a result, $g(u(t))$ remains constant in time and satisfies $g(u(t))=g(u(0))=0$ for all $t\geq 0$ since $u_{0}\in\mathcal{M}$.
∎

###### Proposition 1 .

###### Proof.

## 3 Related Work

##### Hard Constraints.

A number of projection-based approaches have been proposed to constrain the outputs of deep learning models.
Harder et al. [] design constraint layers to ensure mass and energy conservation during downscaling (i.e., super-resolution) of climate model fields from low to high resolution.
Negiar et al. [] and Chalapathi et al. [] add a differentiable implicit constraint layer for finding optimal linear combinations of learned basis functions that approximately satisfy differential constraints at collocation points.
Jiang et al. [] and Duruisseaux et al. [] use a spectral projection layer to enforce linear differential constraints efficiently in Fourier space.
Hard constraints have also been use to take advantage of structured and well-understood dynamical systems, for example, the divergence-free property of incompresible fluid flows [Richter-Powell et al., , Mohan et al., , Xing et al., ] and the symplectic structure of Hamiltonian systems  [Lutter et al., , Greydanus et al., , Zhong et al., , Jin et al., , Burby et al., , Cranmer et al., , Sæmundsson et al., , Valperga et al., , Duruisseaux et al., , ].

The work of Finzi et al. [] is especially relevant to this paper.
They use explicit algebraic constraints to enable the learning of Hamiltonian dynamics in Cartesian coordinates, where the Hamiltonian is “simpler” and therefore easier to learn than in the original generalized coordinates.
They add constraints to the Hamiltonian via Lagrange multipliers, and use a variational approach to derive the corresponding constrained equations of motion.
The result is the standard Hamiltonian equations of motion, multiplied by a projection operator that enforces the constraints in a manner that is consistent with the underlying Hamiltonian structure.
While derived from an entirely different perspective, our approach can be understood as a principled generalization of their method beyond the idealized Hamiltonian setting to a much more general class of NDEs.
In this paper, for example, we include tests on a damped pendulum system, which is not compatible with their method.

##### Stabilized Neural Differential Equations.

White et al. [], whose work is most closely related to ours, constrain NDE solutions using stabilized neural differential equations (SNDEs), given by

$$ $\dot{u}=f_{\theta}(u)-\gamma F(u)g(u),$ (10) $$

where $\gamma\geq 0$ is a scalar, $F:{\mathbb{R}}^{n}\to{\mathbb{R}}^{n\times m}$ is a stabilization matrix, and $g:{\mathbb{R}}^{n}\to{\mathbb{R}}^{m}$ is the same constraint function we defined in Equation 2.
Under mild conditions on $\gamma$ and $F$, the stabilized NDE (10) admits all solutions of the original NDE (1) while rendering the constraint manifold provably asymptotically stable.
In other words, the stabilized NDE (10) is able to learn the correct ground truth dynamics while guaranteeing that the constraints are satisfied.

To understand how SNDEs relate to our proposed method, we first provide the intuition behind Equation 10.
The stabilization term (i.e. the second term) equals zero on the manifold, since $g(u)=0$ there by definition.
Hence, the stabilization is only “activated” when a trajectory leaves the constraint manifold, i.e., when the constraints have already been violated.
The effect of the stabilization is to “nudge” such a trajectory “back towards $\mathcal{M}$”, i.e., back towards satisfying the constraints, with the strength of the nudge determined by the stabilization parameter $\gamma$.
Larger values of $\gamma$ lead to a stronger nudge and stricter enforcement of the constraints.
However, excessively large values of $\gamma$ may require the numerical solver to stop more frequently, with the system eventually becoming stiff.
The key difference with our method is that SNDEs continuously correct violations of the constraints, while our method prevents such violations from occurring in the first place.
Accordingly, we can expect our method to yield a stronger enforcement of the constraints.
We will use SNDEs as a baseline for comparison throughout this paper.

##### Continuous Normalizing Flows on Riemannian Manifolds.

A large body of related work has studied neural ODEs on Riemannian manifolds, primarily in the context of continuous normalizing flows on non-Euclidean geometries.
In contrast to our approach, they typically deal with prototypical Riemannian manifolds, such as spheres and tori, and are not easily adapted to more general manifolds.
For example, the approach of Lou et al. [] requires a local chart for the manifold, the derivation of which is likely to be analytically intractable for embedded submanifolds with arbitrary and nonlinear local defining functions.
Similarly, the approaches of Bose et al. [] and Rezende et al. [] require an expression for the exponential map, while Mathieu and Nickel [] uses a custom geodesic distance layer, none of which are easily applied in our setting.

##### Optimization and Numerical Integrators on Manifolds.

The use of projections is widespread in optimization on manifolds [Absil et al., , Boumal, ], where it is used to project iterates to the constraint manifold and/or gradients onto appropriate tangent spaces (as in our approach).
The idea of projecting gradients onto tangent spaces of the constraint set dates back at least to Luenberger [, ], and is still commonly used. For instance, in Riemannian optimization, Euclidean gradients are typically replaced by Riemannian gradients, which are their orthogonal projections onto tangent spaces.
Explicit formulas and efficient projection algorithms have been derived and used for manifolds of greater practical interest [Edelman et al., , Manton, , Absil et al., , Wen and Yin, , Absil and Malick, , Hauswirth et al., , Zhang and Sra, , Oviedo León and Dalmau-Cedeño, , Oviedo León et al., ].
Where these methods differ crucially from our approach is that they typically seek to constrain the weights or optimization variables themselves to a given matrix manifold, while
our aim is to constrain trajectories of an NDE, the weights of which are not themselves directly constrained.
Regarding numerical integration on manifolds, the available literature is much less extensive [Hairer, ], and most approaches restrict themselves to Lie groups (very structured Riemannian manifolds) and Hamiltonian or mechanical systems, in order to obtain practical algorithms [Iserles et al., , Christiansen et al., , Celledoni et al., , ].
Some of these approaches are based on projection, similarly to optimization on manifolds, while others use local parametrizations or Lagrange multipliers.
Riemannian optimization algorithms have also been obtained by simulating dynamical systems using numerical integrators constrained to the relevant Riemannian manifolds [Alimisis et al., , Tao and Ohsawa, , Lee et al., , França et al., , Duruisseaux and Leok, , , , ].

## 4 Experiments

We now test the proposed approach on several challenging examples, including highly chaotic systems and elaborate power grid models. The details of the data generation, model setup, and training procedure, are provided in Appendices A and B. We evaluate each model’s predicted trajectory $\hat{u}(t)$ versus the ground truth trajectory $u(t)$ using the usual relative error $E(t)=\|u(t)-\hat{u}(t)\|_{2}/\|u(t)\|_{2}$.
For chaotic dynamics, we instead use the bounded relative error $E(t)=\|u(t)-\hat{u}(t)\|_{2}/\left(\|u(t)\|_{2}+\|\hat{u}(t)\|_{2}\right)$, which approaches 1 as trajectories become orthogonal.
All test statistics are averaged over 100 unseen initial conditions.

### 4.1 Conservation Laws

The Fermi–Pasta–Ulam–Tsingou (FPUT) lattice system [Fermi et al., ] describes a vibrating string with nonlinear dynamics.
The displacement $x_{j}$ of a segment $j$ with mass $m$ is governed by the second order equation,

$$ $m\ddot{x}_{j}=k(x_{j+1}+x_{j-1}-2x_{j})\left[1+\alpha(x_{j+1}-x_{j-1})\right],$ (11) $$

where the first factor is simply Hooke’s law with spring constant $k$, and the second factor introduces a nonlinearity governed by the parameter $\alpha$.

In the case of a circular lattice with $N$ segments and periodic boundary conditions $x_{i}=x_{i+N}$, the system conserves energy

$$ $E(x,\dot{x})=\sum_{j=1}^{N}\left[\frac{1}{2}m\dot{x}_{j}^{2}+\frac{1}{2}k(x_{j +1}-x_{j})^{2}+\frac{1}{3}\alpha(x_{j+1}-x_{j})^{3}\right],$ (12) $$

where $x=(x_{1},\ldots,x_{N})$ and $\dot{x}=(\dot{x}_{1},\ldots,\dot{x}_{N})$.
We convert Equation 11 into a first order ODE with state $u=(x,v)\in{\mathbb{R}}^{2N}$ where we introduce the variable $v\coloneqq\dot{x}$.
We suppose that we possess prior physical knowledge of the conservation law in Equation 12, which we enforce as the constraint

$$ $g(u)=E(u)-E(u_{0})=0,$ (13) $$

where $g:{\mathbb{R}}^{2N}\to{\mathbb{R}}$ and $u_{0}$ is the initial condition.

Figure 2 compares PNDEs with SNDEs and an unconstrained NDE on the FPUT system with $N=128$ segments.
The unconstrained NDE is unstable and diverges for long rollouts (longer than the training trajectories and significantly longer than an individual training sample).
In contrast, all of the constrained models are stable, highlighting the value of enforcing known constraints in NDE models.
PNDEs and SNDEs demonstrate similar ability to predict the evolution of the system state, although they both eventually depart from the ground truth trajectory of this nonlinear system.
The PNDE model enforces the constraints several orders of magnitude more accurately than both SNDE models, despite requiring significantly fewer evaluations of the the right-hand side function $f_{\theta}$.

Figure: Figure 2: Performance of PNDEs compared to SNDEs and an unconstrained NDE for the Fermi–Pasta–Ulam–Tsingou (FPUT) lattice system [Fermi et al., ] with $N=128$ segments. Panel (a) shows the ground truth as well as the NDE, SNDE, and PNDE predictions of the system state. Panel (b) shows the relative error. Panel (c) shows the error in the constrained quantity. Note that PNDE and SNDE results in (a) and (b) are visually indistinguishable.
Refer to caption: x1.png

### 4.2 Simplifying Complex Learning Tasks Using Explicit Constraints

Figure: Figure 3: Illustration of a $3$-pendulum.

We now consider an $N$-pendulum system, consisting of $N$ simple pendulums connected end-to-end, as illustrated on the right for $N=3$.
$N=1$ corresponds to the canonical simple pendulum, $N=2$ yields the well-known chaotic double pendulum system, and higher values of $N$ yield systems with progressively more complex and chaotic dynamics.
For realism, and in contrast to much of the related literature, we include damping proportional to the relative angular velocity of each pendulum bob; this could represent, for example, friction in the hinges of the pendulum.

Assuming that each bob has a mass of \qty1 and each arm has a length of \qty1, the equations of motion of the $i$-th bob are then given by

$$ $\sum_{j=1}^{N}a(i,j)\bigl{[}\ddot{\theta}_{j}\cos(\theta_{i}-\theta_{j})+\dot{ \theta}_{j}^{2}\sin(\theta_{i}-\theta_{j})\bigr{]}+b(\dot{\theta}_{i}-\dot{ \theta}_{i-1})+(N-i+1)\mathfrak{g}\sin\theta_{i}=0,$ (14) $$

where $\theta_{i}$ is the angle that the $i$-th arm makes with the vertical, $b$ is the coefficient of friction, $a(i,j)=N-\max(i,j)+1$, and $\mathfrak{g}=\qty{9.81}{m.s^{-1}}$ is acceleration due to gravity.

Given observations of this system, a reasonable approach would be to learn an NDE directly approximating Equation 14, which is in generalized coordinates.
Alternatively, and analogously to the approach taken by Finzi et al. [] for Lagrangian and Hamiltonian systems, we can first convert to Cartesian coordinates, where the equations of motion may be easier to learn.
In Cartesian coordinates, constraints arise naturally from the rigidity of the pendulum arms, i.e., the length of each arm should remain constant in time.
Hence, given the position $q_{i}=(x_{i},y_{i})$ of the $i$-th bob in Cartesian coordinates, we can write down $N$ holonomic constraints on the system,

$$ $g^{(pos)}_{i}(q)=||q_{i}-q_{i-1}||^{2}-1=0,\quad i=1,\ldots,N,$ (15) $$

where $q=(q_{1},\ldots,q_{N})$ and $q_{0}=(0,0)$ is the origin.
We can differentiate these constraints to obtain $N$ additional conjugate constraints on the velocity,

$$ $g^{(vel)}_{i}(q)=(q_{i}-q_{i-1})\cdot(\dot{q}_{i}-\dot{q}_{i-1})=0,\quad i=1, \ldots,N.$ (16) $$

Collecting the position and velocity constraints into a vector, we obtain $2N$ constraints for the $N$-pendulum in Cartesian coordinates:

$$ $g(q)=\left[g^{(pos)}_{1}(q),\ldots,g^{(pos)}_{N}(q),\ g^{(vel)}_{1}(q),\ldots, g^{(vel)}_{N}(q)\right]=0,$ (17) $$

where $g:{\mathbb{R}}^{4N}\to{\mathbb{R}}^{2N}$.
We note that the generalized coordinates of Equation 14 satisfy these constraints automatically.

Figure 4 demonstrates that learning constrained dynamics in Cartesian coordinates produces significantly better predictions than learning the same dynamics expressed in generalized coordinates, where the constraints are satisfied automatically.
The benefits of learning in Cartesian coordinates with constraints becomes especially pronounced as the size and complexity of the system increases.
Amongst constrained models, PNDEs enforce the constraints several orders of magnitude more accurately than SNDEs and do not require tuning of the stabilization parameter $\gamma$.
Unconstrained NDEs in Cartesian coordinates are numerically unstable, again highlighting the importance of constraints.

Figure: Figure 4: Damped $N$-pendulum system. We compare the ability of an NDE to learn the dynamics of the damped $N$-pendulum system in generalized coordinates, where the constraints are satisfied automatically, versus Cartesian coordinates, where we explicitly enforce the constraints. Top: Across all systems, constrained models in Cartesian coordinates better predict the future evolution of the system’s state. The difference in performance becomes more pronounced as $N$ is increased and the corresponding equations of motion become more complex in generalized coordinates. Bottom: PNDEs enforce the constraints exactly, while SNDEs incur large errors, especially for larger values of $N$. For SNDEs, smaller values of $\gamma$ than those shown here do not enforce the constraints, while larger values are computationally prohibitive due to stiffness. In generalized (angular) coordinates, the constraints are satisfied automatically by the choice of coordinate system. Nevertheless, the results show that converting to Cartesian coordinates and adding constraints leads to significantly better generalization. Unconstrained NDEs in Cartesian coordinates (not shown) are numerically unstable, highlighting the necessity of imposing explicit constraints via SNDEs or PNDEs.
Refer to caption: x2.png

### 4.3 Power Grids

The IEEE test systems are standardized power system models widely used in electrical engineering research and education [Christie, ].
Representing key features of power grids such as topologies and operating conditions, these models enable researchers to benchmark algorithms and offer a common framework for analyzing power flow, stability, and reliability in electrical networks.
We will study the dynamics of two test systems in particular: the IEEE 14-Bus system (see Figure 5, right) and the IEEE Reliability Test System 1996 (hereafter RTS-96) [Grigg et al., ].
We adapt the IEEE test cases to represent future grid dynamics by replacing conventional generators with renewable energy sources.
Our renewables-based grids consist of three types of buses, or nodes: (1) renewable energy sources, (2) loads (that is, consumers), and (3) slack buses (a constant voltage bus representing, for example, a large power plant or a connection to another grid).

Figure: Figure 5: The IEEE 14-Bus system. Circles represent generators and arrows represent consumers. Diagram extracted from Leon et al. [].
Refer to caption: extracted/5967760/figures/ieee14bus.png

The state of each bus in the grid is determined by a complex voltage, the real and imaginary part of which we learn separately.
That is, for a grid consisting of $N$ buses, we learn the coupled dynamics of $2N$ voltages.
Normally, a power grid operates at a steady-state equilibrium, called its *operation point*.
We apply random perturbations to each grid (see Appendix A) and learn the dynamics from the resulting transients as it returns to its operation point.

In our grid models, loads are modeled using *PQ-buses*, meaning that the complex power at the PQ-bus $j$ is set to a fixed, constant value given by

$$ $P_{j}+\mathrm{i}Q_{j}=v_{j}\cdot i_{j}^{*}$ (18) $$

where $P_{j}$ and $Q_{j}$ are the fixed active and reactive power, respectively, $v_{j}$ is the complex voltage, and $i_{j}$ is the complex current [Machowski et al., ].
These fixed power values will be the constraints in our experiments.
The complex currents in Equation 18 are determined from the voltages by $i=\mathrm{LY}\times v$, where $i\in\mathbb{C}^{N}$, $v\in\mathbb{C}^{N}$, and the admittance Laplacian $\mathrm{LY}\in\mathbb{C}^{N\times N}$ describes the line admittances, which are well understood and assumed here to be known.
Given $N$ total buses and $M$ PQ-buses, the $2M$ constraints (one real and one imaginary for each PQ-bus) are given by

$$ $g(u)=\left[P_{1}\!-\!\operatorname{Re}{\!(v_{1}\cdot i_{1}^{*})},Q_{1}\!-\! \operatorname{Im}{\!(v_{1}\cdot i_{1}^{*})},\ldots,P_{M}\!-\!\operatorname{Re} {\!(v_{M}\cdot i_{M}^{*})},Q_{M}\!-\!\operatorname{Im}{\!(v_{M}\cdot i_{M}^{*} )}\right]=0,$ (19) $$

where $u=\left(\operatorname{Re}{(v_{1})},\operatorname{Im}{(v_{1})},\ldots,
\operatorname{Re}{(v_{N})},\operatorname{Im}{(v_{N})}\right)\in{\mathbb{R}}^{2N}$ and $g:{\mathbb{R}}^{2N}\to{\mathbb{R}}^{2M}$.
For the IEEE 14-bus system, $N=14$ and $M=9$, while for the IEEE RTS-96 system, $N=74$ and $M=40$.

Table 1 shows the mean-squared error for both IEEE test systems, calculated over the entire test set, and Figure 6 shows the relative error as a function of time, along with some illustrative trajectories.

Unconstrained NDEs incur large errors in the PQ-bus power constraint.
For SNDEs, increasing the stabilization hyperparameter $\gamma$ enforces the constraints more accurately, but larger values lead to stiffness and, ultimately, prohibitively expensive simulations.
Meanwhile, PNDEs enforce the constraints close to machine precision and several orders of magnitude more accurately than the best SNDEs without stiffness.

**Table 1: MSE of the predicted voltage $V$ (all buses) and power constraint (PQ-buses only), along with the number of function evaluations required by the adaptive ODE solver. For both systems, the constrained models predict the voltages at least as well as the NDE. Unsurprisingly, the NDE incurs large errors in the PQ-bus power. SNDEs enforce the constraints more accurately for larger values of $\gamma$, but the largest values lead to stiffness, requiring many more function evaluations. PNDEs efficiently enforce the constraints near machine precision.**
|  | IEEE 14-Bus | IEEE RTS-96 |  |  |  |  |
| --- | --- | --- | --- | --- | --- | --- |
| Model | $V$ ($\times 10^{-5}$) | PQ Power | $f_{\theta}$ Evals | $V$ ($\times 10^{-4}$) | PQ Power | $f_{\theta}$ Evals |
| NDE | $1.8\!\pm\!1.8$ | $(3.0\!\pm\!2.0)\!\times\!10^{-3}$ | 171 | $2.1\!\pm\!0.8$ | $(0.9\!\pm\!1.1)\!\times\!10^{-1}$ | 405 |
| SNDE ($\gamma\!=\!0.1$) | $1.7\!\pm\!1.8$ | $(1.5\!\pm\!1.3)\!\times\!10^{-3}$ | 171 | $2.0\!\pm\!0.8$ | $(5.0\!\pm\!6.0)\!\times\!10^{-2}$ | 405 |
| SNDE ($\gamma\!=\!1$) | $1.6\!\pm\!1.8$ | $(2.0\!\pm\!3.0)\!\times\!10^{-4}$ | 177 | $2.0\!\pm\!0.8$ | $(0.5\!\pm\!1.1)\!\times\!10^{-2}$ | 405 |
| SNDE ($\gamma\!=\!10$) | $1.6\!\pm\!1.8$ | $(0.5\!\pm\!1.3)\!\times\!10^{-5}$ | 297 | $1.9\!\pm\!0.8$ | $(1.7\!\pm\!7.0)\!\times\!10^{-4}$ | 411 |
| SNDE ($\gamma\!=\!100$) | $1.6\!\pm\!1.8$ | $(0.5\!\pm\!1.6)\!\times\!10^{-7}$ | 1929 | $1.9\!\pm\!0.8$ | $(2.0\!\pm\!8.0)\!\times\!10^{-6}$ | 1815 |
| PNDE | $1.6\!\pm\!1.8$ | $\mathbf{(3.5\!\pm\!1.6)}\!\times\!\mathbf{10^{-11}}$ | 171 | $1.9\!\pm\!0.8$ | $\mathbf{(1.9\!\pm\!1.5)}\!\times\!\mathbf{10^{-10}}$ | 405 |

Figure: Figure 6: Test statistics for the IEEE 14-Bus System (top) and IEEE RTS-96 System (bottom). Panels (a) and (d) show a representative voltage for a single trajectory. All models accurately capture the shape of this voltage for the 14-bus system, but appear to struggle with the high frequency components of the RTS-96 system. Panels (b) and (e) show the relative error of the voltages at all buses, calculated over the entire test set. Constrained models outperform unconstrained models, with the gap growing over time. Panels (c) and (f) show the relative error in the power at the PQ-buses, with the PNDE outperforming the best SNDEs by several orders of magnitude.
Refer to caption: x3.png

## 5 Conclusion

We have introduced projected neural differential equations (PNDEs), a new method for enforcing arbitrary algebraic constraints in neural differential equation models.
Compared to state-of-the-art baselines, our approach requires fewer hyperparameters and enforces the constraints several orders of magnitude more accurately.
In experiments on complex nonlinear and chaotic systems, as well as two challenging examples based on state-of-the-art power grid models, PNDEs are efficient, numerically stable, and generalize well to unseen data.
We have also demonstrated the use of explicit constraints to simplify and significantly improve the learning of complex dynamical systems in alternative coordinate systems.
Many promising directions remain for future research.
For example, graph neural networks are also compatible with our method and a natural fit for modeling grid dynamics, thereby enabling critical out-of-distribution tests such as changing grid topologies.

#### Acknowledgments

This work was funded by the Volkswagen Foundation.
The authors gratefully acknowledge the European Regional Development Fund (ERDF), the German Federal Ministry of Education and Research, and the Land Brandenburg for supporting this project by providing resources on the high performance computer system at the Potsdam Institute for Climate Impact Research.

## Appendix A Datasets

For each system, we generate ground truth trajectories from random initial conditions.
The exact number of initial conditions used for training, validation, and testing, as well as the duration of each trajectory and the timestep used, is determined based on the individual characteristics of each system, and summarized in Table 2.
For training and validation, each trajectory is split into contiguous segments consisting of 4 timesteps, which are then randomized and combined along the batch dimension.
The first observation in each segment is the initial condition for the prediction.

With the exception of the power grid experiments, all trajectories are generated using the 9(8) explicit Runge-Kutta algorithm due to Verner [], implemented in DifferentialEquations.jl [Rackauckas and Nie, ] as Vern9, with absolute and relative solver tolerances of $10^{-12}$.

The power grid models, which are differential-algebraic equations (DAEs), are implemented in the Julia package PowerDynamics.jl [Plietzsch et al., ].
Given a random initial condition, we then generate a ground truth trajectory by solving the DAE for 10 seconds using a 4th order L-stable Rosenbrock method [Hairer and Wanner, ], implemented in the Julia package DifferentialEquations.jl [Rackauckas and Nie, ] as Rodas4.

**Table 2: Summary of datasets and hyperparameters**
| System | $n_{\mathrm{train}}$ | $n_{\mathrm{valid}}$ | $n_{\mathrm{test}}$ | $T$ | $dt$ | Hidden Layers | Hidden Width | Epochs |
| --- | --- | --- | --- | --- | --- | --- | --- | --- |
| FPUT | 120 | 60 | 100 | 25 | 0.1 | 3 | 512 | 1000 |
| $N$-pendulum | 160 | 80 | 100 | 10 | 0.02 | 3 | 1024 | 5000 |
| IEEE 14-Bus | 100 | 50 | 100 | 10 | 0.1 | 3 | 512 | 5000 |
| IEEE RTS-96 | 100 | 50 | 100 | 10 | 0.1 | 3 | 512 | 5000 |

### A.1 Initial Conditions

##### Fermi–Pasta–Ulam–Tsingou System.

Given a string defined on the range $x\in[0,1)$, we generate a symmetric initial displacement according to $x_{0}\sim\exp{\left(-(x-\mu)^{2}/\sigma^{2}\right)}$, where $\mu=0.5$ and $\sigma\sim\mathrm{Unif}(0.1,0.3)$.
We release the string from rest, i.e., we set $\dot{x}_{0}=0$.

##### N 𝑁 N italic_N -Pendulum.

The initial displacements of the pendulum arms are drawn independently from $\mathrm{Unif}(0,2\pi)$ and the initial velocities are drawn from $\mathrm{Unif}(-1,1)$.

##### Power Grids.

Given a power grid at its operation point (i.e., a stable equilibrium with all variables at their desired steady-state values) we obtain perturbed initial conditions using the ambient forcing method of Büttner et al. [], which ensures that the perturbed state satisfies the PQ-bus constraints.
The transient behavior of the grid following the perturbation produces rich and complex dynamics that we aim to emulate using NDEs.

## Appendix B Training

All models are trained on NVIDIA H100 GPUs.
We minimize the mean squared error loss using the AdamW optimizer [Loshchilov and Hutter, ] with a constant learning rate of $10^{-3}$, weight decay of $10^{-4}$, and a batch size of 1024.
All systems are modeled using standard multilayer perceptrons with GELU activation functions to encourage smoothness.
The number of hidden layers, hidden width, and training epochs are summarized in Table 2.
All models are initially trained without constraints, and subsequently fine-tuned with the constraints for 100 epochs.

During training, all models are solved using the 5(4) explicit Runge-Kutta algorithm of Tsitouras [], implemented in DifferentialEquations.jl [Rackauckas and Nie, ] as Tsit5, with absolute and relative tolerances of $10^{-4}$.
Gradients of ODE solutions are obtained using the adjoint sensitivity method [Pontryagin, , Chen et al., ], implemented in SciMLSensitivity.jl as BacksolveAdjoint [Rackauckas et al., ].
For testing, we use absolute and relative tolerances of $10^{-6}$.

## Peer Reviews

confidence: 4
contribution: 3
initial_rating: 8
presentation: 4
questions: Can the authors add some "cartoons" to ensure that readers less inclined with differential geometry can still understand the main (pictorially easy to be fair) intuition behind the paper? 

Can the authors expand the discussion on relevant literature (maybe with some table comparison?) as per the "Weaknesses" section? 

Can you provide insights into the computational complexity of the projection operation and discuss potential scalability issues (needs not be FOCS-style theory). Include analysis on the numerical stability and error propagation introduced by the projection maybe if possible?

Re the comment above "The experiments focus on systems where the constraint manifold is relatively straightforward to compute. It would be valuable to test PNDEs on "less trivial" systems with high-dimensional constraints or where the constraint manifold has nontrivial topology maybe? " can you maybe design such a larger instance and maybe harder topology problem? I would not decline the paper based on this but I do think it would massively strengthen the paper. 

Note the unique typos I found was that no $ sign was used in a couple of T_uM instances, just make sure you fix this. 

By addressing these points, the paper would be strengthened in terms of positioning within the existing literature, theoretical rigor, empirical validation, and clarity, ultimately enhancing its significance and impact on the field and would make it a super strong paper for ICLR 2025.
soundness: 3
strengths: The paper presents a substantial advancement in the field of neural differential equationss given that the authors address a, indeed, crucial limitation of standard NDEs—the inability to enforce known constraints in the learned dynamical systems—which often leads to poor generalization and numerical instability. While many solutions to this have been discussed, e.g. within Hamiltonian neural nets, the strengths of the paper are several and go beyond the literature, to the best of my knowledge. c

The introduction of PNDEs is novel contribution and the authors provide a principled method to incorporate hard constraints directly into the learning process. Their approach differs from e.g., Stabilized Neural Differential Equations by ensuring that the constraints are satisfied exactly, rather than asymptotically or approximately. Using projection operators is creative and less common within "deep learning" as opposed to more traditional convex optimization.

The paper demonstrates rigorous, but bried, theoretical development and provides clear mathematical derivations. The authors derive the projection operator using the Jacobian of the defining constraint functions, and in detail explain and ensure that the projected vector field remains within the tangent space of M and they further go to propose in detail that solutions to the PNDE remain on the constraint manifold is well-founded, and the proof is succinct yet thorough. Some extra discussion as well as graphical illustration would be beneficial here! The experimental section is robust, and covers a range of systems.

In terms of clarity, I am happy to see a quite clear and well-written paper which manages to convey complex ideas effectively. That said, certain sections would be harder to follow for people with more ML/DL background but I won't count this as a limitation. Note that the motivation behind enforcing hard constraints in NDEs is very clearly articulated, and the limitations of existing methods are adequately discussed. The derivation of the projection operator is presented step-by-step (again, a graphical illustration would do miracles here), making it accessible to mathematically include readers with a background in differential geometry and dynamical systems. The experimental figures and tables are informative and enhance the understanding of the results. The experimental setup is described in sufficient detail, allowing for reproducibility (although I could not locate a link with a repo).


As mentioned earlier, incorporating hard constraints into NDEs has significant implications for modeling realistic dynamical systems that inherently possess constraints, such as conservation laws and algebraic relationships. The ability of PNDE to enforce these constraints exactly enhances the reliability and accuracy of the models, which is crucial in safety-critical applications like power grid management. However, many practicioners, especially for this example, would claim that the lack of rigorous guarantees is a problem. Returnign to the main ideas of the paper, when suitable examples are considered, improving generalization and numerical stability, PNDEs contribute to advancing the state-of-the-art in data-driven modeling of dynamical systems. This work arguably opens up new possibilities for applying NDEs to a broader class of problems where constraints play a vital role.
summary: The paper addresses the challenge of learning constrained dynamical systems in the context of neural differential equations (i.e., NDEs herafter). This term of NDEs includes the 2018 class of NODEs and generalization therefore such UDEs and bagging them together the authors are interested since indeed they allow for flexible modeling of dynamical systems by parameterizing the vector field f_\theta with neural networks. The authors observe however, that they do not inherently enforce possible known constraints that the system may have, such as conservation laws, holonomic constraints (both applying quite well in Hamiltonian systems for example), or algebraic relationships and this can lead to learned models that to different levels of severity may violate essential properties of the system, resulting in poor generalization and numerical instability.

To overcome this limitation, the authors introduce the "Projected Neural Differential Equations" (PNDEs) whose key idea is to enforce constraints by projecting the neural network vector field f_\theta onto the tangent space T_M of the constraint manifold M, which is defined by (algebraic equations) g(u) = 0. Specifically, they define the projected vector field as Proj_u (f_\theta) \in T_uM
where \mathrm{Proj}_u is the orthogonal projection operator of from T_uE onto T_uM. By integrating this projected vector field, the solutions remain on the manifold M, ensuring that the constraints are satisfied for all time. So in some sense, they try to "get rid" of the components of the vector field/neural net that would be learned but would live outside the constrained submanifold the physical system actually lives on.

The authors provide a detailed derivation of the projection operator using common decomposition techniques. They demonstrate that for an embedded submanifold M defined by smooth constraint functions g(u) the projection can be explicitly computed using the Jacobian of the latter, i.e., the constraints. This allows for efficient computation of the projected vector field during numerical integration.

To validate their approach, the authors conduct experiments on several challenging dynamical systems: the Fermi–Pasta–Ulam–Tsingou lattice system, a damped pendulum system with holonomic constraints and power grid models that incorporate power flow equations. 

Compared to various existing methods they cite, such as SNDEs,  the experiments seem to verify claims that the proposed method offers exact constraint enforcement without introducing additional hyperparameters and/or suffering from stiffness issues in numerical integration which arguably distinguishes PNDEs from penalty-based methods or those that incorporate constraints as soft losses during training, which may not guarantee constraint satisfaction during inference.

Overall, the paper presents a principled and general framework for incorporating hard constraints into neural differential equations by projecting the neural network vector field onto the tangent space of the constraint manifold. This approach enhances the modeling of constrained dynamical systems, improving accuracy, generalizability, and numerical stability.
weaknesses: While the paper, as discussed, presents a novel and effective method for incorporating hard constraints into neural differential equations there are several areas where the work could be improved.

One of the (few) main weaknesses lies in the discussion of related work and positioning of the proposed method within the existing literature. The paper focuses primarily on comparing PNDEs to SNDEs. However, there is a rich body of research on incorporating physical constraints and conservation laws into neural network models of dynamical systems that is not adequately addressed.

For instance, the (indeed) cited HNNs [Greydanus et al., 2019] and Symplectic ODE-Net [Zhong et al., 2020] (since the author do mention inductive bias in the intro) are significant contributions that leverage the symplectic structure of Hamiltonian systems to enforce conservation of energy and other invariants. These methods learn the Hamiltonian function directly and ensure that the learned dynamics preserve the symplectic form, inherently satisfying certain physical constraints. Therefore, it's not clear to as wether the PNDEs would be relevant in systems where HNNs seem to perform very well. As a matter of fact, in recent work on learning generalized Hamiltonians using fully symplectic mappings [Choudhary et al. 2024], addresses the challenge of modeling non-separable Hamiltonian systems using implicit symplectic integrators within neural networks which should be a class of problems where previously I would had assumed PNDEs to be prime candidates to work on but its just not clear to me as to what the best approach would be in such situations. So, overall, I would prefer a more thoorough discussion here. Finally, I would be keen for the authors to portray further unerstanding of the literature of projections. For example, it is known that such projections, introduce certain symmetries. These symmetries ideally can be quotioned out in order to fascillitate easier training, see for example a similar construction in convex optimization and SDPs where the tangential projection symmetries need be addressed [Bellon et. al. 2210.08387].


While the paper provides a clear derivation of the projection operator and proves that solutions to the PNDE remain on the constraint manifold, the theoretical analysis could be strengthened. Specifically, the paper lacks a discussion on the computational complexity and scalability of the projection operation in high-dimensional systems or with complex constraints. Maybe it's too hard? From the practicioner's point of view this is important too. Given the experiments discuss power grid (we would use BnC methods normally and not gradient based methods for a number of reasons) this is important. Also, computing the projection onto the tangent space requires solving a system involving the Jacobian of the constraints, which can be computationally intensive for large-scale systems.

Moreover, the paper does not provide theoretical guarantees on the convergence or stability of the PNDEs beyond the preservation of the constraints. Are there some assumptions that can be made that would allow for an analysis of the numerical errors introduced by the projection and their impact on the overall solution accuracy? Additionally, insights into how the method performs under approximate constraints or in the presence of noise would enhance the understanding of its robustness.


Re the experimental section, while demonstrating the effectiveness of PNDEs on several systems, could be expanded to provide a more comprehensive evaluation. The experiments focus on systems where the constraint manifold is relatively straightforward to compute. It would be valuable to test PNDEs on "less trivial" systems with high-dimensional constraints or where the constraint manifold has nontrivial topology maybe? 

Furthermore, the comparison is primarily with SNDEs and unconstrained NDEs. Including additional baseline methods, such as HNNs, Symplectic Neural Networks, or other constraint-enforcing techniques, would strengthen the empirical evaluation. This would provide a clearer picture of the advantages and limitations of PNDEs relative to existing approaches.

While the paper is generally well-written, certain sections could be clarified for better accessibility. The derivation of the projection operator, although mathematically rigorous, might be challenging for readers not deeply familiar with differential geometry. Providing more intuitive explanations or illustrative examples could help bridge this gap.

Additionally, the notation used in some equations, such as the use of adjoints and pseudoinverses, could be explained in more detail. Ensuring that all symbols and operations are clearly defined would improve the readability of the paper.


Another point is that the proposed method assumes that the constraints can be expressed as explicit algebraic equations and that the Jacobian of the constraints is full rank. In practice, many systems might have constraints that are implicit, differential-algebraic, or have singular Jacobians. What happens then? Discussing how PNDEs could be extended or adapted to handle such cases would enhance the significance and applicability of the work.

---

confidence: 5
contribution: 1
initial_rating: 1
presentation: 2
questions: My concerns regarding this paper are as explained above.
soundness: 2
strengths: The paper is well-written and easy to read. Some experiments were conducted to support the effectiveness of the proposed method.
summary: In this paper, a method for learning differential equations while preserving conservation laws is proposed. Specifically, to preserve conservation laws, the authors project the learned vector field onto the tangent bundle of the manifold defined by the conserved quantities.
weaknesses: The proposed method is not novel because this method has already been proposed; the continuous-time model is shown in [1], and also the discrete-time model is shown in [2]. 

[1] Kasim, M.F. and Lim, Y.H. (2022) Constants of motion network, NeurIPS2022

[2] Matsubara, T. and Yaguchi, T. (2023) FINDE: Neural Differential Equations for Finding and Preserving Invariant Quantities, ICLR2023

In [1], the learned vector field is designed to be orthogonal to the gradient vectors of the conserved quantities. Precisely, the learned vector field is projected onto the tangent space at each point of the manifold defined by the conserved quantities, which is the same as the approach proposed in this paper. In [1], the QR decomposition is used for orthogonalization and hence the method of computing the projection operator is a little different from that of this paper, which uses the pseudo inverse.

In [2], exactly the same approach as this paper is proposed; in [2], the manifold defined by conserved quantities is first introduced. Then, they consider tangent bundles of this manifold, and project the leaned vector field onto the tangent space at each point. More precisely, in [2], a continuous-time model is first considered. Equation (6) in [2], which represents the continuous-time model, is completely identical to (7) in this paper. The pseudo-inverse matrix is used for projection in [2], the pseudo-inverse matrix is specifically computed though (so the model looks a little different.) In addition, in [2] a discrete-time model is also discussed. In the discrete-time model, the discrete gradient, which is a discrete version of the gradient, is considered, and the discrete tangent space is defined using the discrete gradient. The discrete-time model is essentially the projection onto this discrete tangent space.

In addition, it seems that the conserved quantities are assumed to be given in this paper; however, the methods shown in the above papers can handle cases where these quantities are unknown. 

Considering the above, the contributions of this paper is quite limited.

---

confidence: 3
contribution: 1
initial_rating: 5
presentation: 2
questions: Please address the concerns mentioned in the Weaknesses above.
soundness: 2
strengths: Constrained dynamics are present in real-world problems, and existing NODEs that do not consider these constraints run the risk of failing to satisfy them. In contrast, the proposed method achieves hard constraint satisfaction through projection.
summary: This paper presents a projection-based approach to ensure hard constraint satisfaction in constrained dynamics, which is important for several real-world applications. However, several concerns outlined in the weaknesses section remain, and further investigation of these limitations is needed to improve its practical applicability.
weaknesses: 1.	The methodology is based on the assumption that the analytic form of the constraint function $g$ is known, which seems impractical. In real situations where the dynamics are unknown and the states $u$ are given only by the data, there are often many cases where the analytic form of $g$ is also unknown.
2.	What is the difference and advantage of the proposed method of projecting the forcing function $f$ onto the tangent space compared to projecting the predicted state $u$ onto the constraint manifold? Projecting $u$ onto the constraint manifold seems simpler than projecting $f$ onto the tangent space, while still satisfying the constraints.
3.	The authors suggest restricting $f$ to the tangent space to satisfy the constraints. This leads to a question related to the one above: Is the original problem Eq. (1) with constraints on the state equivalent to the constrained problem Eq. (8) considered by the authors with the forcing term in the tangent space? Eq. (8) may limit the range of expressible dynamics.
4.	There is a concern that the computation of the adjoint differential and pseudoinverse in Eq. (7) would be quite difficult for general $g$.
5.	Instead of enforcing hard constraints, constraints could be incorporated into the loss function by penalizing it. For instance, this could involve adding the $L^2$ norm of $g(u)$=$(g(u)-0)$ as a regularization term to the existing NODE loss. An experimental comparison with this approach seems necessary.
6.	While hard constraints ensure that the constraints are satisfied, they are not necessarily superior to soft constraints. Hard constraints can limit the representational power of the network and may negatively impact training because of their complex computational structure. It is crucial to understand and experimentally verify the trade-off between satisfying constraints and the model's capacity.

---

confidence: 4
contribution: 3
initial_rating: 5
presentation: 4
questions: 1. Could the paper provide a detailed formula of computing the "Relative State Error" and "Constraint Error"?

2. Could the paper explain in detail how the trajectory of the figure was selected?

3. What does the vertical axis of the leftmost subgraph in Figure 2 mean？

4. Section 4.2: Why does NDE using generalized coordinates that satisfy constraints perform poorly, and does this mean that preserving constraints cannot directly indicate better predictions? Can the results highlight the importance of constraints?

5. Section 4.3: Do we know the governing function for this example, or are the dynamics learned from the given data? Additionally, the statement 'apply random perturbations to each grid (see Appendix A) and learn the dynamics from the resulting transients' feels unclear. Could the paper clarify this? Also, what are the sizes of the training and test datasets used for this example?

6. For all examples, the paper only presents the state error for a single test trajectory. Could the paper provide more comprehensive quantitative evaluations? 

7. The assumption of known constraints may be too restrictive. It would be helpful to discuss scenarios where the constraints are known, but the governing function is unknown and data is provided.
soundness: 2
strengths: 1. The paper presents an approach to integrating known constraints into neural differential equations, which is a relatively unexplored area in the field of machine learning and dynamical systems. The introduction of the projection method to enforce constraints on the vector field is innovative.

2. The empirical results demonstrate that PNDEs outperform existing methods in terms of accuracy and stability. The experiments conducted on challenging examples, including chaotic systems and power grid models, further validate the robustness of the proposed method.

3. The paper is well structured and clearly written, making complex concepts accessible to readers.
summary: The paper introduces Projected Neural Differential Equations (PNDEs), a novel method designed to incorporate known constraints into neural differential equations. Specifically, PNDEs project the vectors of the vector field onto the tangent spaces of the manifold defined by known constraints. This allows for the enforcement of various constraints without the need for specific coordinate systems. The paper provides empirical evidence through experiments showing that PNDEs satisfy constraints more accurately than state-of-the-art baselines.
weaknesses: 1. There is a potential ambiguity in the notation used in the paper. Specifically, the definition of $\mathcal{E}$ is lacking; is it $\mathbb{R}^n$? The notations $f_{\theta}$ and $\bar{f}_{\theta}$ in equations 1, 3, and 8 are not entirely consistent.

2. The paper primarily focuses on known constraints. In many cases—such as the first example presented—it shows that if the constraints are known, the corresponding total differential equation can be determined.

3. The effect of the ODE solver is not discussed in this paper. Ideally, the constraints should be preserved exactly, as stated by Proposition 1 in the paper. However, the numerical results indicate that they are only preserved approximately, albeit with small errors. These discrepancies may arise from numerical errors introduced by the ODE solver.

4. The setup of the first example differs from that in Section 2. Is it possible to explicitly write out the manifold for the example?

5. The empirical results demonstrate that the proposed PNDEs outperform existing methods in terms of the consistent error. However, the improvement in terms of prediction error is less pronounced. I am not certain whether the ability to preserve the given constraints is the most critical indicator. If we aim to improve this, the simplest and most straightforward approach would be to project the predicted state onto the known manifold during the prediction phase.