Stable Port-Hamiltonian Neural Networks

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.

It would be insightful to provide early in the paper an example or figure of a physical system [...]

We 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.

While 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.

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.

While 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.

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.

Negative 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.

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.

All 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.

To clarify this in the manuscript, we will revise line 246 to: “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)$.”

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?

Lyapunov 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. However, 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.

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)?

Indeed, 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. Compared 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.

Could the same stability guarantees not be obtained more simply by enforcing a global Lyapunov function, as done in sNODEs?

While 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:

• 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.
• 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.

In summary, even for non-physical systems, sPHNNs offer a structured, projection-free framework for learning stable dynamics, combining theoretical guarantees with practical advantages.

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?

We agree that a more detailed comparison of the computational performance will enhance the paper.

For training the models, we measured the following average training times:

Table 1: Representative training times

For inference, we recorded the following times using data from Section 4.3:

Table 2: Representative evaluation times

Here, “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.

In 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.

We will include these results and the accompanying discussion in the appendix of the revised manuscript.

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.

[...] The only existing method that seems to be used as comparison in the experiments is the sNODE, could other methods be applied here?

We 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.

Please can you more clearly describe the contribution of the proposed model in relation to prior work?

We 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.

In 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. While the structure of the evolution equation of the OnsagerNet is quite similar to sPHNNs, there are several key differences between both approaches:

• 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.
• 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.
• 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.

Can you explain in more detail which real word problems this type of model can be applied to?

The 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.

Specifically, 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.

In 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.

How easy is it to determine if a real system is globally stable within the region of interest?

Without 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..

How well does the model deal with e.g. noisy data?

In 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.

To 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.

Table 1: Test RMSE

To 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:

Table 2: Training RMSE

While 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.

We will incorporate the additional results and corresponding discussion in Section 4.3 of the revised manuscript.

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.

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.

To 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. Certainly, 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.

It seems that all the data is generated by simulation. Are there any real-world datasets to validate the effectiveness of the proposed method?

The 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.

Although 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.

To 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.

The difference between Stable Port-Hamiltonian Neural Networks and Stable NODEs should be carefully discussed to show the novelty of the proposed method.

The 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.

As 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.

To 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.”

"Why port-Hamiltonian dynamics?" seems not to be the contribution of the paper. This part should move to the preliminary section instead.

We agree with the reviewer and will integrate Section 3.1 (“Why port-Hamiltonian dynamics”) into the background Section 2.2.

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.

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?

The proposed model incorporates the equilibrium position $x_0$​ directly into its architecture, but its value does not need to be known in advance. In 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. While 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.

To 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.

Does one need to retrain the model if the point needs to be changed?

For 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. However, 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.

Can this framework be extended to parametrized systems?

Yes. 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$. To 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$.

Why is the Lyapunov function introduced, if it is not required in the Theorem? It appears the function is only necessary in the proof?

We 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.” With this change we aim to clarify the connection to Lyapunov theory and improve the intuition behind the theorem.

Is there another baseline comparison beyond deep learning models (that also incorporate classical methods like [structure-preserving] SINDy: [...])?

While 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.