Stable Port-Hamiltonian Neural Networks
Authors: Fabian J. Roth
Dominik K. Klein
Maximilian Kannapinn
Jan Peters
Oliver Weeger
Venue: NeurIPS 2025 poster ()
Keywords: dynamic system identification
port-Hamiltonian systems
Lyapunov stability
Primary area: machine_learning_for_sciences
TL;DR: The present work proposes stable port-Hamiltonian neural networks for accurate, robust, and reliable identification of stable nonlinear dynamic systems.

Abstract: In recent years, nonlinear dynamic system identification using artificial neural networks has garnered attention due to its broad potential applications across science and engineering. 
However, purely data-driven approaches often struggle with extrapolation and may yield physically implausible forecasts. 
Furthermore, the learned dynamics can exhibit instabilities, making it difficult to apply such models safely and robustly. 
This article introduces stable port-Hamiltonian neural networks, a machine learning architecture that incorporates physical biases of energy conservation and dissipation while ensuring global Lyapunov stability of the learned dynamics. 
Through illustrative and real-world examples, we demonstrate that these strong inductive biases facilitate robust learning of stable dynamics from sparse data, while avoiding instability and surpassing purely data-driven approaches in accuracy and physically meaningful generalization. 
Furthermore, the model's applicability and potential for data-driven surrogate modeling are showcased on multi-physics simulation data.

Stable
Neural Networks
Fabian J. Roth, Dominik K. Klein, Maximilian Kannapinn, Jan Peters, Oliver Weeger Technical University of Darmstadt, Darmstadt, Germany kroth,klein, kannapinn, weeger &ecp tu-darmstadd .ee,, peters0ias. cereeiias.ttudrrree tu-darmstadt.c. de
Abstract
In recent years, nonlinear dynamic system identification using artificial neural networks has garnered attention due to its broad potential applications across science and engineering. However, purely data-driven approaches often struggle with extrapolation and may yield physically implausible forecasts. Furthermore, the learned dynamics can exhibit instabilities, making it difficult to apply such models safely and robustly. This article introduces stable nort-Hamiltonian neural networks, a machine learning architecture that incorporates physical biases energy conservation and dissipation while ensuring global Lyapunov stability of the learned dynamics. Through illustrative and real-world examples, we demonstrate that these strong inductive biases facilitate robust learning of stable dynamics from sparse data, while avoiding instability and surpassing purely data-driven approaches in accuracy and physically meaningful generalization. Furthermore, the model's ss applicability and potential for data-driven surrogate modeling are showcased on multiphysics simulation data.
Introduction
Recent years have seen strongly growing interest in using machine learning (ML) identify the dynamics of physical systems directly from observations [5, 28, 31, 32, 40, 42, 53]. This approach has diverse applications, ranging from modeling phenomena with partially or entirely unknown underlying physics to creating surrogate and reduced-orda models for rapid simulations and digital twins [23, 53, 54]. In the latter application, data-driven methods can address the drawbacks of traditional numerical methods, which often require tremendous computational resources and expertise [1, 53J. However, in turn, classical ML algorithms typically require large amounts training data and lack robustness, as they often suffer from poor extrapolation and are prone to making physically implausible forecasts [1, 13, 53, 54). To overcome these limitations, a general consensus is to avoid discarding centuries" worth of established knowledge in physics, which is implicit consequence of relying solely data-driven methods. Instead, there is a growing emphasis on incorporating physical priors and inductive biases into ML approaches [33, 54]. This physics-guice guided ML promises improved accuracy, robustness, and interpretability while reducing data requirements [33, 53).
Several physics-guided architectures have been proposed for learning dynamical systems. These approaches often parameterize an ordinary differential equation (ODE) identify continuous-time dynamics. To ensure adherence the laws of thermodynamics, many architectures build upon frameworks such as Hamiltonian mechanics [9, 16, 20], Lagrangian mechanies [10, 34], Poisson systems [19], GENERIC [18, 56), port-Hamillonian systems (PHS) [11, 13, 37, 38, 57), or the generalized Onsager principle [55]. The latter three allow for the modeling of dissipative dynamics, while the former are generally limited to mnergy-conserve
systems. Additionally, PHS explicitly incorporate external excitations, making them particularly suited for modeling controlled systems.
The mentioned physical modeling frameworks have implications for the stability of the learned dynamics. For example, Hamiltonian neural networks, stable equilibria correspond extrema in
39th Conference on Neural Information Processing Systems (NcurIPS 2025).

---

the Hamiltonian energy function. However, these effects are often not discussed in the referenced works. Nevertheless, stability is an important and fundamental property of any dynamic system. Guarantees on the stability of the dynamics may improve the model's 's robustness, help trust their predictions, and can themselves be viewed as physically motivated biases [14].
Given the importance of the stability property, multiple other works have proposed methods to enhance the stability of learned dynamics, often using neural ordinary differential equations (NODEs) [7] as a basis. Massaroli et al. [36] propose a provably stable NODE variant that guarantees the existence of locally asymptotically stable equilibria, Multiple studies have combined Lyapunov's 's theory with neural networks (NNs) to guarantee the stability of discrete-ting dynamics [14, 29] or to estimate the region of attraction for given dynamical systems [3, 44]. Kolter and Manek [27] describe another Lyapunov-stable architecture based on NODEs. They concurrently learn unconstrained dynamics with a NODE and a convex Lyapunov function. Global stability around an equilibrium point is then ensured by projecting the NODE-dynamics onto the space of stable dynamies as given by the Lyapunov function. Takeishi and Kawahara [50] extend this approach from stable equilibrium to stable invariant sets. Similar projection-base based techniques have also been proposed to learn input-output stable dynamics [26] and dissipative systems [39]. However, these projection methods can introduce discontinuities, hindering training strategies that rely on integrating trajectories rather than matching derivvevivv pairs [46]. Other approaches employ physical modeling frameworks to formulate the learning problem directly in the subspace of stable dynamics. Rettberg et al. [43] propose identifying a linear latent PHS, which makes enforcing global stability relatively straightforward. Yu et al. [55] propose model based on the generalized Onsager principle that can guarantee bounded solutions and the existence of locally asymptotically stable equilibria. While their evolution equation is similar to PHS, their approach does not consider arbitrary, time-dependent external input signals.
This work exploits the connections between stability and physical frameworks and proposes stable port-Hamiltonian neural network (sPHNN) for learning nonlinear and globally stable dynamical systems. The method is projection-free and can incorporate input signals. We show that the approach is capable of learning equilibrium points, though our numerical evaluations focus on practical applications where the equilibrium is known a priori. In such cases, this prior knowledge can be seamlessly integrated into the proposed architecture, improving prediction accuracy while simultaneously reducing data requirements.
2 Background
2.1 Stability of dynamical systems
We consider autonomous dynamical systems first-order ODEs, where the state x(t) E R" evolves according to a vector-valued function f RT + R":
i (t)) = f(x()))
(I)
Suppose fi is sufficiently smooth to ensure unique solution æ(f) for any initial condition æ(to). A solution (tt) is stable in the sense of Lyapunov if, for any small perturbation in the initial condition, the perturbed solution æ(t) remains close for all times, that is, if for any > 0, there exists 5 0
such that (to) - =(to)ll < implies [(t(t) - =(t)l < c for all t >> to [35, 52]. A stronger version of stability, namely asymptotic stability, is obtained if additionally the perturbed solution (t)) converges to á() for ( + do for initial conditions ([)(()) #(to)] < A with some A > 0 [52].
Lyapunov stability is local property, in the sense that for sufficiently large perturbations, the per- turbed solution may stray arbitrarily far even from asymptotically stable solutions [49]. Quantifying the permissible extent of perturbations on the initial condition such that the perturbed solution still converges to the unperturbed solution leads to the concept of a basin of attraction [30]. Asymptoti- cally stable equilibria, for which the basin of attraction extends to the entire phase space, are called globally asymptotically stable [30]. A system with a globally stable equilibrium cannot have any other equilibria, as their existence would imply that there solutions (namely, the other equilibria) that do not converge to the globally stable equilibrium.
In the following, we consider the case of an equilibrium solution ==4) = 0, i.e., f(0) = 0. Placing
the equilibrium at the origin is done for ease of notation and without loss of generality, as arbitrary equilibria can be translated to the origin with dime-independent transformation. The stability of the

---

equilibrium can be shown through the use of Lyapunov's stability theory by finding a suitable scalar Lyapunov function V [24, 52]. Let V: D + R. be continuously differentiable and positive definite, that is, (()) = 0 and ((x) > 0 in D (0), where D C IR"h is a neighborhood of the origin. If the value of V is nonincreasing along trajectories of the system D, then the equilibrium is stable the sense of Lyapunov. Mathematically, this is the case if (ax) is negative
Tv 0V
V(x)=
dæ x dx
f(x) S
Vx D (0).
(2)
In the case of a negative definiteness (v(x) < in D {0)), asymptotic stability can be concluded.
To obtain global asymptotic stability, these conditions must be satisfied globally, that is, for all
E D = IR.", and the Lyapunov function must be radially unbounded, i.e., V(x) - oo for | -
oo [24]. The latter condition avoids the existence of solutions that do not converge to the equilibrium by diverging to infinity without violating the decrease condition in Equation (2).
2.2 Port-Hamiltone systems
PHS theory provides a framework for modeling physical systems that generalizes the underlying mathematical structure of Hamiltonian dynamics. It incorporates the modeling of energy-diss beeggy-disippating elements, which are typically absent classical Hamiltonian systems [51]. Additionally, PHS theory has roots in control theory, emphasizing the interactions of dynamic systems with their environment. In this work, we consider PHS of the form
0H. Φ = (1)) - R(x)]
(x) + G(z))(())
dx
Here, x(t)= E IR' describes the system's 's state, u(t E R. its input, and H. : R" + E denotes the Hamiltonian, which typically represents the total energy of the system [4, 51]. The skew-symmetric structure matrix J : R -> R (x) (xx) describes the conservative energy flux within the system. The symmetric positive semi-definite matrix R: R R R(x:) (():) R (a) > 0 is called dissipation matrix and describes the energy losses. Finally, the input matrix G R - 1
describes how energy enters and exits the system via the inputs. PHS extend the energy conservation of Hamiltonian systems to the energy dissipation inequality:
TH. OH
THL
OH.
THe
H. =
R
+
Gu(t) S s(x,, u).
dæ
dx
dx
8x
(4)
20
s(a,,u)
In the language of system theory, the system is dissipative with respect to the supply rate s(x, u). The formal definition of dissipativity requires H 2 O, which will be ensured later. Since the supply rate s(x, u) is linear in the input it vanishes for zero-input s(x, 0) O. The inequality thus ensures that the system cannot gain energy in the absence of an external excitation (u = O). Besides this physical motivation, the dissipativity of PHS plays crucial role in stability analysis. For the unforeed system with O, the Hamiltonian inherently fulfills the decrease condition H < O. Thus, by using H as a Lyapunov function, it can be shown that any strict minimum in the Hamiltonian implies a stable equilibrium of the unforced dynamics [51]. Furthermore, the equilibrium is asymptotically stable if the system continuously dissipates energy (R > 0) for states near the minimum.
The use of the nort-Hamiltonian formulation as basis for an ML model has multiple favorable consequences. When learning system dynamics from data, the structure of Equation (3) provides both a physical bias and the opportunity for interpreting the resulting trained model, while enabling the modeling of a wide range of relevant systems. Equation (4) highlights the distinctive roles played by the three matrix-values functions. The structure matrix J (xx) results in conservative dynamics, the dissipation matrix R(x) contributes only to the dissipative dynamics, and the input matrix G(x) captures the energy flow over the system boundaries due to the inputs. This knowledge improves the interpretability of the trained model components, Furthermore, it allows new ways of incorporating prior knowledge into the system. For example, the learned system can be constrained to be conservative by fixing R = O. Furthermore, the connection of PHS to the decrease condition from Lyapunov's theory provides the opportunity to introduce further desirable stability constraints.
To illustrate the above concepts with a simple example, consider the dynamics a three-dimensional spinning rigid body described by Euler's 's rotation equations. Expressing the angular velocity vector
3

---

Bo
E NODE ESNODE
Energy
E PHNNN OPHNNN
E aPHNN-LM
200
400
200
400
Time t/8
Time 8/8
Figure Spinning rigid body: Interquartile mean (lines) and range (shaded regions) the energy computed from the predicted states. Left: Models without stability bias; Right: Models with stability bias.
(t) E R in the principal-axis coordinate frame, the evolution equation given by
Iw + w x (Iw) = pp Vt O.
(5)
Here, I d1ag(4,, 12, 13) denotes the inertia matrix and μ is a damping coefficient. By introducing
1 H(w) 2 Fl 12wي + 2 13w3.
(6)
which describes the kinetic energy the system can be rewritten in the form of Equation (3), with state-dependent structure matrix J, constant dissipation matrix R, and vanishing input matrix G:
La Ws
0
W1
112
Ish.
W2
OH
W2
Ia
W3
1213
w1
1
0
72
0
ow
W
W2
Tals W1
0
1 75
J(W)
In anticipation of the numerical experiment in Section 4.1, Figure 1 shows the predicted energy E from various models trained on trajectories from the spinning body system. Those incorporating both, physical biases from the PHS formuulamnt and stability (bPHNN, sPHNN, and sPHNN-LM) perform best, highlighting the importance of these constraints for learning system dynamics.
3 Stable port-Hamiltonian neural networks
Our goal is learn the dynamics of physical systems in the form of Equation (3) with desirable stability properties from observations. To this end, we parameterize the components H. J, R and G of PHS with NNs.
3.1 Stability of Hant-Hamiton dynamics
To guarantee the existence of stable equilibria of port-Hamiltonian systems, is sufficient to enforce the existence of strict local minima in the Hamiltonian. However, these stability implications only local, meaning the size and shape of the corresponding basin attraction are not constrained by Equation (3). Small perturbations might be enough to stray far from the stable equilibria. Furthermore, the general PHS can represent dynamics such as unbounded trajectories, which must be considered unstable from a practical standpoint. The following describes the necessary requirements for achieving global stability. Essentially, the Hamiltonian H. is required to be a suitable Lyapunov function, which includes positive definiteness and radial unboundedness. These properties are achieved by constraining the Hamiltonian to be convex and ensuring the existence of strict minimum. This approach is similar the construction of Lyapunov function by Kolter and Manek [27], but differs in the way the minimum is enforced.
Theorem 3.1. Consider the PHS Equation (3) in the unforced case ul)
0H dæ (x),
== (x)) R(x))
with J J=
R 0.
(8)
Suppose the Hamiltonian H(xx) is convex, twice continuously differentiable, and fulfills:
0H
8 x 0.
H(0)= = 0,
= 0.
(9)
dæ ==0
dx ==

---

Then, H. is suitable Lyapunov function for showing stability of the equilibrium at xxt) = 0, and all solutions are bounded. Furthermore, the equilibrium is globally asymptotically stable if R(
A proof of Theorem 1.1 is provided in Appendix A. However, its implications also follow clear physical intuition: If a dynamical system continuously loses energy via internal dissipation (R > 0), will eventually settle into energy minimum corresponding an asymptotically stable equilibrium. The convexity requirement, together with Equation (9), guarantees the existence of unique strict and global minimum at the origin. Additionally, ensures that the energy for states infinity is unbounded, i.e., H(x) + 30 as ||||| + oo. Thus, no finite state has more energy than state at infinity, and all trajectories must be bounded. As result, all solutions eventually converge the global energy minimum. While convexity is sufficient to ensure global asymptotic stability, it is not necessary condition. To enhance model expressiveness when needed, techniques such as input warping [27] can be employed, but were not required in the applications considered in this work.
For conservative systems with R 0, energy is not dissipated, and attractive equilibria cannot exist. Therefore, a conservative system cannot possess asymptotically stable equilibria. Nevertheless, the minimum in the energy still represents stable equilibrium, and all trajectories stay bounded. The mixed case, where ł O, only allows the same conclusions as the conservative case. However, unless the dissipation matrix vanishes in a large region in phase space, it is unlikely that system with R ž 0 is not globally asymptotically stable.
3.2 Neural network model architecture
Having established the theoretical requirements for stability of PHS, we now outline how these conditions be met with dedicated NN model. We call the resulting architecture sPHNN. Its computation graph is illustrated in Figure 2a. Figure 2b shows the decomposition of a sPHNN with randomly initialized NNs into conservative and dissipative dynamics, highlighting the model's physical bias and inherent interpretability.
Hamiltonian To fulfill the constraints of Theorem 3.1, 1, the Hamiltonian H(xx) is represented by fully input convex neural network (FICNN) [2]. FICNNS can approximate any convex function [8]. They are defined by the recurrence relation
do (Wox bu), Zi+1 σ; (Uiz; + W. + by). i 1,2,
(10)
f(x) Zk:
with the pass-through weight matrices W, that map the input æ directly the activation of the k: hidden layers, bias vectors b;, and layer outputs 2,. If the weight matrices U; have non-negative
Input
Purameters
Operations
Output
J J8H/0a
Conservative
3
Total dynamics
f(z))
from (an, x²)
H
FICNN
R) UM
Gu
R ROHlda
LLA
3
G
Dissipative
FENN reshape
(a) Computation graph.
(b) Composition of the total dynamics.
Figure 2. sPHNN model architecture: (a) Computation graph sPHNN. The FICNN parameterizing normalized to obtain H. and the outputs of the FENNS for J. L, and G are reshaped to be skew-symmetric. lower triangular, and rectangular matrices, respectively. (b) Dynamics of randomly initialized sPINN for = 2, showing the built-in interpretability through the separation into conservative and dissipative dynamics.
5

---

components, the first activation function do is convex, and all subsequent o, are convex and non- decreasing. then the output f(z) is convex in the input æ. Since the Hamiltonian should be twice continuously differentiable, the same must hold for the activation functions. One suitable choice is the softplus activation. The conditions for stability in Equation (9) are then fulfilled by normalizing the FICNN f:
0 dx
H(()) f((=) -f(x))
(x x
+ I x."
(11)
a0
fwom (x).a=)
Fing (an
Consequently, the Hamiltonian H. and its gradient vanish at x' due to the normalization term fnom- Moreover, the regularization term frg with an arbitrarily small c >> O ensures positive definiteness of the Hamiltonian's 's Hessian. However, in practice, Jry can be omitted since the FICNN provides enough bias towards local strict convexity. After training. the stability guarantee can then be recovered by verifying the Hamiltonian's positive definiteness at x". Choosing x 0 fulfills all requirements of Theorem 3.1. However, the presented approach is more general. As the assumption that the equilibrium is located at the origin is only made for ease of notation, the normalization approach in Equation (11) can position the stable equilibrium anywhere in phase space and is not restricted to the origin. If the true equilibrium position is known beforchand, fixing x' introduces an additional inductive bias. While this information can often be obtained from prior physical knowledge, it may not always be available. In such cases, the equilibrium can be inferred directly from data by treating x² as trainable parameter during optimization.
Structure, dissipation and input matrices The matrix- functions J (x), R(x) and G(x)) are parameterized by feedforward neural networks (FFNNs), where the outputs are suitably reshaped into matrices. As constrained to be skew-symmetric, the respective FFNN's output vector mapped to the space of skew-symmetrie matrices via a linear transformation. Similarly, the positive definiteness (or semi-definiteness) of R is ensured via the Cholesky factorization R LLT, where the reshaped FFNN output L is a lower triangular matrix with positive (or non-negative) elements along the main diagonal. For some applications. J. R and G are independent of æ, and thus, learning constant matrices is sufficient. In particular, this is the case when the conservative dynamics of the system under consideration are Hamiltonian in the given coordinates. Furthermore, by fixing J to the symplectic matrix (and choosing R G 0), the sPHNN can be reduced to Hamiltonian NN.
3.3 Training stable port- port-Hamiltonian neural networks
There are fundamentally different approaches for training time-continuous dynamic models such as NODES and sPHNNs: derivative and trajectory filting. The former directly compares the predicted state derivatives à from Equation (3) the true ones. While this approach is efficient, it requires the true time derivatives to be available as training data, which is typically not the case for real-world measurement data. Furthermore, derivative fitting is not applicable when augmented states are used [12], i.c., when x is padded with additional dimensions compared to the observable states to enrich the representable dynamics. Further, model proficient in predicting á does not necessarily produce accurate predictions for trajectories ztt). as small derivative errors can accumulate exponentially when integrating an ODE [17]. On the other hand, trajectory fitting directly compares integrated model trajectories az(4 with ground-truth trajectories, ensuring long-term accuracy. Here, the optimization process necessitates either the propagation of the loss gradients through the ODE solver or the use of the adjoint sensitivity method [7], which computes the gradients by solving modified ODE. In both cases, trajectory fitting incurs significant computational costs due to gradient propagation and model integration, making trajectory fitting much slower than derivative fitting.
distinct advantage of sPHNNs is the favorable influence the stability constraint exerts on the training process when using trajectory fitting. NODEs and similar architectures can suffer from exploding gradients, where small parameter changes lead to large loss gradients. This can slow the training process or prevent it from converging. A potential cause is the crossing between the boundaries of local basins of attraction during training, leading to significant differences in long-term behavior [41]. Since sPHNNs have only one global basin of attraction, this source of training difficulties is eliminated. In our numerical experiments, we observe good convergence of the training error of sPHNNs where NODEs struggle without the use of regularization terms gradient clipping.

---

4 Experiments
We evaluate the proposed model on multiple example problems, including measurement data from a real physical system and multiphysics simulation data. The performance of sPHNNs is compared to that of NODEs, PHNNs, and bounded PHNNs (bPHNNs). The latter two are PHS-based NN models that do not enforce or relax the convexity constraint of sPHNNs, and thus do not guarantee global stability. In PHNNS, the Hamiltonian is modeled by an unconstrained FFNN, whereas in bPHNNs, it is constructed to be positive and radially unbounded. which ensures bounded trajectories. The latter architecture is intentionally designed to resemble OnsagerNet [55], but with the key difference that it supports arbitrary dimendependent inputs, making it suitable for our experimental setting. For examples that purely rely on derivative fitting, we further include stable NODEs (sNODEs) [27]. Using trajectory fitting, we could not achieve convergence of sNODEs, excluding them from most experiments. This is likely due to difficulties caused by discontinuous dynamics stemming from the baased stability enforcement, also mentioned by [46]. Finally, we include variant of sPHNNs with learnable minima (sPHNN-LM), In contrast to s sPHNNs, which have x; fixed to an a-priori known equilibrium position, APHNN-LMs have to infer the equilibrium position from the training data and start out with randomly initialized x". All models are trained by minimizing the mean squared error using the ADAM optimizer [25]. Model hyperparameters and computational costs are listed in Appendix B. Predictions are generated by rollout with numerical integrator using the Runge-Ko Kutta scheme Tsit5, with adaptive step size and known input function u(t). Code and model weights are available on GitHub'.
4.1 Spinning rigid body
We generate data of damped spinning rigid body by numerically integrating Euler's rotation equations, which were already introduced in Section 2.2, Equation (5). Here, accurately capturing the dynamics of the angular velocity vector w with PHS requires state-dependent structure matrix (zx) We train ten instances per model type using derivative fitting with (w, w) pairs sampled from the trajectories. Figure shows the rotational energy F!(w) computed from model predictions w via Equation (6). The PHNN, NODE, and sNODE models show deviations from ground truth, with the former two not converging to zero. Furthermore, the NODEs show unphysical behavior with at times increasing energy and exhibit high variance between instances. In contrast, sPHNNs and bPHNNs match the true energy well, highlighting the importance of providing both stability and physical biases. The sPHNN-LM variant recovered the actual equilibrium position accurately, with an interquartile mean distance between the true and inferred minimum across all instances of 0.010. As result, its predictive performance is nearly identical to that of the sPHNN with fixed x".
4.2 Cascaded tanks
Next, the sPHNN model is evaluated using the cascaded tanks dataset [47, 48). It contains measure- ments of a physical fluid level control system consisting of a pump and two tanks, see Appendix D for details. The training data consists of single trajectory with 1024 pairs of pump voltages u and water level measurements of the lower tank. We model the system using vector x [21, T, where y = x2. Assuming the augmented state 21 describes the water level in the upper tank, it is straightforward to determine the equilibrium of the autonomous system. With zero input, the tanks will eventually drain; therefore. the stable equilibrium is given by x We thus fix the minimum of the sPHNN's Hamiltonian to the origin. To investigate the effects of a potential lack of this knowledge, we also train the sPHNN-LM variant. Its minimum location x is randomly initialized and optimized during training. Apart from this, both sPHNN versions are identical and use a constant fixed symplectic matrix as the structure matrix, as well as constant but learnable dissipation and input matrices. The same structure applies to PHNN and bPHNN.
wwo-dimensional state
We train 20 instances per model type and show a statistical evaluation of the root mean squared eтrors (RMSEs) in Figure 3a. Overall. sPHNN, sPHNN-LM, and bPHNN achieve comparable accuracy, while NODE and PHNN perform worse on average across both trajectories. In addition to their better performance, the models with stability bias exhibit lower variance across model instances, indicating more consistent results. To evaluate the models' extrapolation capabilities, we extend the evaluation data with 400 of zero-input The resulting predictions are depicted in
7

---

sPHNN
sPHNN-LM
OPHNN
PHNN
NODE
10
Ground Truch
Ground Truhh
PHNN
sPHNN &PENNN-M SPINN
NODE
0.8
RMSE/V
20/ P//V
0.6
0.4
0.2
Train Test
Train Test
Train Text
Trin Test
Train Test
3,500
4,096
3,500
4,096 Time f/5
Time t/a
(a) Training and test RMSE.
(b) Relaxation.
Figure 3. Cascaded tanks: (a) RMSE of the models on training and test trajectory. (b): Predictions for the extended test trajectory. At 4096s, s. the pump is turned off. Lines correspond to the interquartile mean and shaded areas represent the interquartile range of the predictions from the 20 model instances.
Figure 3b. Assuming perfectly calibrated sensors, the true system response is an eventual zero output as both tanks drain. By construction, this behavior is guaranteed for sPHNNs. The 20 sPHNN-LM instances reliably identified the systems equilibrium to lie between -0.196 V to 0.432 V, which aligns with physical intuition. In contrast the 6PHNN, NODE and PHNN models exhibit poor physical extrapolation, predicting equilibria far from zero, diverging water levels, or oscillatory behavior. In an additional experiment presented in Appendix D, the sPHNN models are also shown to have successfully learned the hard saturation effect that occurs when a water tank overflows task at which both NODE and PHNN fail.
4.3 Thermal food processing surrogate
With the applicability of sPHNNs real-world data established, this section focuses on one of the promising applications of data-driven dynamic system identification: the efficient construction of surrogate models. Here, data from conjugate heat transfer simulation of convection oven coupled with soft matter model for meats from Kannapinn et al. [21, 23] is employed. The data consists of trajectories representing the temperature histories TA and TB at two probe points within the meat. These result from predefined excitation signal Towen which controls the oven temperature. The surrogate's 's task predict the probe temperatures given the oven temperature as input.
We fit the models to initially no = trajectories consisting of 280 samples each (see Appendix E). Per model, we vary the dimensionality of the state Rn with 2 + nA by adding 0 3 augmented dimensions YI.A The minimum x of the sPHNNs Hamiltonian is set to the point in phase space corresponding to the ambient temperature, motivated by physical intuition about the stable thermal equilibrium. The initial values of the augmented states are fixed to zero. In total, 20 instances are trained per model type and augmented dimension.
The RMSEs evaluated on 15 test trajectories (see Appendix E) are depicted in Figure 4a. The sPHNN predictions tend to improve with more augmented dimensions, achieving errors comparable to those obtained by Kannapinn et al. [21] using commercial software. This appears intuitive, as additional hidden state variables may be required accurately describe the evolution Tf and Tр that results from the multiphysics, conjugate heat transfer process, see also [12]. However, for bPHNNs and sPHNN-LMs no clear trend emerges. The NODEs and PHNNs even tend to deteriorate in performance with more augmented dimensions, likely due to instabilities in the learned dynamics. This becomes apparent Figure 4c (top row), which depicts the models' predictions for extended test case. While the sPHNNs consistently converge to temperatures close to thermal equilibrium, the predictions of NODES and PHNNs tend become unstable when one or more augmented dimensions are used. Although all models improve with more training data (Figure 4b), these instabilities persist when extrapolating in time (Figure 4c, bottom row > 1395 s). In contrast, the stability constraint of the sPHNNs and sPHNN-LMS allows safe use of the additional flexibility provided by a higher- dimensional state space, avoiding the pitfalls observed in the unconstrained models. Given sufficient training data, the sPHNN-LM variant performs on par with the sPHNN, whereas in the scarce-date regime its errors increase and are comparable to those of the bPHNNS. Although the bPHNNs also remain stable, they are consistently outperformed by the sPHNNs for 3. regardless of dataset size, and especially when only a single training trajectory available up 1, see Figure 4b and Figure 4c (bottom left).
8

---

10a
SPINK MPHNN
RMSE/KK
RMSE/K
10!
PEINN
109
10
Number of augmented dimensions 19
Nomber training trajectories R
(a) Test RMSEs for un 2 training trajectories.
(b) Test RMSES for nr 3 augmented dimensions.
2
3
400
Tow
&PHNN &PHNN-LM
350
1D
PIINN
300
NODE
nD 1
TD
YED 15
100
350
71,
300
1,395 3,000 Time
1,395 3,000 Time f/e
1,395
Time l/s
Time t/s
(c) Model predictions for different combinations augmented dimensions and training data.
Figure 4. Thermal food processing surrogate: (a) and (b) RMSEs evaluated on 15 test trajectories for various numbers of augmented dimensions and training trajectories np. (c) Interquartile mean and range of the Ta predictions for custom test case. Top row: Varying number augmented dimensions with fixed number nn 2 of training trajectories. Bottom row: Varying number of training trajectories with fixed number n. = 3 of augmented dimensions. The time = 1395 s marks the length the training trajectories.
Table 1. Median RMSE of 20 instances per models using n, = augmented dimensions, trained on np = noisy trajectories and evaluated on noise-free training and test sets.
Training
Test
Model
% noise
15 % noise
25 % noise
5 noise
15 % noise
25 % noise
sPHNN
0.332
0.854
1.354
1.170
2.595
3.177
SPHNN-LM
0.638
1.096
1.622
2,420
2,888
3.677
bPHNN
0.500
1.049
2.207
1.572
2.828
5.231
PHNN
2,828
6,750
5.779
28.233
25.631
24.224
NODE
0.993
1.974
2.423
3.577
3.482
4.362
To assess the models' robustness to noise, we conducted an additional experiment using n, augmented dimensions and n 2 training trajectories, introducing zero-mean Gaussian noise for both the input Toven and the outputs T'A and T3. The noise amplitudes were varied between % and 25% % of the original signals' standard deviations. Table 1 reports the average RMSE over 20 model instances per model type, evaluated on noise-free test data. To isolate the effect of noise and eliminate generalization error, also evaluated the RMSE on the noise-free free training data. All models demonstrate a degree of robustness to noise, since the integration step involved in trajectory fitting acts low-pass filter, effectively smoothing the learned dynamics. As result, high-frequency Gaussian noise is attenuated, which helps prevent overfitting and improves generalization. Nevertheless, the sPHNN consistently outperforms the baselines on both the training and test sets.
4.4 Additive manufacturing surrogate
Finally, we evaluate the applicability of sPHNNs to higher-dimensional state spaces and build reduced order model (ROM) surrogate of 3D partial differential equation problem. To this end, generate field data by numerically solving the heat conduction equation for moving heat

---

1 0.
t= 2.58
1 4.0s
1 19.0s
-1040
Source
1,000
1.7
True
SPHNN
T/K
sPHNN-LM
650
0/N0
bPHNN
PHNN
0000000
NODE
Figure 5. Heat source field and temperature field predictions on the cuboid domain for a test case with 12.5 mms Q = 400 W. The instances selected for this evaluation resulted in the median test error for the corresponding model type. Colors clipped remain in the legend's range.
source on a cuboid domain with convective and radiative thermal boundary conditions. The finite element simulation is designed to model the evolution of the temperature field of metal additive manufacturing process, see [22] for details. total, 25 trajectories are obtained by varying the heat source speed from 10 mmss to 20mm and power Q from 300 W 500 W. All trajectories span 20 s, with a non zero heat source lasting for the initial 2.5 s to 6 depending on v. However, directly using the discretized field data would lead very high-dimensional models and render the training and inference inefficient. We thus apply a proper orthogonal decomposition (POD) to map the temperature and heat source fields onto 40-dimensional latent spaces, respectively (see Appendix F for details). The models are trained on the latent representations of the two trajectories with (v, Q) (10 mms s1 300 W) and (20 mm s 500 W) using derivative fitting and subsequent fine-tuning with trajectory fitting.
Figure 5 depicts the heat source field (top row) and true and predicted temperature fields from the different model types. While the sPHNN's total test RMSE peaks below 30 K (see Appendix F, Figure F.2) and decreases as the system approaches thermal equilibrium, the errors for NODE and PHNN increase rapidly, driven by instabilities the learned dynamics. These instabilities are likely caused by the sparsity of the training data, which leaves gaps in the phase space, allowing unstable dynamics to emerge. The stability constraint of sPHNN and sPHNN-LM mitigates these issues, ensuring robust performance even with scarce training data, while outperforming the bPHNN thanks to their stronger inductive bias.
5 Conclusion
This work proposes stable port-Hamiltoiial neural networks (sPHNNs) for accurate, robust, and reliable identification of stable nonlinear dynamics. They leverage the nort-Hamiltonian framework, which provides a basis for model interpretation and ensures adherence to the laws of thermodynamics in terms of energy conservation dissipation balances with respect learned Hamiltonian energy function. Additionally, the approach ensures global asymptotic stability of the identified dynamics, without requiring projection, by constraining the Hamiltonian be convex, positive definite Lyapunov function. While this limits sPHNNs to modeling systems that appear globally stable within the region of interest, provides a strong inductive bias that enables robust learning and reasonable generalization from sparse data. We demonstrate the approach's 's viability using data from synthetic examples, real-world measurements, and complex multiphysics simulation models. The resulting models, built with small neural network architectures, exhibit low variance across instances, accurate prediction, and robust generalization. Furthermore, the proposed method can successfully exploit augmented dimensions where unconstrained alternative models struggle with instabilities or data scarcity. While the approach can identify the equilibrium point from data, we find that providing the equilibrium as prior can further enhance prediction accuracy and reduce data requirements. Future work will extend the approach to multistable systems by allowing Hamiltonians with multiple minima, while retaining favorable stability properties.
10

---

Acknowledgments and Disclosure of Funding
The authors acknowledge the financial support provided by the Deutsche Forschungsgemeinschatt (DFG, German Research Foundation, research grant number 4927701 17 and the CRC/TRR 361 CREATOR), the Hessian Center for Artificial Intelligence (Hessian.AI, project EnDyRo), and the Graduate School Computational Engineering at TU Darmstadt. Furthermore, we would like to thank An Thai Le, Daniel Palenicek, and Joe Watson from TU Darmstadt (IAS) for their constructive feedback on the manuscript.
References
[1] A. Abbasi, P. N. Kambali, P. Shahidi, and C. Nataraj. "Physics-Infonnn Machine Learning for Modeling Multidimensional Dynamics". In: Nonlinear Dynamics (2024). DOI: 10. 1007/ 10001-002--0063--3
[2] B. Amos, L.. Xu, and Z. Kolter. "Input Convex Neural Networks". In: Proceedings the 34th International Conference on Machine Learning. Vol. 70. Proceedings of Machine Learning
Research. PMLR, 2017, pp. 146-155.
[3] M. Barreau and N. Bastianello. "Learning and Verifying Maximal Taylor-neural Lyapunov
Functions". In: arXiv e-prints (2024). DOI: 10.
arxiv. 2408. 17246.
[4] C. Beattie, V. Mehrmann, H. Xu, and H. Zwart. "Linear Port- Port-Hamiltonian Descriptor Systems". In: Mathematics of Control, Signals, and Systems 30.4 (2018), p. 17. DOI: 10. 1007 1007000049-- 018 018-0223-3.
[5] L. Böttcher. "Gradient-Firee Training of Neural ODEs for System Identification and Control Using Ensemble Kalman Inversion". In: ICML Workshop on New Frontiers Learning, Control, and Dynamical Systems. 2023.
[6] S. P. Boyd and L. Vandenberghe. Convex Optimization. Version 29, Cambridge New York Melbourne New Delhi Singapore: Cambridge University Press, 2023. 716 pp.
[7] R. T. Q. Chen, Y. Rubanova, Bettencourt, and D. K. Duvenaud. "Neural Ordinary Differ- ential Equations". In: Advances in Neural Information Processing Systems. Vol, 31. Curran Associates, Inc., 2018.
[8] Y. Chen, Y. Shi, and B. Zhang. "Optimal Control via Neural Networks: A Convex Approach". In: International Conference on Learning Representations. 2019.
[9] Choudhary, F. Lindner, E. G. Holliday, S. T. Miller, S. Sinha, and W. L. Ditto. "Forecasting Hamiltomian Dynamics without Canonical Coordinates". In: Nonlinear Dynamics 103.2 (2021),
pp. 1553-1562. DOI: 10.
[10] M. Cranmer, S. Greydanus, S. Hoyer, P. Battaglia, D. Spergel, and Ho. "Lagrangian Neural Networks". In: ICLR 2020 Workshop on Integration of Deep Neural Models and Differential
Equations. 2019.
[11] S. Desai, M. Mattheakis, D. Sondak, Protopapas, and S. Roberts. "Port-Hamiltniian Neural Networks for Learning Explicit Time-Dependen Dynamical Systems". In: Physical Review E
104.3 (2021), 1) DOL: 10. 1 1101 1003//haaee....... 104. 034312.
[12] E. Dupont, A. Doucet, and Y. W. Tch. "Augmented Neural Odes". In: Advances in Neural Information Processing Systems. Vol. 32. Curran Associates, Inc., 2019.
[13] S. Eidnes, A. J. Stasik, C. Sterud, E. Bøhn, and S. Riemer-Sprennen
Hameltonamm Neural Networks with Sttate-Dependee External Forces". In: Physica D: Nonlinear Phenomena 446 (2023), p. 133673. DOI: 10. 1016 physd. 2023. 133673.
[14] N. B. Erichson, M. Muehlehach, and M. W. Mahoney. IPforres-Informed Autoencoders for Lyapunov-stable Fluid Flow Prediction" In: Machine Learning and the Physical Sciences (NeurIPS 2019). Neural Information Processing Systems. Vancouver, Canada, 2019.
[15] X. Glorot and Y. Bengio. "Understanding the Difficultty of Training Deep Feedforward Neural Networks". In: Proceedings of the Thirteenth International Conference on Artificial Intelli- gence and Statistics. Vol. 9. Sardinia, Italy: PMLR, 2010, pp. 249-256.
[16] S. Greydanus, M. Dzamba, and J. Yosinski. "Hamiltonian Neural Networks". In: Advances in Neural Information Processing Systems. Vol. 32. Curran Associates, Inc., 2019.
[17] E. Hairer, S. P. Nørsett, and G. Wanner. Solving Ordinary Differential Equations I: Nonstiff
Problems. 2nd rev. ed. Springer Series Computational Mathematics Heidelberg London: Springer, 2009. 528 pp.

---

[18] Q. Hernández, A. Badias, D. Gonzalez, F. Chinesta, and E. Cueto. Structure-Preservvii Neural Networks". In: Journal of Computational Physics 426 (2021), p. 109950. DOI: 10. 1016/j jcp. 2020. 109950.
[19] P. Jin, Z. Zhang, I. G. Kevrekidis, and G. E. Karniadakis. "Learning Poisson Systems and Trajectories of Autonomous Systems via Poisson Neural Networks". In: IEEE Transactions on Neural Networks and Learning Systems 34.11 (2023), pp.
DOI: 10. 1199/TNNLS.
2022. 3148734.
[20] P. Jin, Z. Zhang, A. Zhu, Y. Tang, and G. E. Karniadakis. "SympNets: Intrinsic Structure- Preserving Symplectic Networks for Identifying Hamiltonian Systems". In: Neural Networks 132 (2020), pp. 166-179, DOI: 10. 1016/j neunet coaunet.020...... 017.
[21] M. Kannapinn, M. K. Pham, and M. Schäfer. "Physics-Based Digital Twins for Autonomous Thermal Food Processing: Efficient, Non-Intrusive Reduced-Order Modeling". In: Innovative Food Science & Emerging Technologies 81 (2022), 103143. DOI: 10. 1016/ ifset. 2022. 103143.
[22] M. Kannapinn, F. Roth, and O. Weeger. "Digital Twin Inference from Multi-Physical Simula- tion Data of DED Additive Manufacturing Processes with Neural ODEs". 2024. [23] M. Kannapinn, M. Schäfer, and O. Weeger. "TwinLab: A Framework for Data-Elficient Train- ing of Non-Intrusive Reduced-Ordes Models for Digital Twins". In: Engineering Computations
(2024). DOI: 10. 1100/EC23-000000 0855.
[24]
H. K. Khalil. Nonlinear Systems. 3rd ed. Upper Saddle River, N.J: Prentice Hall, 2002. 750 pp.
[25] D. Kingma and J. Ba. "Adam: A Method for Stochastic Optimization". In: International Conference on Learning Representations, 2015.
[26] R. Kojima and Y. Okamoto. "Learning Deep Input-Output Stable Dynamics". In: Advances in Neural Information Processing Systems. Vol. 35, Curran Associates, Inc., 2022, pp. 8187-8198,
[27] J. Z. Kolter and G. Manek. "Learning Stable Deep Dynamics Models", In: Advances in Neural Information Processing Systems. Vol. 32. Curran Associates, Inc., 2019.
[28] Z. Lai, C. Mylonas, S. Nagarajaiah, and E. Chatzi. "Structural Identification with Physics- Informed Neural Ordinary Differential Equations". In: Journal of Sound and Vibration 508
(2021), P. 116196. DOr: 10. j016/j- jsv 116196.
[29] N. Lawrence, P. Loewen, M. Forbes, J. Backstrom, and B. Gopaluni. "Almost Surely Sta- ble Deep Dynamies". In: Advances in Neural Information Processing Systems 33 (2020), pp. 18942-18953.
[30] G. C. Layek. An Introduction to Dynamical Systems and Chuos. 1st ed. 2015. New Delhi: Springer India: Imprint: Springer, 2015. 1 p. DOI: 10. 1000 11-97778--222
[31] K. Lee and F. J. Parish. "Parameterized Neural Ordinary Differential Equations: Applications to Computational Physics Problems". In: Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 477. 7255 (2021), p. 20210162
DOI: 10 1098 /rsaa.
2021.0162.
[32] C. Legaard, T. Schranz, G. Schweiger, J. Drgoña, B. Falay, C. Gomes, A. Iosifidis, M. Abkar, and P. Larsen. "Constructing Neural Network Based Models for Simulating Dynamical Sys- tems". In: Acm Computing Surveys 55.1 11 (2023). DOI: 10. 1145/3567591.
[33] D. Liu and Y. Wang. "Multi-Fidelity Physics-Cond Constrained Neural Network and Its Application in Materials Modeling". In: Journal of Mechanical Design 141.12 (2019), p. 121403. DO1:
10.
[34] M. Lutter, C. Ritter, and J. Peters. "Deep Lagrangian Networks: Using Physics as Model Prior
for Deep Learning". In: International Conference on Learning Representattinss. 2019.
[35] J. G. Malkin. Theorie der Stabilität einer Bewegung. München: R. Oldenbourg, 1959. 400 pp-
[36]
S. Massaroli, M. Poli, M. Bin, J. Park, A. Yamashita, and H. Asama. "Stable Neural Flows".
In: arXiv e-prints (2020), S.Kiv:2003.080662 DOI: 10. 48544/8550/00 /ar/vv/rr/vvv... 2003 08063.
[37] H. Nakano, R. Ariizumi, T. Asai, and S.-i. Azuma. "Model Estimation Ensuring Passivity by Using Port- Hor--Hamltonian Model and Deep Learning", In: 61st Annual Conference of the Society of Instrument and Control Engineers (SICE). Kumamoto, Japan: IEEE, 2022, pp. 886-891. 2022. 9905555..
DO1: 23919/SICESE
[38] C. Neary and U. Topcu. "Compositional Learning of Dynamical System Models Using Port- Hamiltonian Neural Networks". In: Proceedings of the 5th Annual Learning for Dynamics and Control Conference. Vol. 211. PMLR, 2023, pp. 679-691.
12

---

[39] Y. Okamoto and R. Kojima. "Learning Deep Dissipative Dynamics". . In: Proceedings of the AAAI Conference on Artificial Intelligence. Vol. 39. 2025, pp. 19749-19757. DOI: 10. 1609/ aaai. v39118....... 34175.
[40] J. S. R. Park, S. W. Cheung, Y. Choi, and Y. Shin. "tLaSDI: Thermodynamiis-Infomoe Latent Space Dynamics Identification". In: Computer Methods in Applied Mechanics and Engineering 429 (2024), p. 117144. DOI: 10. 1016/j j. cma. 2024. 117144.
[41] Pascanu, T. Mikolov, and Y. Bengio. "On the Difficulty of Training Recurrent Neural Networks". In: Proceedings of the 30th International Conference on Machine Learning. Vol. 28. Proceedings of Machine Learning Research 3. Atlanta, Georgia, USA: PMLR, 2013, PP. 1310-1318.
[42] S. Ramasinghe, H. Saratchandran, V. Shevchenko, and S. Lucey. "On the Effectiveness of Neural Priors in Modeling Dynamical Systems". In: JCML Workshop on New Frontiers in Learning. Control, and Dynamical Systems. 2023.
[43] J. Rettberg, Kneifl, J. Herb, P. Buchfink, J. Fehr, and B. Haasdonk. "Data-Driven Identifica- tion of Latent Hott-Hamiltonian Systems". In: arXiv e-prints (2024). DOI: 10. 48550//55/// arxiv.
2408. 08185.
[44] S. M. Richards, Berkenkamp, and A. Krause. "The Lyapunov Neural Network: Adaptive Stability Certification for Safe Learning of Dynamical Systems". In: Proceedings of the 2nd Conference on Robot Learning. Vol. 87. PMLR, 2018, pp. 466-476.
[45] R. T. Rockafellar. Corvex Analysis. Princeton University Press, 1970.
[46] A. Schlaginhaufen, P. Wenk, A. Krause, and F. Dorfler. "Learning Stable Deep Dynamics Models for Partially Observed or Delayed Dynamical Systems". In: Advances in Neural Information Processing Systems. Vol. 34. Curran Associates, Inc., 2021, pp. 11870-11882.
[47] M. Schoukens, P. Mattson, T. Wigren, and J.-P. Noël. Cascaded Tanks Benchmark Com- bining Soft and Hard Nonlinearities. .mat .csv. Version 1. (License: CC BY-SA 4.0). 2020. DOI: 10.
[48] M. Schoukens, Mattsson, T. Wigren, and J.-P. Noel. "Cascaded Tanks Benchmark Combining Soft and Hard Nonlincarities". In: Workshop on Nonlinear System Identification Benchmarks
2016. Brussels, Belgium, 2016.
[49] R. Seydel. Practical Bifurcation and Stability Analysis. 3rd ed. Interdisciplinary Applied Mathematics 5. New York: Springer, 2010. 483 pp.
[50] N. Takeishi and Y. Kawahara. "Learning Dynamics Models with Stable Invariant Sets". In: Proceedings of the AAAI Conference on Artificial Intelligence 35.11 (2021), pp. 9782-9790.
DOI: 10. 1 /aaaiaai.... www.w11....11............1...111
[51] J. van der Schaft and D. Jeltsema. Port-Hamiltonian Systems Theory: An Introductory Overview. Foundations and Trends in Systems and Control 1, 2/3 (2014). Boston Delft: Now,
2014. 216 pp.
[52] F. Verhulst. Nonlinear Differential Equations and Dynamical Systems. Universitext. Berlin New York Paris [etc.]: Springer, 1990.
[53] R. Wang and R. Yu. "Physics-Guided Deep Learning for Dynamical Systems: A Survey". In:
arXiv e-prints (2 DOI: 10. 48550/ arXiv .2..7.01277 01272.
[54] J. Willard, X. Jia, S. Xu, M. Steinbach, and V. Kumar. "Integrating Scientific Knowledge with Machine Learning for Engincering and Environmental Systems". In: Acm Computing Surveys
55.4 (2022). DOI: 10.
[55] H. Yu, X. Tian, W. E, and Q. Li. "OnsagerNet: Learning Stable and Interpretable Dynamics Using a Generalized Onsager Principle". In: Physical Review Fluids 6.11 (2021), p. 114402.
DO1: 10. 11003/PhyRevFluiie
114402.
[56] Z. Zhang, Y. Shin, and G. Em Karniadakis. "GFINNs: GENERIC Formalism Informed Neural Networks for Deterministic and Stochastic Dynamical Systems". In: Philosophical Transac- tions of the Royal Society A: Mathematical, Physical and Engineering Sciences 380.2229
(2022), P. 20210207 DOI: 1............111 0098/rsta.200 2021.0207.
[57] Y. D. Zhong, B. Dey, and A. Chakraborty. "Dissipative SymODEN: Encoding Hamiltonian Dy- namics with Dissipation and Control into Deep Learning", In: ICLR Workshop on Integration of Deep Neural Models and Differential Equations. 2020.
DOI: 10. 48550/ arxiv. 2002. 08860.
13

---

A Proof of stability
This section provides proof of Theorem 3.1, which describes the essential requirements for achieving stable dynamics with port-Hamillonian systems. Due to the relevance of convexity to the following discussion, the following Lemma first provides some widely known equivalent formulations of convexity for scalar- scalar-valued valued multivariate functions.
Lemma A.1. Consider a twice continuously d f 12 - R, x + f(x). If convex, the following statements hold:
1. (x + (1 - (x) 2 f()x1 + (1 2))22) V0 22 E R", A [0,
2. 0/1
(x2 x1)= - f(x2) - f(x1) V01, E R,,
3.
8' Dada
( = 0
Ve E IR.
For proof, the reader is referred [6]. Note that for strict convexity, the inequalities in statements 1 and 2 must be strict for all æ * Ø2, and the Hessian in statement 3 must be positive definite. With the notion of convexity established, attention can now be directed the proof of Theorem 3.1. For convenience, the statement of the theorem is repeated verbatim below:
Theorem A.2. Consider the PHS Equation (3) in the unforced case u = 0:
0H, (æ)) with J= R R 2.
a -= (xx) R(x))
(A.1)
dæ
Suppose the Hamiltonian H(xx) is convex, twice continuously differentiable, and fulfills:
OH = =0
8² > O.
H(0) =
(A.2)
dx ==
dxθx ==
Then, the system in Equation (A.1) has a stable equilibrium at x=()) = 0, and all solutions are bounded. Furthermore, the equilibrium is globally asymptotically stable if R(x) > 0.
Proof. Since the gradient of the Hamiltonian vanishes and its Hessian is positive definite at æ = 0, has strict local minimum at the origin. This implies the existence of an R > 0 such that
H(x)) H()) = 0 Vo E {||=||| S R, x + 0}.
(A.3)
Consider the function g(x) ||< for some E (0, R with a = inf ||m -r H(x) > 0. We claim
that g is a lower bound for H. i.e. g(x) < H(xx) for all ||=|| >
Suppose there existed ø E Rn that violated the claim, that is | > r and g(Ø) H()). Due to the convexity of H, it holds
XH(#) + ()H(O) > H(A (1 A) 0) AA (0, 1]
(A.4)
r and thus Ag(Ø)= > XH( 2 H(x)). Choosing À = /
and applying the definition of g we obtain
g( H
||ق|||
za
a
(A.5)
This is contradiction, and therefore is lower bound.
With the bound established, showing that H is radially unbounded is straightforward For any path (t) with 4o ||x(t]] = oc, once [[=( >> r we have
lim H(x(t)) 2 lim g(x(()) D0. t->oc
(A.6)
t-000
We would like to remark that the same conclusion could alternatively be reached by showing that H has at least one bounded sublevel set due its strict minimum and then concluding radial unboundedness (i.e., coercivity) using Corollary 8.7.1 in [45] together with the standard equivalence for proper closed convex functions. However, the presented approach using a lower bound can directly be extended to show global positive definiteness of H. For this, note that for all ||>||| > the lower bound implies < ((x) S H(x)). Together with Equation (A.3), this implies that H is indeed positive definite for all R"
14

---

From the not-HHmiltonian structure of Equation (A.1) we immediately obtain the decrease condition H. & O (see Equation (4)). Using the H. as Lyapunov function, we obtain all requirements to conclude local stability of the equilibrium. Furthermore, the decrease condition implies that solution starting at x" cannot leave the sublevel set Ux {x [?(() I H(x")). Since OH S S2g fæ IR" (g(2)) S H(x*)} and f2 is clearly bounded, all solutions of Equation (A.1)
must be bounded well.
Finally, show that global asymptotic stability follows from R > 0, it remains be demonstrated that H < 0 for all æ: + 0. This is implied by Equation (4) if 0H for all x: o O. Proving that the gradient of H. vanishes only at the origin is done by contradiction. Assume that there exists an ~ + 0 Oz O. Furthermore, let x= = 0. Since H. is convex it follows (see Lemma A. that
the
such that 0
0H
H(x*) - H()) 2
(x*
(A.7)
dr
07|
H( - H(x")
(Ø x" 0.
(A.8)
dx
This implies H(x') H() 0. However, since H is globally positive definite and thus only vanishes at the origin, we have x² = x. This is contradiction as ø + O, which concludes the proof.
B Experiment hyperparameters
In each experiment presented within this manuscript, the different models were trained for the same number of steps using identical learning rates. The mean squared error loss function and the ADAM [25] optimizer were used throughout. The model and training hyperparameters, along with the number of model instances trained for statistical evaluation, are provided in Table Throughout the experiments presented in Sections 4.1 4.3, all NNs consisted of two hidden layers, each containing 16 neurons. The network width was increased to 32 neurons for the additive manufacturing surrogate in Section 4.4. All network weights were initialized using the Glorot uniform [15] scheme.
The PHS beSe- models (sPHNN, sPHNN-LM, PHNN, bPHNN) share the same structure, differing only in how the Hamiltonian H. is parametrized. Specifically, the sPHNNs use an FICNN with normalization term, as described in Section The only difference for the sPHNN -LM variant is that the minimum xx of the Hamiltonian is learned use an unconstrained FFNN without any normalization. Finally, the bPHNNS ensures positivity and radially unboundedness of the Hamiltonian by adopting the method described by Yu al. [55] for the OnsagerNets Energy function. Here, the output of an unconstrained network is squared and combined with quadratic regularization term. However, unlike OnsagerNet, we used standard FFNNs instead of ResNets.
additional parameter. The PHNNs, by contrast,
The models were trained Windows machine with Intel 17-13700K (13th Generation) and 32 GB of memory. Table 3 provides representative computational costs for training and prediction. Here, Derivative evaluation refers the average time taken by each model compute the state derivative t
Table 2. Hyperparamcters per experiment. TF and DF refer to trajectory and derivative fitting, respectively. The number of instances trained per model type listed in the "# inst." column, and and denote the number dimensions of state and input. The final three columns list parametrization SPHNN-LM, bPHNN, and PHNN as either state-dependent matrices (s.) or constant matrices (c.).
matrices sPHNN.
Experiment
training steps
learning rate
inst.
n
m
matrices R
G
Randomly initialized model
1
2
8.
8.
Spinning rigid body Cascaded tanks
00000 DF
1 X 10
10
3
50 000 TF
5 10
20
2
1
C.
C.
C.
Food processing surrogate
30 000 TF
1 10
20
2 to
5
1
Additive manufacturing surrogate
20 000 DF
1 10
20
40
40
and 10 000 TF
1 10
15

---

Table 3. Representative computational costs for training and prediction. The top report the training times for derivative and trajectory fitting, respectively. bottom two rows show the time single derivative evaluation and full trajectory integration using data and models from Section 4.3.
sPHNN
SPHNN-LM
SNODE
bPHNN
PHNN
NODE
Derivative fit Section 4.1 (min)
0.577
0.644
0.672
0.635
0.539
0.346
Trajectory fit Scction 4.3 (min)
37.663
39.160
38.329
19.111
17.951
Derivative evaluation (mus)
0.007
0.007
0.007
0.007
0.005
Model integration (ıus)
0.365
0.393
0.422
0.237
0.221
for given state and input u, which needs be evaluated for time integration. In contrast, Model integration denotes the time required to generate a single trajectory prediction via integration. cases, the reported integration time also depends on the
While the latter better reflects practical integration scheme and chosen step size. Since we use an adaptive step size controller, the step size turn depends on the learned dynamies.
C Spinning rigid body system
The governing equations of the diree-dimeniiona spinning rigid body used in the experiment in Section 4.1 are described in Section 2.2. We generated (w, (w,w) pairs for training by evaluating
Equation (5) for varying w. To simulate a measurement process, we selected w values from 50s long
numerically integrated trajectories, using 10 randomly chosen initial conditions (0) E [0, 1l with
parameters diag(1. [1, 2, 3) and μ = 0.01. A representative trajectory of w is shown in Figure C.I.
NODE
PHNN
SPHNN
4PHNN-LM
NODE
0.5
BRHNNN
0
-0.5
50 Time t/s
100
0
50
100
Time t/s
Figure C.I. Representative, extended training trajectory for spinning rigid body. including interquartile mean and range of model predictions for wa. Left: Models without stability bias; Right: Models with stability bias.
D Cascaded tanks data
The cascaded tanks dataset [47, 48] used in Section 4.2 consists of measurements of fluid level control system. Liquid flows from an upper tank through an outlet into a lower tank, exits via similar outlet into a reservoir, and is pumped back to the upper tank, see Figure D. The system's 's input is the voltage applied to the pump, with higher values resulting in a faster volume flow. When the pump speed is large enough, an overflow fluid level the lower tank, and the resulting voltage represents the system's output y. The dataset consists of two trajectories comprising 1024 pairs (u, y) each and Figure D. depicts the training trajectory. The initial water level in the upper tank is unknown, but the same for both trajectories.
occur in the upper tank. A capacitive sensor measures the
Besides the experiments presented in Section 4.2, we further investigated the models' capabilities to learn the hard saturation that the cascaded tanks data system exhibits at an output value of V This occurs when the tanks overflow in response to large input voltage over an extended period of time (see Figure Dbl [48]. We explore the models' behavior in this regime by constructing pump voltage signal that linearly increases from 1 V initially to 8 V at s800s. The resulting predictions
16

---

10
Water level sensor
Pump
5
1,000
2,000
3,000
4,000
t/s
(a) Experimental setup
(b) Cascaded tanks training data (overflow highlighted).
NODE
SPUNN 4PHNN.LM NPHNN
20
33/V
10
1,000
2,000
1,000
2,000
Time
Time t/s
(c) Saturation.
Figure D. Cascaded tanks: (a) System setup, (h) Training trajectory. Predictions for linearly increasing pump voltage. Due to an overflow, true system exhibits hard saturation 222 = 10 V. Lines correspond to the interquartile mean, whereas shaded areas represent the interquartile range of the predictions from the 20 model instances.
shown in Figure D. For pump speeds up ~ V (~2400 the sPHNN, SPHNN-LM, and bPHNN instances capture the saturation effect, while the NODE and PHNN instances fail to do so.
E
Thermal food processing surrogate data
This section provides additional details about the data used Section 4.3, which is subset of the data generated and used by Kannapinn et al. [21] in the creation of ROMs. The dataset describes the thermal processing chicken breasts in a convection oven and was generated
the simulation of soft-matte model for meats implemented in COMSOL Multiphysics. Each trajectory in the dataset describes the temperature history at two points inside the full order model. These result from predefined excitation signal that controls the oven temperature. The probe points are situated in the center and close to the surface of the food item, with the corresponding temperatures labeled as Ta and TB, respectively. To facilitate cross-referencingg the alphanumerie identifiers used by Kannapinn et al. [21] for the individual trajectories are retained. The training data consists of trajectories 745 and 795. They use an amplitude-modulated pseudo-random binary sequence (APRBS) and multi-sine as the excitation signal Towen- Figure E. (left and middle) illustrates Toven along with the and surface temperatures for both trajectories. This selection provides fair evaluation of the models across the relevant temperature spectrum. For studying the impact of the number of training trajectories on the models, we additionally included 13 trajectories with identifiers 835, 851, 853, 861, 868, 873, 874, 888, 899, 902, 905, 906 and 917. In the experiment, the test group AP15 [21] utilized as test data. It comprises 15 trajectories with APRBS excitations, specifically selected to ensure an even distribution of amplitudes and median values. Figure E. (right) shows an exemplary test trajectory, including the predictions for Ta from all models using augmented dimensions. All trajectories span 1395 s and consist of 280 samples. The data is normalized for training to ensure the numerical values are of magnitude one. For this, an affine transformation is applied to Towen TA, and Tn individually. It shifts the ambient temperature to a value of zero and rescales the signal to have unit variance. The scaling and shifting factors are calculated using only the training data to ensure that the test data does not influence the training procedure.
17

---

Trajectory 745
Trajectory 795
Trajectory 462
450
Ta Tyy
Teees
Tp
sPHNN NPHNN
APHNN-LM PIINN
400
NODE
350
300
G00
1,000
500
1,000
500
1,000
Time t/s
Time f/s
Time f//s
Figure E.I. Left and middle: Training trajectories. Right: Exemplary test trajectory with interquartile mean and range of the model predictions for TA-
F Thermal field data
F1 Dimensionality reduction
The following describes the POD method for dimensionality reduction of the field data applied in Section 4.4. The process is first described for general field data given as a matrix A = [ao..... aM J of M snapshots a; E RRY, where M is the number of snapshots and N is the number of finite element nodes. To identify the dominant modes of the data, we perform a singular value decomposition (SVD) of A:
= UEVT E RRNN
(F.I)
where U E RRMM and V E N are unitary matrices and E E R MxN is a rectangular diagonal matrix containing the singular values 01 2 ... 2 оM on its diagonal. Dimensionality reduction is achieved by truncating the SVD A UEVT, where U E Mxx and V e R contain the first n columns of U and V. and Σ diag ( (TI...... . .,, ) is the truncated matrix of the dominant singular values. The dominant modes in the data are then given by M <V. Here, the scaling factor e the standard deviation of the elements in US. This scaling ensures that the latent trajectories will have unit variance. To map a latent state E IR to the full-order field X E R.V we compute:
X=Ma
(F2)
Conversely, to map from the full-order field to the latent space, we solve the least-squares regression problem:
x arg min || Ma'I MTX.
(F3)
In Section 4.4, we consider two fields: The thermal field describing the state of the system and the excitation given by the heat source field. For each, a separate SVD is performed to compute My and M, using the respective snapshots from the two training trajectories. We thus obtain two pairs of encoders & and decoders 6z
(T) = arg min IT - MTæel,
Ør ΦT(#) MTæ,
(F.4)
b (q) arg min l|q Maul,
0g(u) Mau. u.
(F.5)
To fully leverage the benefits of sPHNNs, it is essential to ensure that the absence of a heat source encoded the zero-vector in the corresponding latent space. Since the mapping from the latent space Equation (F.2) is linear, and q O corresponds to the absence of a heat source in the full-order space, this requirement is automatically fulfilled. However, for applications where a nonzero state describes the absence of an external energy source or for general nonlinear latent mappings (such as autoencoders), the condition can be met by applying a constant translation in the latent space, which shifts the encoded state ã +(()) to the origin. In Section 4.4, we use this method to ensure that the thermal equilibrium corresponds to the origin of the temperature latent space. While this not strictly necessary, knowledge of the cquilibrium location in the latent space enables its integration
18

---

into the sPHNN's architecture by fixing the stable equilibrium position. Consequently, we arrive at the representation of the latent state
(T) 1 (Tey);
T = (r t/T (Tey))
(F.6)
where Tag is the full-order state vector describing thermal equilibrium. The complete model archi- tecture represented Figure F. There, the function f is used as a placeholder for an arbitrary evolution equation, to generalize to all model architectures used in Section 4.4.
Field input
Latent space
Field output
T(u)
f(x,
J
T(1)
initial conditon
((1)
4ig
Figure F.I. Model setup used Section 4.4. The initial temperature field (fo) and input field g(t) are mapped via encoders Φr, Øy placeholder for a specific model, e.g., NODE or sPHNN. Finally, the predicted latent trajectory (t) is mapped to the full-order field via the decoder dr.
latent space, where the latent dynamics f are integrated. Here, used as
Finally, we note that alternative dimensionality reduction techniques could be applied. In preliminary experiments, we explored the of autoencoders for this purpose. Although autoencoders offer advantages over the linear SVD regarding achievable compression rates, they exhibited poor extrapo- lation performance when trained on sparse data. This issue could potentially be mitigated by training the autoencoders with additional synthetic data, which would not necessitate costly simulations; however, further research is needed to explore this approach.
F.2 Additional results
We briefly present additional results from Section 4.4 to demonstrate that the differences in predictions across the various model architectures remain consistent across multiple instances. Figure F.2 depicts the RMSEs of each model from Section 4.4 computed over all 25 trajectories and the spatial domain. The behavior of the models with stability bias can be distinguished from the NODE and PHNN, which suffer from instabilities that result in rapidly increasing errors. While both sPHNN variants and bPHNNs exhibit decreasing errors the system returns thermal equilibrium, sPHNNs achieve a peak RMSE nearly an order of magnitude lower, due to their strong inductive bias from enforcing global stability. The sPHNN-LM performs almost identically to the sPHNN initially, while showing slightly larger errors in the convergence to equilibrium. This is due to slight errors in the learned equilibrium position, which can be addressed by increasing the amount of training data, as demonstrated in Section 4.3.
10
APINN
RMSE/K
SPHNN LM PHNN
1
10a
10a
Time t/s
Figure F.2. Interquartile mean and range of the RMSEs evaluated with all 25 trajectories of the thermal field data.
19

---

NeurIPS Paper Checklist
1. Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
Justification: Our claim, that the enforcement of global stability and physics-based condi- tions results in a strong inductive bias and increased performance in applicable cases, is highlighted in the abstract and introduction, and demonstrated by our empirical results.
Guidelines:
The answer NA means that the abstract and introduction do not include the claims made in the paper.
The abstract and/or introduction should clearly state the claims made, including the contributions made the paper and important assumptions and limitations. A No or NA answer to this question will not be perceived well by the reviewers.
The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings.
is fine to include aspirational goals as motivation as long as is clear that these goals not attained by the paper.
2. Limitations
Question:
Does the paper discuss the limitations of the work performed by the authors?
Answer: [Yes]
Justification: Although do include dedicated "Limitations" section, the paper explicitly discusses the assumptions and limitations of our approach, such as the global stability requirement and its implications for model expressiveness, relevant points throughout the text.
Guidelines:
The answer NA means that the paper has limitation while the answer No means that the paper has limitations, but those are not discussed in the paper.
The authors are encouraged to create separate "Limitations" section in their paper.
The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g.. independence assumptions, noiseless settings, model well-specification, asymptotic approximations only holding locally). The authors should reflect on how these assumptions might be violated in practice and what the implications would be.
The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on a few datasets or with few runs. In general, empirical results often depend on implicit assumptions, which should be articulated.
The authors should reflect the factors that influence the performance of the approach. For example, facial recognition algorithm may perform poorly when image resolution low or images are taken in low lighting. Or speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails to handle technical jargon.
The authors should discuss the computational efficiency the proposed algorithms and how they scale with dataset size.
If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and fairness.
While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, worse outcome might be that reviewers discover limitations that aren't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an impor- tant role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed not penalize honesty concerning limitations.
3. Theory assumptions and proofs
20

---

Question: For each theoretical result, does the paper provide the full set of assumptions and complete (and correct) proof?
Answer: [Yes]
Justification: Our theoretical result provided in Theorem 3.1 are accompanied by intuitive discussion Section and rigorous proof the appendix Appendix A.
Guidelines:
The answer NA means that the paper does not include theoretical results. All the theorems, formulas, and proofs in the paper should be numbered and cross- referenced.
All assumptions should be clearly stated or referenced in the statement of any theorems.
The proofs can either appear in the main paper or the supplemental material, but if they appear in the supplemental material, the authors are encouraged provide short proof sketch to provide intuition.
Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix or supplemental material.
Theorems and Lemmas that the proof relies upon should be properly referenced.
4. Experimental result reproducibility
Does the paper fully disclose all the information needed to reproduce the main perimental results of the paper to the extent that it affects the main claims and/or conclusions of the paper (regardless whether the code and data are provided or not)?
Question:
Answer: [Yes]
Justification: We provide extensive description of the proposed model architecture as well as the experimental setup. This includes the model hyperparameters, which are listed in Appendix B.
Guidelines:
The answer NA means that the paper does not include experiments.
If the paper includes experiments, a No answer to this question will not perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not.
If the contribution is dataset and/or model, the authors should describe the steps taken make their results reproducible verifiable.
Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is a novel architecture, describing the architecture fully might suffice, or if the contribution be necessary to either make it possible for others to replicate the model with the same dataset, or provide access to model. In general. releasing code and data often one good way to accomplish this, but reproducibilityy can also be provided via detailed instructions for how to replicate the results, access to a hosted model (e.g.. in the case of large language model), releasing of a model checkpoint, other means that are appropriate to the research performed.
specific model and empirical evaluation, may
While NeurlPS does not require releasing code, the conference does require all submis- sions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) If the contribution is primarily a new algorithm, the paper should make it clear how reproduce that algorithm. (b) If the contribution is primarily the architecture clearly and fully.
model architecture, the paper should describe
(c) the contribution is a new model (c.g., a large language model), then there should either be way to access this model for reproducing the results or way to reproduce the model (e.g., with an open- open-source dataset or instructions for how to construct the dataset).
(d) We recognize that reproducibility may be tricky in some cases, in which case authors are welcome describe the partieular way they provide for reproducibiliity In the case of closed-source models, it may be that access to the model is limited in
21

---

some way (e.g., to registered users), but it should be possible for other researchers have some path reproducing or verifying the results.
5. Open access to data and code
Question: Does the paper provide open access to the data and code, with sufficient instruc- tions faithfully reproduce the main experimental results, as described in supplemental material?
Answer: [Yes]
Justification: Code and model weights are provided in this GitHub repository: https: //github. com/ CPSCcb CCPShubssphnn httpo//sphnnnuuiicattteeeee
The repository includes instructions on how to set up the environment, access the data, and run experiments to reproduce results of the paper.
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurIPS code and data submission guidelines (https://nips cc/
htt.guudddee/oeeeubbbieeeee sonPolicy) for more details.
While we encourage the release of code and data, we understand that this might be possible, so "No" is an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (e.g., for new open-source benchmark).
The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurlPS code and data submission guidelines (https: // httppss.cc/pullcc///e/e// honoo///cgiddess/deeuubbiiiiiiii CodeeSubissiiii stt.uubmiiissoooooooccoooo for more details.
The authors should provide instructions data access and preparation, including how to access the raw data, preprocessed data, intermediate data, and generated data, etc.
The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. only subset of experiments are reproducible, they should state which ones are omitted from the script and why.
At submission time, to preserve anonymity, the authors should release anonymized versions (if applicable).
Providing as much information as possible in supplemental material (appended to the paper) recommended, but including URLs to data and code permirted.
6. Experimental setting/details
Question: Does the paper specify all the training and test details (e.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary understand the results?
Answer: [Yes]
Justification: We supply the essential training and test details necessary understand the results in the main paper, and provide additional information about the data and model hyperparameters the appendix and code repository.
Guidelines:
The answer NA means that the paper does not include experiments.
The experimental setting should be presented in the core of the paper to level of detail that is necessary to appreciate the results and make sense of them.
The full details can be provided either with the code, in appendix, or as supplemental material.
7. Experiment statistical significance
Question: Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Answer: [Yes]
Justification: In the discussion of our experiments, we provide suitable measures on the variance between individual initializations of the different models.
Guidelines:
22

---

The answer NA means that the paper does not include experiments.
The authors should answer "Yes" if the results are accompanied by enror bars, confi-
dence intervals, or statistical significance tests, at least for the experiments that support the main claims of the paper.
The factors of variability that the error bars capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall run with given experimental conditions).
The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.)
The assumptions made should be given (e.g., Normally distributed errors).
should be clear whether the error bar is the standard deviation the standard error of the mean.
It is OK to report I-sigma error bars, but should state it. The authors should preferably report a 2-sigma error bar than state that they have 96% CI, if the hypothesis of Normality of errors is not verified.
For asymmetric distributions, the authors should be careful not show in tables figures symmetric error bars that would yield results that are out of range (c.g. negative error rates).
error bars are reported in tables or plots, The authors should explain the text how they were calculated and reference the corresponding figures or tables in the text.
8. Experiments compute resources
For each experiment, does the paper provide sufficient information on the com- puter resources (type of compute workers, memory, time of execution) needed to reproduce the experiments?
Question:
Answer: [Yes]
Justification: Given that the target application involves small neural networks trained on sparse data, the computational requirements are relatively modest. We have specified the system used generate the reported results, along with approximate training times, in Appendix B.
Guidelines:
The answer NA means that the paper does not include experiments.
The paper should indicate the type of compute workers CPU or GPU, internal cluster, cloud provider, including relevant memory and storage.
The paper should provide the amount of compute required for each of the individual experimental rums as well as estimate the total compute.
The paper should disclose whether the full research project required more compute than the experiments reported the paper (e.g., preliminary or failed experiments that didn't i make it into the paper).
9. Code of ethics
Question: Does the research conducted the paper conform, in every respect, with the NeurlPS Code of Ethics https: http://eurrrrr.cccoo httprr/ss...cppbpl/e//eeeeeteedd
Answer: [Yes]
Justification:
The present research conforms with the NeurIPS Code of Ethics.
Guidelines:
The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.
If the authors answer No, they should explain the special circumstances that require a deviation from the Code of Ethics.
The authors should make sure to preserve anonymity (c.g.. if there special consid- eration due to laws or regulations in their jurisdiction).
10. Broader impacts
Question: Does the paper discuss both potential positive societal impacts and negative socictal impacts of the work performed?
23

---

Answer: [NA]
Justification: This work focuses on foundational research involving small-scale, physically structured neural networks trained on sparse data for modeling dynamical systems. It is not tied to specific applications or deployments, and we do not foresee any direct societal impacts arising from its current form.
Guidelines:
The answer NA means that there is no societal impact of the work performed.
If the authors answer NA or No. they should explain why their work has no societal impact or why the paper does not address societal impact.
Examples of negative societal impacts include potential malicious unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g.. deploymentt of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations.
The conference expects that many papers will be foundational research and not tied to partieular applications, let alone deployments. However, if there is direct path to any negative applications, the authors should point it out. For example, it is legitimate point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster.
The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology.
If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition attacks, mechanisms for monitoring misuse, mechanisms to monitor how system learns from feedback over time, improving the efficiency and accessibility of MI.).
11. Safeguards
Does the paper deseribe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Question:
Answer: [NA]
Justification: This paper does not involve the release of models or datasets that pose a high risk of misuse.
Guidelines:
The answer NA means that the paper poses such risks.
Released models that have high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters.
Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images.
We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make best faith effort.
12. Licenses for existing assets
Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?
Question:
Answer: [Yes]
Justification: All external assets used in this work are properly cited in the bibliography, including relevant version information, URLs, and licenses where applicable. We ensured
24

---

that the license terms and conditions of use were respected and named the license explicitly in the bibliography. Guidelines:
The answer NA means that the paper does not use existing assets.
The authors should cite the original paper that produced the code package or dataset.
The authors should state which version of the asset used and, possible, include a URL.
The name of the license (e.g., CC-BY 4.0) should be included for each asset.
For scraped data from a particular source (c.g.. website). the copyright and terms of service of that source should be provided.
assets are released, the license, copyright information, and terms use in the package should be provided. For popular datasets, paperswi thcode com com/datasets has curated licenses for some datasets. Their licensing guide can help determine the license of dataset.
For existing datasets that re-packaged, both the original license and the license of the derived asset (if has changed) should be provided.
If this information is not available online, the authors are encouraged to reach out to the asset's ss creators.
13. New assets
Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer: [NA]
Justification: This work does not release any new datasets, models, other assets.
Guidelines:
The answer NA means that the paper does release new assets.
Researchers should communicate the details the dataset/codeemoodl as part their submissions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent was obtained from people whose asset is used.
At submission time, remember to anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file.
14. Crowdsourcing and research with human subjects
Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well as details about compensation (if any)?
Answer: [NA]
Justification:
This work does not involve crowdsourcing or research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing research with human subjects.
Including this information the supplemental material is fine, but if the main contribu- tion of the paper involves human subjects, then much detail as possible should be included the main paper.
According to the NeurlPS Code of Ethics, workers involved in data collection, curation, other labor should be paid at least the minimum wage in the country of the data collector.
15. Institutional review board (IRB) approvals or equivalent for research with human subjects
Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or an equivalent approval/review based on the requirements of your country or institution) were obtained?
25

---

Answer: [NA]
Justification: This work does not involve human subjects, user studies, or any form of data collection from individuals.
Guidelines:
The answer NA means that the paper does not involve crowdsoureing nor research with human subjects.
Depending on the country in which research is conducted, IRB approval (or equivalent) may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.
We recognize that the procedures for this may vary significantly between institutions and locations, and expect authors adhere the NeurIPS Code of Ethics and the guidelines for their institution.
For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review.
16. Declaration of LLM usage
Question: Does the paper describe the usage of LLMs if it is important, original, or non-standard component of the core methods this research? Note that the LLM is used only for writing, editing, formatting purposes and does not impact the core methodology, scientific rigorousness, or originality of the research, declaration is not required.
Answer: [NA]
Justification: This work does not involve large language models in any part of the core research methodology.
Guidelines:
The answer NA means that the core method development in this research does not involve LLMs as any important, original, non-standard components.
Please refer to L.LM policy (https: http:///eurrspeeeee http://urripsscccoo/eeeeeeeeeeee ecc/oooffeencces/22255 for what should or should not be described.
26


Reviews:
[{"review_id": "Kz0wcwbFK3", "reviewer": "Reviewer_FvJb", "summary": "This work concerns an introduction of new neural network architecture for data-driven dynamical system identification based on port-Hamiltonian systems. The architecture is based on the port-Hamiltonian system formulation, in which (1) the structure matrix and the dissipation matrix are each modeled by feed forward neural networks and their structure is enforced (skew-symmetry and symmetric positive semi-definite), (2) the port-Hamiltonian is modeled by a fully input convex neural network (FICNN). The authors show that this new formulation is asymptotically stable, based on a theorem developed in the work. The new model is tested on a few dissipative systems, and shows improved asymptotic behavior as well as superior overall accuracy.", "questions": "- Can this framework be extended to parametrized systems?\n\n- Why is the Lyapunov function introduced, if it is not required in the Theorem? It appears the function is only necessary in the proof?\n\n- Is there another baseline comparison beyond deep learning models (that also incorporate classical methods like SINDy: Lee, Kookjin, Nathaniel Trask, and Panos Stinis. \"Structure-preserving sparse identification of nonlinear dynamics for data-driven modeling.\" Mathematical and Scientific Machine Learning. PMLR, 2022.)?", "limitations": "- The explicit specification using $\\mathbf{x}_0$ should be elaborated upon. More examples where equilibrium points are not known a priori would be enlightening.", "rating": 5, "confidence": 5, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 4, "significance": 2, "originality": 2, "strengths_and_weaknesses": "Strengths\n\n- The presentation is clear and straightforward. The new architecture as well as the main theorem is explained clearly.\n\n- The specific formulation of the port-Hamiltonian is new to this reviewer (although the use of input convex neural networks for Hamiltonian systems is not new).\n\nWeaknesses\n\n- The developed theorem requires a specification of a special point in the state space $\\mathbf{x}_0$. This requirement is a bit peculiar, given that generally stability conditions are often formulated without specifying special points. This point is not sufficiently elaborated upon in the paper; how do you select this point in general situations? Does one need to retrain the model if the point needs to be changed?", "was_revised": true, "final_justification": "The reviewers have provided corrections and have addressed my confusion regarding the knowledge of the equilibrium point $x_0$ and provided evidence that the methodology can still apply when $x_0$ is unknown. The authors agreed to incorporate clarifying comments regarding the key scope of the paper.", "rebuttal": "We thank the reviewer for the insightful questions and the opportunity to clarify important aspects of our approach. Below, we address each point, providing explanations and outlining planned revisions to improve clarity and completeness in the manuscript.\n\n\n> The developed theorem requires a specification of a special point in the state space $x_0$. [...] This point is not sufficiently elaborated upon in the paper; how do you select this point in general situations? \n\nThe proposed model incorporates the equilibrium position $x_0$​ directly into its architecture, but its value does not need to be known in advance.\nIn the theoretical parts of the paper, such as Theorem 3.1, we assume $x_0=0$ to simplify notation. This is without loss of generality and standard practice in stability theory literature [Verhulst 1990, Khalil 2015]. However, in a general application where $x_0$ is unknown, we treat it as a trainable parameter that is optimised during training. This way, arbitrary equilibrium positions can be inferred from observations of the system. We demonstrate this in section 4.2 for the cascaded tanks example, where $x_0$ is learned from the measurement data. \nWhile learning $x_0$ is possible, we argue that for many relevant and globally stable systems the equilibrium can be inferred a priori from physical reasoning, as demonstrated in our other numerical examples. Including this knowledge as an inductive bias into the model can significantly improve model quality. \n\nTo emphasise the above discussion more explicitly, we will revise line 197 to read: “Fixing $x_0$ to the true equilibrium position introduces an additional inductive bias. While this information can often be obtained from prior physical knowledge, it may not always be available. In such cases, the equilibrium $x_0$ can be inferred directly from data by treating it as a trainable parameter during model training.\n\n\n> Does one need to retrain the model if the point needs to be changed?\n\nFor a specific, non-parametric dynamic system, the equilibrium point is fixed and does not change. As described above, it can be considered as a trainable parameter and thus be identified from data. \nHowever, for a parameterised system, also $x_0$ could depend on the parameter set $\\mu$. In ongoing work, we are extending the presented approach to parameterised systems, allowing parameter-dependent equilibrium points $x_0(\\mu)$, see also the next answer below.\n\n\n> Can this framework be extended to parametrized systems?\n\nYes. By including a set of parameters $p$ as inputs to the neural networks describing the Hamiltonian $\\mathcal{H}$ and the matrices $\\boldsymbol{J}$, $\\boldsymbol{R}$, and $\\boldsymbol{G}$, the parametric dependence of the dynamics can be learned from data. The equilibrium position $x_0$ can similarly be modelled as a neural network function of $p$.\nTo ensure sufficient flexibility and to avoid restricting $\\mathcal{H}(x, p)$ to be convex in $p$, partially input convex neural networks (PICNNs) [Amos et al., 2017] can be used to parameterise the Hamiltonian. We are currently developing this extension and could include it in the final version for one example, such as the spinning rigid body, which can be parameterised by its principal moments of inertia $I_1, I_2, I_3$ or a damping coefficient $\\mu$.\n\n\n> Why is the Lyapunov function introduced, if it is not required in the Theorem? It appears the function is only necessary in the proof?\n\nWe acknowledge that the dependency of the theorem on Lyapunov theory and especially Lyapunov functions is rather implicit and subtle. To address this, we will change line 155 from “Then, the system in Equation (5) has a stable equilibrium at $x(t) = 0$, and all solutions are bounded.” to “Then, the Hamiltonian $\\mathcal{H}$ is a suitable Lyapunov function for showing stability of the equilibrium at $x(t)=\\boldsymbol{0}$, and all solutions are bounded.”\nWith this change we aim to clarify the connection to Lyapunov theory and improve the intuition behind the theorem.\n\n\n> Is there another baseline comparison beyond deep learning models (that also incorporate classical methods like [structure-preserving] SINDy: [...])?\n\nWhile we have not made a direct comparison with structure-preserving SINDy yet, we expect that the method will struggle with complex dynamics in higher-dimensional settings (in the mentioned reference, the examples have only up to 3 states). In particular, ensuring a convex Hamiltonian or energy for achieving stability will be difficult with SINDy, as multivariate polynomials are generally non-convex (e.g., $x^2 y^2$ is not convex). In our opinion, input-convex NNs are much more versatile for this purpose."}, {"review_id": "Jgju3rCnkB", "reviewer": "Reviewer_1bj2", "summary": "This paper aims to incorporate physical prior into dynamical system learning. The authors introduce a framework named stable port-Hamiltonian neural networks. This framework aims to previous the physical biases while ensuring global Lyapunov stability. They also establish the requirement for stability of port-Hamiltonian systems and utilize this for model design. The authors evaluate the proposed method on different systems in comparison to NODEs, PHNNs, and bounded PHNNs (bPHNNs).", "questions": "Please see weaknesses.", "limitations": "Yes", "rating": 3, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 2, "clarity": 3, "significance": 3, "originality": 2, "strengths_and_weaknesses": "Strengths\n\n1. This paper is well organized and clearly written. \n\n2. The technical details are also easy to follow.\n\n3. The authors evaluate the proposed method on different systems to validate the effectiveness of the proposed method.\n\nWeaknesses\n\n1. The major concern is that the compared methods are out-of-date and limited. I do suggest that authors should include more state-of-the-art and up-to-date methods in performance comparison given that there are so many works for dynamical system modeling.\n\n2. It seems that all the data is generated by simulation. Are there any real-world datasets to validate the effectiveness of the proposed method? \n\n3. The difference between Stable Port-Hamiltonian Neural Networks and Stable NODEs should be carefully discussed to show the novelty of the proposed method. \n\n4. \"Why port-Hamiltonian dynamics?\" seems not to be the contribution of the paper. This part should move to the preliminary section instead.", "was_revised": true, "final_justification": "More baselines are needed.", "rebuttal": "We appreciate the reviewer’s thoughtful comments and the opportunity to clarify the scope and contributions of our work. Below, we address each point in detail, explaining our methodological choices and planned revisions to strengthen the manuscript.\n\n\n> The major concern is that the compared methods are out-of-date and limited. I do suggest that authors should include more state-of-the-art and up-to-date methods in performance comparison given that there are so many works for dynamical system modeling.\n\nTo include further baselines, we would appreciate it if the reviewer could provide more details on which methods are considered as state-of-the-art and should be included. \nCertainly, deep learning architectures such as transformers have shown great potential for time series forecasting in recent years, but are computationally demanding, data hungry, and requiring a lot of fine-tuning of the architectures (see, e.g., Chen et al., “A Closer Look at Transformers for Time Series Forecasting: Understanding Why They Work and Where They Struggle”). Thus, we focus on comparing to methods that model dynamic systems as ODEs, which is close to classical, physics-based modeling and allows to enforce properties such as thermodynamic consistency and stability, which we deem as highly beneficial for obtaining reliable predictions. In this context, we regard (stable) NODEs and energy-based concepts such as PHNNs (or similarly GENERIC, which exhibits similar performance, see Urdeitx et al., “A comparison of single and double generator formalisms for thermodynamics-informed neural networks”.) as reasonable baselines. \n\n\n> It seems that all the data is generated by simulation. Are there any real-world datasets to validate the effectiveness of the proposed method?\n\nThe data used in the cascaded tanks example in Section 4.2 originates from measurements of a real-world fluid level control system. This setup reflects an actual physical process and involves real sensor data. \n\nAlthough the data in Section 4.3 is generated through simulation, it is based on a high-fidelity multiphysics model that closely captures the physical processes of an actual thermal food processing system, where temperatures are measurable using standard temperature sensors. \n\nTo further address the potential impact of measurement noise, which may be absent from simulated data, we have included additional evaluations showing that the proposed method remains highly robust to random additive perturbations in the training data. For further details, please see our response to Reviewer n6dB.\n\n\n> The difference between Stable Port-Hamiltonian Neural Networks and Stable NODEs should be carefully discussed to show the novelty of the proposed method.\n\nThe main difference between our approach and stable NODEs (\"learning stable deep dynamics models\" [Kolter 2019]) or similar methods [Kojima 2022, Takeishi 2021, Okamoto 2024] is that our approach is projection-free. The named methods learn unconstrained nominal dynamics and then project them onto the subspace of stable systems, as determined by a concurrently learned Lyapunov function. In contrast, our approach leverages the port-Hamiltonian formalism to formulate the learning problem directly within the subspace of stable dynamics, eliminating the need for projection and avoiding associated numerical challenges. Additionally, the port-Hamiltonian formulation of our model offers interpretability (e.g., the split into conservative and dissipative dynamics) that other approaches lack. \n\nAs we show with the spinning rigid body experiment in Section 4.1 (see Fig. 2 and C.1), the sNODEs show stable behavior (they converge to 0 energy and angular velocities), but their accuracy and variance are much worse than the proposed sPHNNs.\n\nTo complement the discussion of the above points in the introduction, we will revise line 342 to read: “Additionally, the approach ensures global asymptotic stability of the identified dynamics without requiring projection, by constraining the Hamiltonian to be a convex, positive definite Lyapunov function.”\n\n\n> \"Why port-Hamiltonian dynamics?\" seems not to be the contribution of the paper. This part should move to the preliminary section instead.\n\nWe agree with the reviewer and will integrate Section 3.1 (“Why port-Hamiltonian dynamics”) into the background Section 2.2."}, {"review_id": "XqQLSAFnBb", "reviewer": "Reviewer_n6dB", "summary": "This paper discusses a type of physics guided neural network based on port-Hamiltonian systems, and in particular how global Lyapunov stability can be enforced in such a network. Lyapunov stability is desirable for some problems, because it is often a realistic inductive bias for physical systems, can improve extrapolation, and makes model training more robust.  The model works by representing each of the constituent terms of a port-Hamiltonian system (Hamiltonian, structure matrix, dissipation matrix, and input matrix) with a neural network. The key contribution seems to be the choice of architecture for the Hamiltonian network; by ensuring it is convex (and has a minimum of value 0 at the origin) global stability is enforced. The proposed model is validated on a number of different examples, both with real and synthetic data, and is shown to perform better than alternatives.", "questions": "- Please can you more clearly describe the contribution of the proposed model in relation to prior work? \n- Can you explain in more detail which real word problems this type of model can be applied to? How easy is it to determine if a real system is globally stable within the region of interest? How well does the model deal with e.g. noisy data?", "limitations": "yes", "rating": 5, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 4, "clarity": 3, "significance": 3, "originality": 2, "strengths_and_weaknesses": "## Strengths\n\n- The paper is technically sound to the best of my knowledge (Note: I have not checked the proofs in detail, and do not have a deep knowledge of the theory employed). The main claims made by the paper are adequately supported by the experiments. The experiments are standard for papers in this area.\n- The clarity of the paper is excellent. The writing is clear and easy to follow, and the mathematical exposition precise. The paper has a very polished feel, and the figures and plots are well designed and communicate the experimental results well. I enjoyed reading this paper!\n- While impact is difficult to judge, the proposed method seems to solve the problem in a more comprehensive way than past work (e.g OnsagerNet).\n- The paper is original in the sense that it offers a new combination of existing techniques.\n\n## Weaknesses \n\n- For me, the main weakness of this paper is the lack of clear explanation of the relationship of the proposed method to prior work. This makes it difficult for the reader to understand exactly what the contribution is. A number of related methods are discussed in the introduction, but it is unclear which of these proposed methods could be used to tackle the problems in the experiments section. The only existing method that seems to be used as comparison in the experiments is the sNODE, could other methods be applied here?\n- If I am correct in my understanding that the key contribution is constraining the form of the Hamiltonian network, then the contribution of the paper is modest in relation to prior work.", "was_revised": true, "final_justification": "Authors addressed all my concerns in response, and so I am raising my score to \"Accept\".", "rebuttal": "We appreciate the reviewer’s thorough evaluation and thoughtful comments. We are glad the clarity and rigour of the paper resonated well. Below, we aim to clarify the specific contributions of our approach, explain the choice of baselines in our experiments, and elaborate on how our method complements and extends existing techniques in the field. We also outline the planned changes in response to the review.\n\n\n> [...] The only existing method that seems to be used as comparison in the experiments is the sNODE, could other methods be applied here?\n\nWe want to clarify that the bounded PHNN (bPHNN) baseline model included in all our numerical experiments represents a version of OnsagerNet [Yu 2021], extended to support arbitrary time-dependent inputs to make them more versatile and suitable for our experimental setting.\n\n\n> Please can you more clearly describe the contribution of the proposed model in relation to prior work?\n\nWe can differentiate prior work that enforce some notion of stability into projection based methods such as stable NODEs (\"learning stable deep dynamics models\" [Kolter 2019]) or similar methods [Kojima 2022, Takeishi 2021, Okamoto 2024] and constrained methods for dissipative dynamics such as the OnsagerNet proposed by [Yu 2021]. The former learn unconstrained nominal dynamics and then project them onto the subspace of stable systems, as determined by a concurrently learned Lyapunov function. In contrast, our approach leverages the port-Hamiltonian formalism to formulate the learning problem directly within the subspace of stable dynamics, thereby eliminating the need for projection and avoiding associated numerical challenges, such as discontinuities in the dynamics resulting from the projection. These can hinder the training via trajectory fitting, which is required for modeling dynamics with augmented states. This also effectively excluded stable NODEs from most of our numerical examples, see also the discussion at the beginning of Section 4.\n\nIn contrast, constrained methods such as our stable port-Hamiltonian neural networks (sPHNNs) and OnsagerNet (or here bPHNNs)  formulate the learning problem directly in the subspace of stable dynamics, sidestepping the mentioned issues. \nWhile the structure of the evolution equation of the OnsagerNet is quite similar to sPHNNs, there are several key differences between both approaches:\n- **Inputs**: OnsagerNets do not consider arbitrary, time-dependent external input signals. This is possible with sPHNNs and is even required for most of the numerical examples we present.\n- **Energy sources**: OnsagerNets use state-dependent external forces to model autonomous systems with internal energy sources (necessary for systems with, e.g., limit-cycles or strange attractors). This introduces an ambiguity for decomposing the dynamics into the potential $V$ and the forcing term $f$, potentially weakening the induced physical bias and interpretability. In contrast, for sPHNNs, energy sources necessarily enter through the external inputs $u$, clearly distinguishing between internal dissipative dynamics and external influences.\n- **Stability**: OnsagerNets apply a potential that is bounded from below by zero and coercive (i.e., radially unbounded). Under certain conditions on the forcing term, OnsagerNet guarantees bounded solutions and at least one stable equilibrium, though multiple stable and unstable equilibria may exist. In contrast, sPHNNs enforce global asymptotic Lyapunov stability and guarantee a single equilibrium. While this is more restrictive, it introduces a stronger inductive bias, leading to improved performance in applicable cases, particularly in the small-data regime.\n\n\n> Can you explain in more detail which real word problems this type of model can be applied to? \n\nThe proposed model is based on the port-Hamiltonian framework, which is widely used to describe physical systems in various fields, including mechanics, electromechanics, quantum mechanics, control theory, acoustics, fluid dynamics, and thermodynamics. This broad applicability carries over to our model, as demonstrated by our experiments for problems in rigid body dynamics, fluid mechanics, conjugate heat transfer, and thermal diffusion.\n\nSpecifically, the stable port-Hamiltonian neural networks (sPHNNs) we introduce can be applied to identifying the dynamics of such systems, provided the system is globally stable. While global stability is a strong assumption, our numerical examples show that it holds in many practically relevant scenarios. In such cases, sPHNNs offer a strong inductive bias, which can significantly improve learning efficiency and generalisation compared to more general-purpose models.\n\nIn practical applications, we envision sPHNNs being used as efficient yet accurate surrogate models for complex simulations, particularly in engineering contexts where stability and physical consistency are critical. We demonstrate this applicability in Sections 4.2-4.4.\n\n\n> How easy is it to determine if a real system is globally stable within the region of interest? \n\nWithout prior knowledge, determining the stability properties of a given system is a challenging task. While in many cases, even basic physical insights can be helpful, such knowledge may not always be readily available. Under these circumstances, we see our model as a first step in an iterative process: after training an sPHNN, its performance can indicate the validity of the global stability assumption. This approach is low-cost, requiring minimal data, and is fast due to the model’s small network sizes. If performance is insufficient, one can then iteratively explore more general models with weaker inductive biases such as bPHNNs, NODEs, etc..\n\n\n> How well does the model deal with e.g. noisy data?\n\nIn our numerical experiments, the model demonstrated strong robustness due to the incorporation of the global stability constraint. This was especially evident with respect to variations in the initial trainable parameters and the quantity of training data. While the real-world measurement data employed in Section 4.2 is noisy, it has a relatively large signal-to-noise ratio. \n\nTo more directly assess the model’s robustness to noise, we conducted an additional experiment based on the setup presented in Section 4.3. We have trained the models using $n_A=3$ augmented dimensions and $n_D=2$ training trajectories, introducing zero-mean Gaussian noise for both the input $T_{oven}$ and the outputs $T_A$ and $T_B$. The noise amplitudes were varied between 5% and 25% of the original signals’ standard deviations. The table below reports the average root mean squared error (RMSE) over 20 model instances per model type, evaluated on noise-free test data.\n\n*Table 1: Test RMSE*\n| Model   |   5% noise |   15% noise |   25% noise |\n|:--------|-----------:|------------:|------------:|\n| sPHNN   |      1.170 |       2.595 |       3.177 |\n| bPHNN   |      1.572 |       2.828 |       5.231 |\n| PHNN    |     28.233 |      25.631 |      24.224 |\n| NODE    |      3.577 |       3.482 |       4.362 |\n\nTo isolate the effect of noise and eliminate generalisation error, we also evaluated the RMSE of the models trained on noisy data on the noise-free training data:\n\n*Table 2: Training RMSE*\n| Model   |   5% noise |   15% noise |   25% noise |\n|:--------|-----------:|------------:|------------:|\n| sPHNN   |      0.332 |       0.854 |       1.354 |\n| bPHNN   |      0.500 |       1.049 |       2.207 |\n| PHNN    |      2.828 |       6.750 |       5.779 |\n| NODE    |      0.993 |       1.974 |       2.423 |\n\nWhile all models demonstrate a degree of robustness to noise, the sPHNN consistently outperforms the baselines on both the training and test sets. This robustness can be partially attributed to the training process: the integration step involved in trajectory fitting inherently acts as a low-pass filter, effectively smoothing the learned dynamics. As a result, high-frequency Gaussian noise is attenuated, which helps prevent overfitting and improves generalisation.\n\nWe will incorporate the additional results and corresponding discussion in Section 4.3 of the revised manuscript."}, {"review_id": "TMnMLbDOyP", "reviewer": "Reviewer_BVE6", "summary": "This paper tackles the problem of learning stable dynamical systems from data using neural networks. Building on the framework of Lyapunov stability and port-Hamiltonian systems, the authors propose a novel architecture, stable port-Hamiltonian neural networks (sPHNN), in which the system dynamics are modeled using port-Hamiltonian structure. Neural networks are used to learn an energy function and structured matrices representing energy-conserving flows, dissipative effects, and control inputs. The architecture enforces Lyapunov conditions on the learned energy function, thereby guaranteeing asymptotic stability. The proposed model is evaluated on a set of experiments, including both synthetic and real-world datasets, and demonstrates improved stability and data efficiency compared to several baselines.", "questions": "1. In Theorem 3.1, isn't coercity for a convex function admitting an uniquer minimizer a standard result in convex analysis? See Rockafeller's Convex Analysis book, Corollary 8.7.1. \n\n2. In line 93, the derivative $\\dot{V}(x)$ is scalar-valued. Should it not be referred to as “nonpositive” rather than “negative semi-definite”? Also, I found the second member of the equality chain (2) confusing as $\\dot{x}$  is a velocity, implying a trajectory, while the definition $\\dot{V}(x)$ is a function of the state and can be defined as a function of $f$, $V$ and $x$ regardless of any trajectory or dynamical system.\n\n3. How is the prediction task evaluated? Are the models used autoregressively, with known or fixed control inputs $u$? More details on this setup would help interpret the results.\n\n4.   If the system is known to admit a Lyapunov function vanishing at $x_0 = 0$, can we expect the learned Hamiltonian to converge to this function under the training dynamics?\n\n5.   What is the interpretation or expected benefit of applying a Hamiltonian structure to systems that are not inherently physical or do not follow Hamiltonian dynamics (e.g., temperature fields)? Could the same stability guarantees not be obtained more simply by enforcing a global Lyapunov function, as done in sNODEs?\n\n6.   The proposed architecture involves gradient computations of the learned Hamiltonian, which may incur additional cost. How does the computational complexity compare to standard NODEs or other baseline methods?", "limitations": "The submission does not discuss some important potential limitations such as the computational cost. See my question 6.", "rating": 5, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 4, "clarity": 4, "significance": 3, "originality": 3, "strengths_and_weaknesses": "Overall, the paper is very pleasant to read. In spite of minor clairty issues, the method is soundly and clearly presented, and the results are strong.\n \n### Strengths\n\n- The paper is clearly written and mathematically well-grounded. The theoretical framework is well-motivated, and all necessary background is introduced with appropriate references.\n- The proposed approach addresses an important and practical problem: learning dynamical models that are both data-efficient and stable, which is particularly relevant in low-data regimes.\n- The architecture is original and benefits from strong theoretical foundations. Structuring the model around a learned energy function with interpretable terms (conservative and dissipative) adds physical meaning and interpretability.\n- The paper provides a clear and detailed exposition of the model, with Theorem 3.1 nicely summarizing the desirable properties of the architecture.\n- The baseline methods are well-documented and chosen appropriately, enabling a fair comparison.\n- Experimental results on real-world and high-dimensional synthetic systems show a clear advantage of the proposed method over alternatives.\n\n### Weaknesses\n\n  1. It would be insightful to provide early in the paper an example or figure of a physical system illustrating the need for stability, and what an unconstrained architecture might produce. It would also be nice to introduce an example of port-Hamiltonian system to illustrate the framework of section 2. Typically, a torque-actuated pendulum with friction, where all the terms have a clear physical interpretation, would help ground the abstract formulation of (3). Perhaps a slightly more complex system might also be interesting to expose. \n\n2.   Although the model enforces structure via an interpretable decomposition into conservative and dissipative terms, this structure is defined with respect to a learned energy function. This function is optimized jointly with other parameters to minimize prediction error, and does not necessarily align with any known physical energy.\nIn the case where the learned system is known to be close to be Hamiltonian, and some prior knowledge on some energy function of the system is available, I am understanding that it could be incorporated in function $f$, and (3) then mimics the system dynamics. When it is not the case though, I am wondering what the interpretation of modeling the dynamics with (3) is. Why should we expect this structure to learn more efficiently, compared to alternative Lyapounov-based learning architectures? See my questions 4. and 5..\n\n3. Evaluation methodology needs clarification.  The experimental setup, particularly how predictions are computed and evaluated, lacks clarity. See question 3.\n\n4.  The paper does not discuss potential downsides of the proposed approach, especially regarding computational cost. The gradient computations involved in the Hamiltonian structure may significantly increase training time.", "was_revised": true, "final_justification": "The proposed method is solid and the paper is clearly written. The authors properly addressed my questions. I maintain my positive score of 5.", "rebuttal": "We sincerely thank the reviewer for their thoughtful and constructive feedback, as well as the considerable effort invested in this review. The comments were particularly insightful and raised important points that have helped us improve both the clarity and depth of the manuscript. Below, we respond to each suggestion in detail and outline the corresponding revisions we plan to implement.\n\n\n> It would be insightful to provide early in the paper an example or figure of a physical system [...]\n\nWe thank the reviewer for the suggestion and fully agree that such an early example will enhance the accessibility and motivation. We propose to move the spinning rigid body example, especially Equation C.3 and Figure C.1, before Section 3. This example provides a physically interpretable system within the port-Hamiltonian framework and clearly illustrates the benefits of incorporating physically motivated biases when identifying its dynamics.\n\nWhile learning the dynamics of simple systems using unconstrained methods often does not develop the same level of instability as in more practically relevant scenarios (see, e.g., Sections 4.3 and 4.4), we believe that this example would still be insightful and engaging. In particular, it would help demonstrate the interpretability advantages afforded by the port-Hamiltonian structure.\n\n\n> In Theorem 3.1, isn't coercity for a convex function admitting an uniquer minimizer a standard result in convex analysis? See Rockafeller's Convex Analysis book, Corollary 8.7.1.\n\nWhile the result might be standard, we’ve still included it in our proof in the appendix to aid in the intuitive, geometric understanding. We will add a suitable reference in the revision.\n\n\n> In line 93, the derivative $\\dot V(x)$ is scalar-valued. Should it not be referred to as “nonpositive” rather than “negative semi-definite”? Also, I found the second member of the equality chain (2) confusing as $\\dot{x}$ is a velocity, implying a trajectory, while the definition $\\dot{V}(x)$ is a function of the state and can be defined as a function of $f$, $V$ and $x$, regardless of any trajectory or dynamical system.\n\nNegative semi-definiteness refers here to the definition of definiteness of a function, which is related to, but distinct from, the definiteness of matrices. Furthermore, the notation $\\dot{V}(x)$ denotes the change of the function $V$ along the trajectory $x(t)$, which reduces Equation (2) simply to the application of the chain rule. This notation is standard in the literature about Lyapunov theory (see, for example, H. K. Khalil “Nonlinear systems”, or F. Verhulst “Nonlinear differential equations and dynamical systems”). However, to the best of our knowledge, it is rarely used outside of stability theory.\n\n\n> How is the prediction task evaluated? Are the models used autoregressively, with known or fixed control inputs $u$? More details on this setup would help interpret the results.\n\nAll models take the form of an ODE $\\dot{x}=f(x,u(t))$, and predictions are generated by numerically integrating this system forward in time using a known input function $u(t)$. This is not an autoregressive setup in the conventional sense, as the models predict state derivatives rather than directly outputting the next state. For our evaluation, the input trajectory is provided beforehand for the duration of the prediction. Nevertheless, the model can be used in a setting where $u(t)$ is supplied in real time, without requiring access to the whole input trajectory in advance.\n\nTo clarify this in the manuscript, we will revise line 246 to:\n“Model predictions are generated by rollout with a numerical integrator using the Runge-Kutta scheme Tsit5, with adaptive step size and a known input function $u(t)$.”\n\n\n> If the system is known to admit a Lyapunov function vanishing at x_0=0, can we expect the learned Hamiltonian to converge to this function under the training dynamics?\n\nLyapunov functions are generally not unique, so one should not expect the learned Hamiltonian to converge to any particular Lyapunov function, even if one is known for the true system. Importantly, for establishing stability properties, only the existence of a Lyapunov function matters - not its specific form. Our model enforces this by constraining the Hamiltonian always to satisfy the Lyapunov conditions, thereby ensuring global stability regardless of whether it matches a known energy function.\nHowever, since the Hamiltonian not only acts as a Lyapunov function but also generates the port-Hamiltonian dynamics, the best predictions are achieved if the learned Hamiltonian resembles the true system's energy.\n\n\n> What is the interpretation or expected benefit of applying a Hamiltonian structure to systems that are not inherently physical or do not follow Hamiltonian dynamics (e.g., temperature fields)? \n\nIndeed, systems like temperature dynamics are not inherently Hamiltonian in the classical sense, and non-physical systems may lack a physical energy. In such cases, our use of a port-Hamiltonian structure is not necessarily intended to recover a physical Hamiltonian, but rather to enforce a principled, structured decomposition of the dynamics into conservative and dissipative components. \nCompared to models without structure, this decomposition offers interpretability and modularity. Even without physical correspondence, the decomposition can yield meaningful insight into the system’s qualitative behaviour (see, e.g., Figure 1b). Conversely, knowledge about the quantitative behaviour of the system can be integrated into the individual components.\n\n\n> Could the same stability guarantees not be obtained more simply by enforcing a global Lyapunov function, as done in sNODEs?\n\nWhile sNODEs also provide a stability guarantee, we argue this is not achieved more simply, neither in computational terms nor in practical application. Our reasoning is twofold:\n- **Computational complexity**: sNODEs require numerically evaluating gradients of a learned Lyapunov function and enforcing stability via a projection step. This offers no computational advantage over sPHNNs. In fact, for our numerical example 4.1, sPHNNs were faster to train and to evaluate than sNODEs (see tables below). While we did not focus on runtime optimisation (as the models trained in <2 minutes), these results suggest that sPHNNs are at least as efficient, if not simpler, in practice.\n- **Numerics** The projection step in sNODEs introduces discontinuities in the dynamics (see the discussion at the beginning of Section 4 and in the introduction). These discontinuities can make the dynamics more challenging to integrate numerically, often causing adaptive-step integrators to take smaller steps and, in turn, requiring more model evaluations. Furthermore, in combination with trajectory fitting, we observed failure at training convergence for sNODEs, which excluded them from most of our numerical examples. In contrast, methods like our sPHNN, which don't rely on a projection to guarantee stability, do not suffer from these issues. In cases where the target dynamics are not inherently physical or do not follow Hamiltonian dynamics, one can therefore view the port-Hamiltonian formulation of sPHNNs simply as a way of achieving projection-free stability guarantees.\n\nIn summary, even for non-physical systems, sPHNNs offer a structured, projection-free framework for learning stable dynamics, combining theoretical guarantees with practical advantages.\n\n\n> The proposed architecture involves gradient computations of the learned Hamiltonian, which may incur additional cost. How does the computational complexity compare to standard NODEs or other baseline methods?\n\nWe agree that a more detailed comparison of the computational performance will enhance the paper.\n\nFor training the models, we measured the following average training times:\n\n*Table 1: Representative training times*\n| Training Time (min)              |  sPHNN |  bPHNN |   PHNN |   NODE | sNODE |\n|----------------------------------|-------:|-------:|-------:|-------:|------:|\n| Derivative fitting (Section 4.1) |  0.388 |  1.291 |  0.346 |  0.223 | 1.456 |\n| Trajectory fitting (Section 4.3) | 37.663 | 38.329 | 19.111 | 17.951 |     - |\n\nFor inference, we recorded the following times using data from Section 4.3:\n\n*Table 2: Representative evaluation times*\n| Evaluation time (ms)   |   sPHNN |   bPHNN |   PHNN |   NODE |\n|:-----------------------|--------:|--------:|-------:|-------:|\n| Derivative evaluation  |   0.006 |   0.006 |  0.006 |  0.005 |\n| Model integration      |   0.355 |   0.410 |  0.228 |  0.219 |\n\nHere, “Derivative evaluation” refers to the average time taken by each model to compute the state derivative $\\dot{x}$ for a given state $x$ and input $u$, which needs to be evaluated for time integration. In contrast, “Model integration” denotes the time required to generate a single trajectory prediction via integration. While the latter better reflects practical use cases, the reported integration time also depends on the integration scheme and chosen step size. Since we use an adaptive step size controller, the step size in turn depends on the learned dynamics, which explains the evaluation differences between sPHNN and PHNN: the latter predicted wrong but smooth dynamics, which may be faster to integrate.\n\nIn summary, while the proposed sPHNN introduces some additional computational cost relative to unconstrained NODEs, it performs similarly to or better than baseline models that also provide stability guarantees. Importantly, all models exhibit short training times compared to conventional deep learning approaches, due to their ability to capture complex dynamics with relatively small neural networks. Consequently, computational efficiency typically does not present a bottleneck in practical scenarios.\n\nWe will include these results and the accompanying discussion in the appendix of the revised manuscript."}]

Author rebuttal:

Meta-review: Stable port-Hamiltonian neural networks enforce global Lyapunov stability via a convex Hamiltonian parameterization, combining physical interpretability with inductive bias. Strengths are a clear theoretical foundation, projection-free stability guarantees, and solid empirical results. After an extensive discussion period, the authors managed to convince three reviewers fully, who are now clearly in favor of acceptance, highlighting clarity, rigor, and practical relevance. One reviewer questioned the choice of baselines and novelty but did not engage further after the authors’ detailed clarifications, which I find convincing. I recommend acceptance as a poster.
Decision: Accept (poster)