Paper deep dive
Safe-by-Design Learning via Energy-based Neural Networks
Simone Betteti, Morteza Lahijanian, Luca Laurenti
Intelligence
Status: not_run | Model: - | Prompt: - | Confidence: 0%
Entities (0)
Relation Signals (0)
No relation signals yet.
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:Learning neural-network models of dynamical systems with safety guarantees is a fundamental requirement for their deployment in safety-critical settings. Safety is commonly established by proving the invariance of a desired subset in state-space, ensuring that every trajectory initialized in this subset remains confined to it for all time under admissible inputs. Existing frameworks, however, either rely on computationally expensive post-hoc verification or employ safety-enforcing mechanisms without formal correctness guarantees. In this paper, we introduce a novel neural architecture grounded in energy-based modern Hopfield networks to guarantee safety-by-design while retaining sufficient expressiveness to model complex nonlinear dynamics. Specifically, we integrate modern Hopfield networks with a port-Hamiltonian neural ODE, enabling by design the construction of barrier functions yielding explicit admissible-input sets and quantitative robustness radii. Across several benchmarks, including an 12-dimensional nanodrone model, our framework achieves state-of-the-art performance while producing certified invariant sets that are more robust to external solicitations than comparable existing approaches.
Tags
Links
- Source: https://arxiv.org/abs/2609.36942v1
- Canonical: https://arxiv.org/abs/2609.36942v1
Trouble viewing inline? Open PDF directly →
Full Text
131,553 characters extracted from source content.
Expand or collapse full text
Safe-by-Design Learning via Energy-based Neural Networks Simone Betteti Affiliation: RIAS Lab Affiliation: The Italian Institute of Artificial Intelligence for Industry Affiliation: Torino, 10129, ITA Email: simone.betteti@ai4i.it Morteza Lahijanian Affiliation: RECUV Affiliation: University of Colorado Boulder Affiliation: Boulder, CO 80303, USA Email: morteza.lahijanian@colorado.edu Luca Laurenti †thanks: L. Laurenti is also with the Delft Center for Systems and Control, TU Delft, Delft, 2600 A, The Netherlands. Affiliation: RIAS Lab Affiliation: The Italian Institute of Artificial Intelligence for Industry Affiliation: Torino, 10129, ITA Email: luca.laurenti@ai4i.it Abstract Learning neural-network models of dynamical systems with safety guarantees is a fundamental requirement for their deployment in safety-critical settings. Safety is commonly established by proving the invariance of a desired subset in state-space, ensuring that every trajectory initialized in this subset remains confined to it for all time under admissible inputs. Existing frameworks, however, either rely on computationally expensive post-hoc verification or employ safety-enforcing mechanisms without formal correctness guarantees. In this paper, we introduce a novel neural architecture grounded in energy-based modern Hopfield networks to guarantee safety-by-design while retaining sufficient expressiveness to model complex nonlinear dynamics. Specifically, we integrate modern Hopfield networks with a port-Hamiltonian neural ODE, enabling by design the construction of barrier functions yielding explicit admissible-input sets and quantitative robustness radii. Across several benchmarks, including an 12-dimensional nanodrone model, our framework achieves state-of-the-art performance while producing certified invariant sets that are more robust to external solicitations than comparable existing approaches. 1 Introduction Learning dynamical systems from data is increasingly central to robotics, autonomous systems, and scientific modeling, where neural state-space models can reconstruct complex nonlinear trajectories with remarkable accuracy (Karniadakis et al., 2021; Legaard et al., 2023). In safety-critical applications, however, predictive accuracy over observed trajectories is not sufficient. A model that performs well over the finite horizons represented during training may develop unstable or physically undesirable behavior when propagated under new inputs or for longer time intervals. Safety is particularly important when learned dynamics are intended for applications such as flight, autonomous operation, or human-robot interaction, where the behavior of the model outside densely sampled training regions is itself operationally relevant (Ashmore et al., 2021; Brunke et al., 2022; Gyevnár and Kasirzadeh, 2025). The central question is therefore not only whether a neural model can accurately identify nonlinear dynamics, but whether useful properties of its future evolution can be certified before deployment. Safety in learning-enabled dynamical systems is commonly enforced through external mechanisms on top of an unconstrained predictor, including safety filters, runtime monitors, constrained controllers, and post-hoc verification (Brunke et al., 2022; Hewing et al., 2020). Certifiability mechanisms are indispensable components of safety-critical systems, but they leave open a complementary question: can certifiability be made an intrinsic property of the learned dynamics themselves? Neural ordinary differential equations (NODEs) (Chen et al., 2018; Greydanus et al., 2019) have motivated several approaches towards safety-certified dynamics, using Lyapunov functions, dissipativity constraints, or algebraic structure to control behavior beyond the training trajectories (Kolter and Manek, 2019; Lawrence et al., 2020; Kojima and Okamoto, 2022; Kang et al., 2021; Yang et al., 2022). Port-Hamiltonian neural networks are particularly attractive because they decompose the vector field into energy-preserving interconnection, energy dissipation, and external actuation through explicit input ports (Van der Schaft, 2007; Massaroli et al., 2020). Their energy balance provides a natural structural primitive for stability and control without requiring the learned vector field to remain unconstrained. A central tension nevertheless remains between certifiability and expressivity. Strong stability guarantees in recent port-Hamiltonian architectures are commonly obtained through restrictive Hamiltonian parameterizations, including convex energies designed around a single equilibrium (Roth et al., 2025). Yet nonconvex energy landscapes are precisely what allow energy-based models to represent multiple stable regimes and complex nonlinear organization. Modern Hopfield networks provide a particularly expressive realization of this principle: advances in energy-based architectures have shown that rich, nonconvex landscapes can coexist with an explicit energy structure (Krotov and Hopfield, 2020; Hoover et al., 2022; Hoover et al., 2023). For controlled physical dynamics, this suggests a different route to safe learning: rather than simplifying the Hamiltonian until certification becomes tractable, preserve an expressive energy landscape and exploit its geometry directly to determine which external inputs are compatible with safe evolution. We pursue this approach through a port-Hamiltonian energy-based model (pH-EBM), in which a radially unbounded modern Hopfield energy (Fig. 1b) parametrizes the Hamiltonian of a controlled neural ODE (Fig. 1a). Radial unboundedness confines bounded-energy trajectories to compact sublevel sets, while the interior Hamiltonian remains free to develop nonconvex and multi-well structure. The port-Hamiltonian decomposition provides an explicit energy balance between dissipation and external actuation. We then use the same learned Hamiltonian to define energy barrier functions: for a prescribed sublevel component, the barrier condition identifies the external inputs that cannot inject sufficient energy to drive the learned trajectory across its boundary (Fig. 1c). Safety is therefore not introduced as an additional network or post-processing stage, but derived from the energy geometry already governing the neural dynamics. The resulting guarantees concern invariance of prescribed regions of the learned model state space under certified classes of inputs; when bounded disturbances or model mismatch are explicitly characterized, the same construction extends to robust certificates. Our contributions are threefold. First, we introduce a pH-EBM that combines the dissipative structure of port-Hamiltonian NODEs with a coercive, globally nonconvex modern-Hopfield Hamiltonian, separating boundedness at large state norm from the geometry required to represent complex operating regimes. Second, we derive explicit set-valued safety certificates from the learned energy, such as exact state-dependent admissible-input sets, state-independent robustness radii for entire energy regions, and robust extensions under bounded perturbations. We further show how the certified margin is governed by the geometry of the energy shell, dissipation normal to its boundary, and exposure of this direction to the input port, with local curvature providing tractable closed-form bounds. Third, we evaluate the framework across nonlinear system-identification and controlled-mechanics benchmarks spanning classical physical systems, nonconvex multi-attractor dynamics, adversarial out-of-distribution forcing, and high-dimensional flight identification. The experiments show that pH-EBM can offer competitive predictive accuracy, often outperforming state-of-the-art methods, while providing non-trivial safety guarantees by design. Figure 1: Safe-by-design Energy-Based Model. a) Input-output trajectories from physical or simulated systems are used to train the pH-EBM: the external input u drives the neural dynamics, while the measured output y supervises reconstruction. b) The learned Hamiltonian, parameterized by a modern Hopfield energy, organizes the neural ODE into a continuous energy landscape whose low-energy regions encode stable operating regimes. Beyond interpolating between sparsely observed trajectories, this geometry identifies energy sublevel sets that can be certified as invariant for prescribed classes of inputs. c) Inputs lying within the certified admissible set therefore generate trajectories that remain confined to the desired operating region, turning the learned energy landscape into an explicit quantitative safety certificate. 2 Problem formulation We consider the learning of an unknown controlled continuous-time dynamical system described by x˙⋆(t) x (t) =f⋆(x⋆(t),u(t)),y(t)=h⋆(x⋆(t),u(t))+v(t), =f (x (t),u(t) ), y(t)=h (x (t),u(t) )+v(t), (1) where x⋆(t)∈n⋆x (t)∈ ^n_ is the generally unobserved physical state, u(t)∈⊆ku(t) ^k is a known bounded external input, y(t)∈py(t)∈ ^p is the measured output, and v denote measurement perturbations. Throughout the paper, we assume that f⋆f satisfies the standard Carathéodory conditions, i.e., locally Lipschitz continuous in x and locally essentially bounded over the admissible values of u. These assumptions guarantee the existence and uniqueness of forward-complete Carathéodory solutions (Khalil, 2002). Both the physical state x⋆x (including its dimensionality) and the vector field f⋆f are unknown, and the system can be observed only through sampled inputs u and outputs y. The available dataset therefore consists of M experiments, possibly collected from different unknown initial conditions and under different input signals, i.e., :=(u(j)(ti),y(j)(ti))i=0Njj=1M.D:= \ (u^(j)(t_i),y^(j)(t_i) )_i=0^N_j \_j=1^M. (2) Since the physical state x⋆x is unobserved, our objective is not necessarily to recover its original realization, but to learn a latent state-space model that reproduces the observed input-output behaviour and generalizes within a neighbourhood sufficiently covered by the data. Specifically, we seek a continuous-time neural state-space model of the form z˙(t) z(t) =fθ(z(t),u(t)),y^(t)=hϕ(z(t),u(t)), =f_ (z(t),u(t) ), y(t)=h_ (z(t),u(t) ), (3) where z(t)∈dz(t)∈ ^d is the latent model state, θ and ϕ are trainable parameters, hϕh_ is a continuous latent-to-output decoder function, and y^(t)∈p y(t)∈ ^p is the measurement of the latent state z(t)z(t). Note that y y is in the same space as y of System (1). The initial latent state associated with each experiment is inferred jointly with the model parameters, or produced by an appropriate encoder. Our objective is not merely to fit the trajectories in D. Instead, we seek a safe-by-design neural architecture. To give this notion a physical interpretation, let safe⊆pY_safe ^p be a prescribed closed set of outputs that satisfy the relevant operational safety requirements. Safety of the learned model is established by jointly constructing a robustly invariant latent set safeZ_safe whose corresponding outputs remain in safeY_safe. Definition 1 (ρ -Robust Invariant Sets). Let ρ>0 >0 be an input-intensity threshold and define the associated input set u∈ρ:=u∈:‖u‖≤ρu _ :=\u :\>\|u\|≤ \. We say that the compact set of latent states ⊂ℝdC ^d is ρ -Robust invariant if: • for all z∈z and u∈ρu _ it holds that hϕ(z,u)∈safeh_ (z,u) _safe and • for all z(0)∈z(0) and input signals u such that that u(t)∈ρu(t) _ for all t>0t>0, it holds that z(t)∈z(t) for all t>0t>0. In this paper, we will establish results for the latent representations z, and then exploit the continuity of the decoder map hϕh_ to extend safety to compacts in measurement space. Importantly, in our setting the geometry of safeZ_safe is generally unknown before training. Since the latent representation z, its dynamics fθf_ , and the output map hΦh_ are learned jointly from data, an invariant set in latent space safeZ_safe cannot generally be prescribed a priori. It must instead be learned jointly with the neural state-space realization. This leads to the following problem: Problem 1. Learning and certifying neural ODEs. Given a dataset D, a set safe,Z_safe, and an input-intensity parameter ρ>0 >0, learn a neural ODE z˙(t)=fθ(z(t),u(t)) z(t)=f_ (z(t),u(t)) with measured output y^(t)=hϕ(z(t),u(t)) y(t)=h_ (z(t),u(t)) such that for all u(t)∈ρt≥0\u(t) _ \_t≥ 0 returns a maximal ρ -Robust Invariant Set ⊆safeC _safe. Many autonomous systems admit nonempty positively invariant regions, while systems operating in closed loop with stabilising feedback often admit nonempty robust positively invariant regions under bounded disturbances (Blanchini, 1999; Blanchini and Miani, 2015). The objective of Problem 1 is therefore to endow an input-output model learned from data with this structural property: we seek to identify a latent state-space realization while simultaneously guaranteeing that a nontrivial latent safe set is robustly positively invariant. Solving Problem 1 is, however, particularly challenging because generic post-hoc methods for computing or certifying invariant sets for neural dynamics, based for instance on gridding, mixed-integer programming, or SMT solving, are computationally demanding and scale poorly with the state dimension and network size (Abate et al., 2021; Dawson et al., 2023). We therefore pursue a complementary route and consider neural architectures that admit a guaranteed ρ -robust invariant set by construction. In particular, in Section 3 we introduce a port-Hamiltonian neural architecture whose Hamiltonian sublevel components provide barrier certificates under explicit dissipation and input conditions. The barrier certificates ultimately allow us to characterize safe ρ -robustly forward-invariant sets. In Section 4 we instead present the performances of our pH-EBM model in complex dynamics learning tasks, displaying how the model achieves state-of-the-art accuracy and provides explicit safety certificates, using drones’ bounded long time rollouts as a paradigmatic proof of concept. 3 Safe-by-design port-Hamiltonian neural ODEs The solution of Problem 1 with a neural ODE prompts the choice of a vector field structure that simultaneously combines two key ingredients: (i) a stability-inducing mechanism to confine trajectories and (i) an input port. Thereby, following (Roth et al., 2025), we specialize the neural state-space model in (3) to the port-Hamiltonian vector field z˙=fθ(z,u)=[JθJ(z)−RθR(z)]∇ℋθH(z)+GθG(z)u, z=f_ (z,u)= [J_ _J(z)-R_ _R(z) ] _ _H(z)+G_ _G(z)u, (4) where z∈dz∈ ^d, u∈⊆ku ^k, and ℋθH:d→H_ _H: ^d→ is a continuously differentiable Hamiltonian. Moreover, JθJ(z)=−JθJ(z)⊤∈d×dJ_ _J(z)=-J_ _J(z) ∈ ^d× d is skew-symmetric, 0⪯RθR(z)=RθR(z)⊤∈d×d0 R_ _R(z)=R_ _R(z) ∈ ^d× d is symmetric positive semidefinite, and GθG(z)∈d×kG_ _G(z)∈ ^d× k is the input map. The symmetry constraints on JθJ_ _J and RθR_ _R are enforced directly through their parameterizations. We assume that JθJ_ _J, RθR_ _R, GθG_ _G, and ∇ℋθH _ _H are locally Lipschitz in z and that u:[0,∞)→u:[0,∞) is locally essentially bounded. Consequently, (4) admits a unique Carathéodory solution on its maximal interval of existence. Equation 4 is a standard port-Hamiltonian ODE and satisfies the energy balance (Van der Schaft, 2007; Willems, 2007) dtℋθH(z) ddtH_ _H(z) =−∇ℋθH(z)⊤RθR(z)∇ℋθH(z)+∇ℋθH(z)⊤GθG(z)u. =- _ _H(z) R_ _R(z) _ _H(z)+ _ _H(z) G_ _G(z)u. (5) Thus, skew-symmetry of JθJ_ _J makes the interconnection term energy preserving, positive semidefiniteness of RθR_ _R makes the internal dynamics dissipative, and GθGuG_ _Gu accounts for energy exchanged with the environment. Importantly, this balance does not require convexity of the Hamiltonian. Consequently, we now define a non-convex Hamiltonian, which provides the key architectural ingredient for constructing compact confinement regions in the latent space. In particular, we employ the hybrid modern Hopfield Hamiltonian introduced in (Betteti and Laurenti, 2026): ℋθH(z)=12‖z‖22−∑h=1Lgh(z)⊤bh−∑h=2Lℱh(gh(z)),H_ _H(z)= 12\|z\|_2^2- _h=1^Lg_h(z) b_h- _h=2^LF_h\! (g_h(z) ), (6) where L≥2L≥ 2 is the number of layers, with widths N1=d,N2,…,NLN_1=d,N_2,…,N_L. We set g1(z)=zg_1(z)=z, Ψ1(z)=z _1(z)=z and recursively define the hidden feature vectors as gh(z)=Wh(h−1)Ψh−1(gh−1(z)),h=2,…,L.g_h(z)=W_h(h-1) _h-1\! (g_h-1(z) ), h=2,…,L. (7) Here, ℱh:Nh→F_h: ^N_h→ are convex differentiable potentials and Ψh=∇ℱh _h= _h their continuous gradients, Wh(h+1)∈Nh×Nh+1W_h(h+1)∈ ^N_h× N_h+1 are interaction matrices between consecutive layers, with tied weights W(h+1)h=Wh(h+1)⊤W_(h+1)h=W_h(h+1) , and bh∈Nhb_h∈ ^N_h are bias vectors entering the Hamiltonian directly. The parameters θH _H collect the interaction matrices, biases, and any trainable parameters of the potentials. Notice that each gh(z)g_h(z) is obtained by composing all preceding layers, allowing increasingly complex features to be encoded throughout the architecture. In addition, each hidden layer gHg_H is a standard MLP, to which the universal approximation theorem applies (HK-SM-WH:89): consequently, the model can in principle approximate any vector field given sufficient data and parameters. The following lemma guarantees that a bounded first hidden activation is sufficient to preserve the quadratic growth of the Hamiltonian at infinity. Lemma 2. Assume that the first hidden activation Ψ2 _2 has bounded image. Then, ℋθHH_ _H is coercive. That is, there exist constants c0,c1≥0c_0,c_1≥ 0 such that ℋθH(z)≥12‖z‖22−c1‖z‖2−c0.H_ _H(z)≥ 12\|z\|_2^2-c_1\|z\|_2-c_0. Lemma 2 follows by the boundedness of Ψ2=∇ℱ2, _2= _2, which implies that ℱ2F_2 grows at most linearly, while all feature vectors gh(z)g_h(z) with h≥3h≥ 3 remain bounded. Hence, the negative terms in (6) grow at most linearly in ‖z‖2\|z\|_2 and are dominated by its quadratic leading term. Coercivity and continuity of ℋθHH_ _H guarantee that its sublevel sets are compact despite the potentially nonconvex hidden-layer contributions. This is a key property we rely on in the following subsection to construct ρ -Robust Invariant Sets. 3.1 ρ -Robust Invariant Set Computation To characterize ρ -Robust Invariant Sets for (4), we rely on barrier certificates (Ames et al., 2017). In particular, we show that an isolated local minimum of ℋθHH_ _H induces a compact energy sublevel component whose robustness to external inputs can be explicitly quantified. To show that, we let z⋆z_ be an isolated local minimizer of ℋθHH_ _H and define the shifted energy V⋆(z):=ℋθH(z)−ℋθH(z⋆)V_ (z):=H_ _H(z)-H_ _H(z_ ). For ϵ>0ε>0, consider ϵ,⋆:=Compz⋆z:hϵ,⋆(z)≥0,hϵ,⋆(z):=ϵ−V⋆(z),Γϵ,⋆:=∂ϵ,⋆,C_ε, :=Comp_z_ \z:h_ε, (z)≥ 0 \, h_ε, (z):=ε-V_ (z), _ε, := _ε, , (8) where Compz⋆Comp_z_ denotes the connected component containing z⋆z_ . Since ℋθHH_ _H is coercive, ϵ,⋆C_ε, is compact. We assume that ℋθH(z⋆)=ϵH_ _H(z_ )=ε is a regular value of ℋθHH_ _H, so that ∇ℋθH(z)≠0 _ _H(z)≠ 0 for every z∈Γϵ,⋆z∈ _ε, . Given a norm ∥⋅∥\|·\| on the input space, we denote the corresponding dual norm by ∥⋅∥∗\|·\|_* and define dH(z):=∇ℋθH(z)⊤RθR(z)∇ℋθH(z)≥0,aH(z):=GθG(z)⊤∇ℋθH(z)∈k.d_H(z):= _ _H(z) R_ _R(z) _ _H(z)≥ 0, a_H(z):=G_ _G(z) _ _H(z)∈ ^k. Then, in Theorem 3 below, we rely on the fact that along the dynamics (4) it holds that h˙ϵ,⋆(z,u)=dH(z)−aH(z)⊤u≥dH(z)−‖aH(z)‖∗‖u‖, h_ε, (z,u)=d_H(z)-a_H(z) u≥ d_H(z)-\|a_H(z)\|_*\|u\|, (9) to prove that hϵ,⋆(z,u)h_ε, (z,u) is a valid barrier function for System 4. Theorem 3 (ρ -robust invariant energy sets). Let ℋθHH_ _H be defined in (6), and consider the port-Hamiltonian dynamics in (4). For every z∈Γϵ,⋆z∈ _ε, , define ρBF(z):=dH(z)‖aH(z)‖∗,aH(z)≠0,+∞,aH(z)=0. _BF(z):= cases d_H(z)\|a_H(z)\|_*,&a_H(z)≠ 0,\\[5.16663pt] +∞,&a_H(z)=0. cases (10) Then the vector field is inward-pointing or tangent to ϵ,⋆C_ε, at z for every input satisfying ‖u‖≤ρBF(z)\|u\|≤ _BF(z). Consequently, ϵ,⋆C_ε, in (8) is an ρϵ,⋆ _ε, -robust invariant set, where ρϵ,⋆:=infz∈Γϵ,⋆ρBF(z). _ε, := _z∈ _ε, _BF(z). (11) Proof. For every z∈Γϵ,⋆z∈ _ε, , the barrier function condition implies h˙ϵ,⋆(z,u)≥0 h_ε, (z,u)≥ 0 whenever ‖u‖≤ρBF(z)\|u\|≤ _BF(z). The uniform radius in (11) therefore makes the vector field inward-pointing or tangent at every boundary point. Invariance then follows from the standard barrier-certificate conditions introduced in Appendix B.3. ∎ Theorem 3 provides an explicit robustness radius quantifying the tolerance of each energy component ϵ,⋆C_ε, to external inputs. Remark 1. The construction extends directly to bounded additive perturbations. In particular, consider z˙=fθ(z,u)+w z=f_ (z,u)+w, w∈(z)w (z) where (z)W(z) is nonempty and compact, and define its support function as σ(z)(q):=supw∈(z)q⊤w. _W(z)(q):= _w (z)q w. In this case, h˙ϵ,⋆≥dH(z)−‖aH(z)‖∗‖u‖−σ(z)(∇ℋθH(z)) h_ε, ≥ d_H(z)-\|a_H(z)\|_*\|u\|- _W(z) ( _ _H(z) ). Hence, provided that dH(z)≥σ(z)(∇ℋθH(z))d_H(z)≥ _W(z) ( _ _H(z) ) for all z∈Γϵ,⋆z∈ _ε, , the disturbance-robust radius is ρϵ,⋆:=infz∈Γϵ,⋆aH(z)≠0dH(z)−σ(z)(∇ℋθH(z))‖aH(z)‖∗, _ε, ^W:= _ subarraycz∈ _ε, \\ a_H(z)≠ 0 subarray d_H(z)- _W(z) ( _ _H(z) )\|a_H(z)\|_*, (12) with the quotient interpreted as +∞+∞ when aH(z)=0a_H(z)=0. Remark 2. The norm-ball certificate can also be generalized to arbitrary compact convex input sets. In particular, a compact convex set U is robustly admissible if σ(aH(z))+σ(z)(∇ℋθH(z))≤dH(z) _U(a_H(z))+ _W(z)( _ _H(z))≤ d_H(z), for all z∈Γϵ,⋆z∈ _ε, . Equivalently, the largest input envelope certified by this condition is ¯ϵ,⋆rob:=∩z∈Γϵ,⋆u∈k:aH(z)⊤u≤dH(z)−σ(z)(∇ℋθH(z)) U_ε, ^rob:=U _z∈ _ε, \u∈ ^k:a_H(z) u≤ d_H(z)- _W(z)( _ _H(z))\ (13) This set is convex because it is the intersection of U with a family of half-spaces. Thus, the construction accommodates norm balls, boxes, polytopes, ellipsoids, and asymmetric actuator constraints. Remark 3. A consequence of Theorem 15 is that if ℋθHH_ _H has finitely many isolated local minima zi⋆i=1M\z_i \_i=1^M, repeating the construction for suitable regular energy levels ϵi _i gives the robustness atlas :=(ϵi,i,¯ϵi,irob,ρϵi,i)i=1M, A:=\(C_ _i,i, U_ _i,i^rob, _ _i,i)\_i=1^M, (14) which assigns a certified operating envelope to each learned dynamical regime. 3.2 Computable robustness bounds from local geometry The exact robust radius in (11) requires minimizing a state-dependent ratio over the generally nonconvex boundary Γϵ,⋆ _ε, . Although exact, this characterization may be computationally demanding and does not explicitly reveal how the robustness margin depends on the local geometry of the learned model. We therefore derive a conservative lower bound in terms of the steepness of the energy shell, dissipation in its normal direction, and exposure of the normal direction to the input port. We then bound the shell-steepness term using a local Polyak-Łojasiewicz condition and derive an input-to-energy estimate under persistent forcing. For readability, in what follows, we suppress the parameter subscripts in ℋθHH_ _H, RθR_ _R, and GθG_ _G. We restrict attention to energy levels for which ℋ(z⋆)+ϵH(z_ )+ε is a regular value of ℋH, so that ∇ℋ(z)≠0 (z)≠ 0 on Γϵ,⋆ _ε, and the outward energy-normal direction η(z):=∇ℋ(z)‖∇ℋ(z)‖2 (z):= (z)\| (z)\|_2 is well defined. Since Γϵ,⋆ _ε, is compact, regularity also implies that the gradient norm is uniformly bounded away from zero on the boundary. We now define the following quantities rϵ,⋆:=infz∈Γϵ,⋆η(z)⊤R(z)η(z),g¯ϵ,⋆⟂:=supz∈Γϵ,⋆‖G(z)⊤η(z)‖∗,κϵ,⋆:=infz∈Γϵ,⋆‖∇ℋ(z)‖2.r_ε, := _z∈ _ε, (z) R(z) (z), g_ε, := _z∈ _ε, \|G(z) (z)\|_*, _ε, := _z∈ _ε, \| (z)\|_2. (15) Here, rϵ,⋆r_ε, measures normal dissipation, g¯ϵ,⋆⟂ g_ε, measures exposure of the energy shell to the input port, and κϵ,⋆ _ε, measures the minimum steepness of the energy shell. Regularity and compactness imply κϵ,⋆>0 _ε, >0, whereas rϵ,⋆>0r_ε, >0 is an additional dissipativity assumption: positive semidefiniteness of R alone allows the dissipation to vanish in the normal direction. We are now ready to state Theorem 4. Theorem 4 (Geometric lower bound on the robustness radius). Suppose that rϵ,⋆>0r_ε, >0 and g¯ϵ,⋆⟂≠0 g_ε, ≠ 0. Then ρϵ,⋆≥ρ¯ϵ,⋆geo:=rϵ,⋆κϵ,⋆g¯ϵ,⋆⟂. _ε, ≥ _ε, ^geo:= r_ε, _ε, g_ε, . (16) Consequently, ϵ,⋆C_ε, is a ρ¯ϵ,⋆geo _ε, ^geo-robust invariant set. The proof is provided in Appendix B.6. The theorem guarantees that a large overall input gain is not necessarily detrimental: only its component normal to the energy shell consumes the available robustness margin. We next lower-bound κϵ,⋆ _ε, using a local Polyak–Łojasiewicz (PL) condition (Karimi et al., 2016; Pengyun et al., 2023). In our framework, the Hamiltonian is a known scalar, while its gradient is implicitely defined: the local PL condition bounds the growth ∇ℋ on a shell Γϵ,⋆ _ε, with the specified energy budget ϵ>0ε>0, and therefore readily provides a local estimate of the ρϵ,⋆ _ε, -robustness radius. Since ℋθHH_ _H is C1C^1 with locally Lipschitz gradient, the appropriate local curvature object is its Clarke generalized Hessian, ∂C2ℋθH:=∂C(∇ℋθH) _C^2H_ _H:= _C( _ _H). Proposition 5 (Local Clarke-Hessian certificate for PL). Let ∇ℋθH(z⋆)=0 _ _H(z_ )=0, and suppose that ℋθHH_ _H is C1C^1 with locally Lipschitz gradient in a neighbourhood of z⋆z_ . Define m⋆:=minQ∈∂C2ℋθH(z⋆)λmin(Q)m_ := _Q∈ _C^2H_ _H(z_ ) _ (Q). If m⋆>0m_ >0, then, for every 0<μ⋆<m⋆0< _ <m_ , there exists a convex neighbourhood δN_ of z⋆z_ on which 12‖∇ℋθH(z)‖22≥μ⋆(ℋθH(z)−ℋθH(z⋆)). 12\| _ _H(z)\|_2^2≥ _ (H_ _H(z)-H_ _H(z_ ) ). (17) The proposition is strictly local and leaves the global Hamiltonian free to be nonconvex and multi-well. When ℋH is C2C^2 near z⋆z_ , it reduces to the classical positive-definite Hessian test. The proof and activation-specific formulas are provided in Appendix B.9. To translate the PL condition into a robustness bound, the selected component must satisfy ϵ,⋆⊆δC_ε, _ . We also require uniform dissipation along the energy gradient, i.e., the existence of r⋆>0r_ >0 such that dH(z)≥r⋆‖∇ℋ(z)‖22d_H(z)≥ r_ \| (z)\|_2^2 on δN_ , which is satisfied for neighbourhood of isolated local minima with no flat directions. Corollary 6 (PL-based robustness bound). Suppose (17) holds for some 0<μ⋆<m⋆0< _ <m_ . Then κϵ,⋆≥2μ⋆ϵ _ε, ≥ 2 _ ε and rϵ,⋆≥r⋆r_ε, ≥ r_ . Consequently, if g¯ϵ,⋆⟂>0 g_ε, >0, then ρϵ,⋆≥ρϵ,⋆PL:=r⋆2μ⋆ϵg¯ϵ,⋆⟂. _ε, ≥ _ε, ^PL:= r_ 2 _ ε g_ε, . (18) The same local PL induced geometry gives more than invariance of one prescribed shell, and extends the idea of energy-driven convergence to energy-tube confiment of sollicited trajectories. Proposition 7 (Input-to-energy tube). Let the uniform bounds ‖G(z)‖2≤g¯\|G(z)\|_2≤ g, and ‖u(t)‖2≤u¯\|u(t)\|_2≤ u hold in the component ϵ,⋆C_ε, and set c:=g¯u¯c:= g u. If c≤r⋆2μ⋆ϵc≤ r_ 2 _ ε then ϵ,⋆C_ε, is a u¯ u-Robust Invariant set. In particular, lim inft→∞hϵ,⋆(z(t))≥ϵ−c22r⋆2μ⋆. _t→∞h_ε, (z(t))≥ε- c^22r_ ^2 _ . (19) The condition c≤r⋆2μ⋆ϵc≤ r_ 2 _ ε guarantees that the tube centered at z⋆z_ with radius c/(2r⋆2μ⋆)c/(2r_ ^2 _ ) fits inside ϵ,⋆C_ε, , recovering for c=r⋆2μ⋆ϵc=r_ 2 _ ε the robust boundary margin from a dynamic energy estimate. Table 1: Prediction accuracy and certified robustness across nonlinear dynamical systems. Comparative performance of the pH-EBM across benchmarks addressing nonlinear system identification, nonconvex dynamics, and higher-dimensional controlled systems. Benchmark Published portHNN-u pH-EBM Metrics RMSE ↓ \> ρϵ,⋆ _ε, ↑ RMSE ↓ \> ρϵ,⋆ _ε, ↑ RMSE ↓ \> ρϵ,⋆ _ε, ↑ Silverbox [0.2930.293,-] [47.82147.821, 2e-22e-2] [0.4220.422, e-2e-1] CED [0.0540.054, -] [0.2170.217, 3e-23e-2] [0.0640.064, e-6e-1] Duffing double-well [-, -] [0.2540.254, 5e-45e-4] [0.0480.048, e-4e-1] 33-link pendulum [0.0510.051, -] [0.2760.276, 5e-35e-3] [0.0140.014, e-2e-1] NanoDrone 0.5s0.5\>s [13.71213.712, -] [-, -] [12.52112.521, e-6e-2] NanoDrone 5s5\>s [1363.7011363.701, -] [-, -] [369.025369.025, e-6e-2] 4 Experimental validation Benchmark suite and evaluation protocol. The experimental suite spans established nonlinear system identification (Silverbox (Wigren and Schoukens, 2013) and CED (T. and M., 2017))11 1 Dataset specifications and benchmark protocols are available at https://w.nonlinearbenchmark.org., nonconvex controlled mechanics (Duffing (Molero et al., 2012) and the n-link pendulum (Okamoto and Kojima, 2025)), and high-dimensional long-horizon flight dynamics (NanoDrone (Busetto et al., 2026)); the tasks and main results are summarized in Table 1. We compare against benchmark-specific literature references, portHNN-u, and dissipative or black-box NODE baselines where available. The reported pH-EBM models are trained with a Huber observation loss and a mixed AdamW-to-SGD optimization schedule. Prediction is evaluated using each benchmark’s native metric, while robustness is assessed using the common energy-shell and normalized-radius protocol detailed in Appendix C.7. The pH-EBM framework is available at https://github.com/sim1bet/energy-safe-dynamics. Accuracy and robustness across benchmarks. Across the suite, the pH-EBM combines competitive or improved predictive accuracy with substantially larger certified input margins (Table 1). On Silverbox and CED, it remains close to the strongest published errors (0.4220.422 vs. 0.2930.293 and 0.0640.064 vs. 0.0540.054, respectively) while reducing the portHNN-u error and increasing the corresponding robustness radius by approximately 10×10× and 20×20×. The advantage is strongest on the nonconvex mechanical systems: on Duffing, RMSE decreases from 0.2540.254 to 0.0480.048 while the radius increases from 5×10−45× 10^-4 to 4×10−14× 10^-1; on the 33-link pendulum, RMSE decreases from 0.2760.276 to 0.0140.014 while the radius increases from 5×10−35× 10^-3 to 2×10−12× 10^-1. The latter also improves on the published RMSE of 0.0510.051, and retains lower error under the Brownian and step out-of-distribution forcing tests reported in Appendix C.1, while requiring approximately 1/201/20 of the training time of the dissipative NODE reference. Finally, on NanoDrone the pH-EBM improves the reported RMSE at both 0.5s0.5\,s (12.52112.521 vs. 13.71213.712) and, more markedly, at 5s5\,s (369.025369.025 vs. 1363.7011363.701), while retaining a normalized certificate of 6×10−26× 10^-2. These results show that the structural constraints do not impose a systematic accuracy–robustness trade-off: across the tested systems, improved certification is compatible with competitive—and in the nonconvex and long-horizon regimes, substantially improved—prediction. Figure 2: Safe long-horizon reconstruction of NanoDrone flight dynamics with pH-EBM. a) Mean absolute error (MAE) of pH-EBM and the black-box baseline (Busetto et al., 2026) over a 5,s5,s Random rollout, ten times the 0.5,s0.5,s training horizon. The pH-EBM maintains substantially lower error, while the B error rapidly grows beyond the training window. b-c) Random-trajectory reconstruction relative to the projected certified region ϵ,⋆C_ε, at ϵ=−12.926ε=-12.926. The pH-EBM remains inside the certified set and turns inward near Γϵ,⋆ _ε, , whereas the B trajectory exits the region after approximately 2,s2,s. d) Extrapolation to the unseen Melon excitation. Although B is more accurate over the original 0.5,s0.5,s horizon, its long-rollout error quickly increases, while pH-EBM remains bounded, highlighting the distinction between short-horizon accuracy and safe recursive deployment. e-f) Two- and three-dimensional projections of the learned Hamiltonian and shell Γ−12.926,⋆ _-12.926, , revealing a strongly nonconvex yet radially unbounded energy landscape that supports certified long-horizon evolution. Robust long time-horizon rollouts in high-dimensional Nanodrone application. The NanoDrone benchmark considers identification of flight dynamics from collections of short, disjoint 0.5,s0.5,s trajectories acquired under markedly different excitation distributions (Random, Chirp, Square, and Melon). We evaluate the pH-EBM against the reference black-box neural model of Busetto et al. (2026) under two increasingly demanding settings: (i) long-horizon prediction over the training-supported S3 distributions (Random, Chirp, and Square), and (i) long-horizon prediction under the unseen Melon distribution, corresponding to the original extrapolation task. In both cases, models are rolled out for 5,s5,s, i.e., ten times beyond the duration of the trajectories available during training. This setting therefore tests not only whether the learned dynamics fit short flight segments, but whether they remain predictive and operationally well behaved when recursively propagated far beyond the training horizon. On the S3 distributions, the pH-EBM achieves substantially lower prediction error than the reference black-box model, including state-of-the-art short-horizon validation accuracy together with markedly improved long-horizon reconstruction. Figure 2.a shows that its mean absolute error (MAE) remains consistently smaller throughout the 5,s5,s rollout, demonstrating that the structured Hamiltonian dynamics retain the expressivity required for high-dimensional flight identification without sacrificing long-term stability. More importantly, this predictive advantage is accompanied by an explicit safety certificate. Figures 2.(b-c) show the learned trajectory evolving inside the certified energy component ϵ,⋆C_ε, : as it approaches the boundary Γϵ,⋆ _ε, , the pH-EBM dynamics remain inward-pointing and the trajectory does not leave the certified operating region. In contrast, the black-box NODE accurately tracks the Random trajectory over approximately the 0.5,s0.5,s training horizon, but subsequently drifts away and eventually exits the same projected safe region. Figures 2.(e-f) visualize the corresponding two- and three-dimensional projections of the learned Hamiltonian and its certified boundary shell, making the relation between the learned energy geometry and the long-horizon trajectory explicit. The unseen Melon distribution exposes a complementary and particularly relevant trade-off. Over the original short-horizon extrapolation task, the pH-EBM is less accurate than the reference black-box model. Yet this ordering reverses qualitatively as the rollout is extended: while the black-box prediction error grows rapidly and the trajectory develops an unstable mode, the pH-EBM remains confined to the certified region and its average MAE approaches a bounded long-horizon regime (Fig. 2.d). The experiment therefore separates two notions that short-horizon benchmark scores alone cannot distinguish: local predictive fitness and safe recursive deployment. A model may provide the smaller error immediately outside the training distribution while still generating unstable trajectories when propagated for longer horizons. The pH-EBM instead sacrifices some short-horizon Melon accuracy while retaining an a priori forward-invariance certificate and preventing the unstable long-term behavior observed in the unconstrained NODE. For safety-critical flight dynamics, where prediction is repeatedly propagated and instability can be more consequential than a moderate increase in instantaneous error, this distinction is central. Figure 3: Geometry reconstruction, accuracy, and certified robustness in the 2D Duffing oscillator. a) Ground-truth Duffing Hamiltonian (left) and energies reconstructed by portHNN-u (center) and pH-EBM (right). Both recover the characteristic nonconvex double-well geometry. The pH-EBM panel additionally reports the pointwise input tolerance, whose minimum occurs near the saddle separating the two wells. b) Theoretical (purple) and empirical (blue) state-independent radii ρϵ,⋆ _ε, across an energy-level sweep. The certificate decreases as the shell approaches the saddle energy, then increases again after the certified component expands beyond it, revealing reduced input tolerance near transitions between operating regimes. c) Reconstruction error and empirical robustness radius for portHNN-u and pH-EBM. While both models achieve accurate trajectory reconstruction, the pH-EBM yields certified radii approximately 22-33 orders of magnitude larger than portHNN-u. Nonconvex geometry and attractor-dependent certification. The Duffing oscillator (Molero et al., 2012) provides a canonical controlled system with a nonconvex double-well potential, where external forcing drives transitions between distinct operating regimes. We compare the pH-EBM against the port-Hamiltonian portHNN-u neural architecture of Desai et al. (2021), evaluating both predictive reconstruction and the geometry of the resulting ρ -safety certificate. In the non-chaotic regime, both models recover the characteristic double-well energy landscape from input-output trajectories (Fig. 3.a), showing that accurate identification can preserve the underlying nonconvex structure. The pH-EBM further reveals how certification depends on the learnt geometry and NODE specifications. The pointwise radius ρBF(z) _BF(z) decreases near the saddle separating the two wells, and the state-independent radius ρϵ,⋆ _ε, correspondingly weakens as the certified shell approaches the critical energy before increasing again beyond it (Fig. 3.b). Saddle-proximity weakenanig of the certificate directly illustrates the theoretical link between critical energy levels and reduced input tolerance. Importantly, similar predictive reconstructions do not imply comparable safety margins: under the adopted configuration, the pH-EBM attains a state-independent certificate approximately two orders of magnitude larger than portHNN-u (Fig. 3.c). Robustness is therefore governed not by prediction error alone, but by the energy geometry and dissipation structure induced by the learned architecture. Out-of-distribution robustness and mechanical scaling. The n-link pendulum provides a challenging nonlinear mechanical benchmark whose complexity increases rapidly with the number of coupled joints. For n=2n=2 and n=3n=3, we compare the pH-EBM against the naive and dissipative neural ODEs studied in Okamoto and Kojima (2025) and the portHNN-u model. Across these systems, the pH-EBM improves on the reference models in three key respects: predictive accuracy, robustness to adversarial inputs, and training efficiency (Appendix C.1). In particular, under both Brownian and step out-of-distribution forcing, the pH-EBM combines lower prediction error with a certified ρ -robust invariant region. Moreover, training requires approximately 1/201/20 of the time for the dissipative reference, illustrating that stability and safety need not come at the computational cost of externally imposed dissipativity mechanisms. 5 Conclusions Limitations of the pH-EBM framework. The proposed port-Hamiltonian energy-based model combines expressive nonlinear dynamics with explicit bounded-rollout and safety certificates, but introduces two practical costs: training complexity and hyperparameter sensitivity. In our experiments, pH-EBMs generally require longer training and more targeted model selection than unconstrained NODEs, although their computational cost remains comparable to portHNN-u and can be substantially lower than that of alternative stable architectures such as dissipative NODEs. A more fundamental limitation concerns the scope of the safety guarantee: the certificates derived here apply directly to the learned model dynamics, not automatically to the unknown physical plant. Extending invariance guarantees to the true system requires a validated enclosure of model mismatch and external disturbances, which must then be incorporated explicitly into the robust barrier condition. Reducing this model-to-plant gap, together with the training and model-selection burden, is therefore central to future deployment in safety-critical settings. Future directions. The pH-EBM provides a general substrate for learning nonlinear dynamics while retaining stability and quantitative safety guarantees. A natural next step is to extend the framework to physical interactions not explicitly modeled here, including contact and hybrid dynamics, and to integrate the learned model with feedback controllers rather than treating external inputs as prescribed signals. Developing scalable certificate-aware training and closed-loop safety guarantees would move the framework from certified system identification toward end-to-end deployment in safety-critical applications such as autonomous flight, robotic manipulation, and human–robot interaction. AI usage disclosure. Large language models, including ChatGPT and Claude, were used to assist with language polishing, codebase organization and optimization, and the implementation of dataset interfaces for the experimental framework. They were not used for research ideation, formulation of the scientific contributions, development of the theoretical results, or derivation and verification of mathematical proofs. All scientific content, experimental design, implementation choices, and reported results were reviewed and validated by the authors. Reproducibility statement. We provide extensive material to facilitate reproduction of both the theoretical and experimental results. The Appendix states the assumptions underlying the theoretical guarantees and contains complete proofs of the main results. It also documents the datasets, preprocessing, model architectures, training protocols, hyperparameters, evaluation procedures, and certificate computations used throughout the experiments. The accompanying supplementary material includes the implementation of the pH-EBM framework together with the experiment configurations and scripts required to reproduce the reported benchmarks and robustness analyses. References Abate et al. (2021) A. Abate, D. Ahmed, A. Edwards, M. Giacobbe, and A. Peruffo FOSSIL: a software tool for the formal synthesis of lyapunov functions and barrier certificates using neural networks. In Proceedings of the 24th International Conference on Hybrid Systems: Computation and Control, HSCC ’21, p. 1–11. External Links: Document Cited by: §2. Ames et al. (2019) A. D. Ames, S. Coogan, M. Egerstedt, G. Notomista, K. Sreenath, and P. Tabuada Control barrier functions: theory and applications. In 2019 18th European Control Conference (ECC), p. 3420–3431. External Links: Document Cited by: §B.3. Ames et al. (2017) A. D. Ames, X. Xu, J. W. Grizzle, and P. Tabuada Control barrier function based quadratic programs for safety critical systems. IEEE Transactions on Automatic Control 62 (8), p. 3861–3876. External Links: ISSN 1558-2523, Document Cited by: §B.3, §3.1. Ashmore et al. (2021) R. Ashmore, R. Calinescu, and C. Paterson Assuring the machine learning lifecycle: desiderata, methods, and challenges. ACM Computing Surveys 54 (5), p. 1–39. External Links: ISSN 1557-7341, Document Cited by: §1. Betteti and Laurenti (2026) S. Betteti and L. Laurenti Hybrid energy-based models for physical ai: provably stable identification of port-hamiltonian dynamics. arXiv. External Links: Document Cited by: §B.2, §3. Blanchini and Miani (2015) F. Blanchini and S. Miani Set-theoretic methods in control. Springer International Publishing. External Links: ISBN 9783319179339, ISSN 2324-9757, Document Cited by: §2. Blanchini (1999) F. Blanchini Set invariance in control. Automatica 35 (11), p. 1747–1767. External Links: ISSN 0005-1098, Document Cited by: §2. Brunke et al. (2022) L. Brunke, M. Greeff, A. W. Hall, Z. Yuan, S. Zhou, J. Panerati, and A. P. Schoellig Safe learning in robotics: from learning-based control to safe reinforcement learning. Annual Review of Control, Robotics, and Autonomous Systems 5 (1), p. 411–444. External Links: ISSN 2573-5144, Document Cited by: §1, §1. Busetto et al. (2026) R. Busetto, E. Cereda, M. Forgione, G. Maroni, D. Piga, and D. Palossi Nonlinear system identification for a nano-drone benchmark. Control Engineering Practice 172, p. 106871. External Links: ISSN 0967-0661, Document Cited by: Figure 2, §4, §4. Chen et al. (2018) R. T. Q. Chen, Y. Rubanova, J. Bettencourt, and D. K. Duvenaud Neural ordinary differential equations. Vol. 31. Cited by: §1. Clarke (1975) F. H. Clarke Generalized gradients and applications. Transactions of the American Mathematical Society 205, p. 247–247. External Links: ISSN 0002-9947, Document Cited by: §A.1. Dawson et al. (2023) C. Dawson, S. Gao, and C. Fan Safe control with learned certificates: a survey of neural lyapunov, barrier, and contraction methods for robotics and control. IEEE Transactions on Robotics 39 (3), p. 1749–1767. External Links: ISSN 1941-0468, Document Cited by: §2. Desai et al. (2021) S. A. Desai, M. Mattheakis, D. Sondak, P. Protopapas, and S. J. Roberts Port-hamiltonian neural networks for learning explicit time-dependent dynamical systems. Physical Review E 104 (3). External Links: ISSN 2470-0053, Document Cited by: §4. Greydanus et al. (2019) S. Greydanus, M. Dzamba, and J. Yosinski Hamiltonian neural networks. Vol. 32. Cited by: §1. Gyevnár and Kasirzadeh (2025) B. Gyevnár and A. Kasirzadeh AI safety for everyone. Nature Machine Intelligence 7 (4), p. 531–542. External Links: ISSN 2522-5839, Document Cited by: §1. Hewing et al. (2020) L. Hewing, K. P. Wabersich, M. Menner, and M. N. Zeilinger Learning-based model predictive control: toward safe learning in control. Annual Review of Control, Robotics, and Autonomous Systems 3 (1), p. 269–296. External Links: ISSN 2573-5144, Document Cited by: §1. Hoover et al. (2022) B. Hoover, D. H. Chau, H. Strobelt, and D. Krotov A universal abstraction for hierarchical hopfield networks. In The Symbiosis of Deep Learning and Differential Equations I, Cited by: §B.2, §1. Hoover et al. (2023) B. Hoover, Y. Liang, B. Pham, R. Panda, H. Strobelt, D. H. Chau, M. Zaki, and D. Krotov Energy transformer. Advances in neural information processing systems 36, p. 27532–27559. Cited by: §1. Kang et al. (2021) Q. Kang, Y. Song, Q. Ding, and W. P. Tay Stable neural ode with lyapunov-stable equilibrium points for defending against adversarial attacks. In Advances in Neural Information Processing Systems, Vol. 34, p. 14925–14937. Cited by: §1. Karimi et al. (2016) H. Karimi, J. Nutini, and M. Schmidt Linear convergence of gradient and proximal-gradient methods under the polyak-Łojasiewicz condition. In Machine Learning and Knowledge Discovery in Databases, p. 795–811. External Links: ISBN 9783319461281, ISSN 1611-3349, Document Cited by: §B.7, §3.2. Karniadakis et al. (2021) G. E. Karniadakis, I. G. Kevrekidis, L. Lu, P. Perdikaris, S. Wang, and L. Yang Physics-informed machine learning. Nature Reviews Physics 3 (6), p. 422–440. External Links: ISSN 2522-5820, Document Cited by: §1. Khalil (2002) H. K. Khalil Nonlinear systems. Prentice-Hall, Upper Saddle River, NJ. Note: The book can be consulted by contacting: PH-AID: Wallet, Lionel Cited by: §2. Kojima and Okamoto (2022) R. Kojima and Y. Okamoto Learning deep input-output stable dynamics. Vol. 35, p. 8187–8198. Cited by: §1. Kolter and Manek (2019) J. Z. Kolter and G. Manek Learning stable deep dynamics models. Vol. 32. Cited by: §1. Krotov and Hopfield (2020) D. Krotov and J. J. Hopfield Large associative memory problem in neurobiology and machine learning. In International Conference on Learning Representations, External Links: Document Cited by: §B.2, §1. Lawrence et al. (2020) N. Lawrence, P. Loewen, M. Forbes, J. Backstrom, and B. Gopaluni Almost surely stable deep dynamics. Vol. 33, p. 18942–18953. Cited by: §1. Lee (2012) J. M. Lee Introduction to smooth manifolds. Springer New York. External Links: ISBN 9781441999825, ISSN 0072-5285, Document Cited by: §A.1. Legaard et al. (2023) 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 systems. ACM Computing Surveys 55 (11), p. 1–34. External Links: ISSN 1557-7341, Document Cited by: §1. Massaroli et al. (2020) S. Massaroli, M. Poli, M. Bin, J. Park, A. Yamashita, and H. Asama Stable neural flows. arXiv. External Links: Document Cited by: §1. Molero et al. (2012) F. J. Molero, M. Lara, S. Ferrer, and F. Céspedes 2-d duffing oscillator: elliptic functions from a dynamical systems point of view. Qualitative Theory of Dynamical Systems 12 (1), p. 115–139. External Links: ISSN 1662-3592, Document Cited by: §4, §4. Okamoto and Kojima (2025) Y. Okamoto and R. Kojima Learning deep dissipative dynamics. Proceedings of the AAAI Conference on Artificial Intelligence 39 (18), p. 19749–19757. External Links: ISSN 2159-5399, Document Cited by: §4, §4. Pengyun et al. (2023) Y. Pengyun, F. Cong, and L. Zhouchen On the lower bound of minimizing polyak-Łojasiewicz functions. In Proceedings of Thirty Sixth Conference on Learning Theory, G. Neu and L. Rosasco (Eds.), Proceedings of Machine Learning Research, Vol. 195, p. 2948–2968. Cited by: §3.2. Roth et al. (2025) F. J. Roth, D. K. Klein, M. Kannapinn, J. Peters, and O. Weeger Stable port-hamiltonian neural networks. arXiv. External Links: Document Cited by: §1, §3. T. and M. (2017) W. T. and S. M. Coupled electric drives data set and reference models. Technical Report Uppsala Universitet, Uppsala University Sweden. Cited by: §4. Van der Schaft (2007) A. Van der Schaft Port-hamiltonian systems: an introductory survey. EMS Press. External Links: ISBN 9783985475384, Document Cited by: §B.1, §1, §3. Wigren and Schoukens (2013) T. Wigren and J. Schoukens Three free data sets for development and benchmarking in nonlinear system identification. In 2013 European Control Conference (ECC), p. 2933–2938. External Links: Document Cited by: §4. Willems (2007) J. C. Willems Dissipative dynamical systems. European Journal of Control 13 (2-3), p. 134–151. External Links: ISSN 0947-3580, Document Cited by: §3. Yang et al. (2022) A. Yang, J. Xiong, M. Raginsky, and E. Rosenbaum Input-to-state stable neural ordinary differential equations with applications to transient modeling of circuits. In Proceedings of The 4th Annual Learning for Dynamics and Control Conference, Proceedings of Machine Learning Research, Vol. 168, p. 663–675. Cited by: §1. Appendix A Preliminaries Notation. 1d 1_d denotes the d-dimensional vector with all ones, 0d 0_d the d-dimensional vector with all zeros, while ℐdI_d denotes the identity matrix in d×d. We denote compact subsets of d as sets ℬ⊂dB⊂ ⊂ ^d, and their boundary as ∂ℬ . For two real vectors x,yx,y of the same dimension, x⊤yx y denotes the standard inner product. Let f:d→f: ^d→ ; for f∈Ck(d)f∈ C^k( ^d), the function is k-times continuously differentiable. For f∈C2(d)f∈ C^2( ^d), the gradient of f is denoted as ∇f∇ f, and the Hessian as D(∇f)=D2fD(∇ f)=D^2f. The partial derivative of f with respect to the variable xix_i is denoted as ∂f/∂xi∂ f/∂x_i. Θ bounds the growth of f∼Θ(g)f (g) through the existence of constants a,b>0a,b>0 such and a function g such that ag(x)<f(x)<bg(x)a\ g(x)<f(x)<b\ g(x). Given functions g:n→dg: ^n→ ^d and f:d→f: ^d→ , we denote the composition of the two functions as f∘g(y)f g(y), for y∈ny∈ ^n. The abbreviation a.e. stands for almost everywhere, that is for all ℬ⊂dB⊂ ^d except zero measure sets. Given a matrix A∈d×dA∈ ^d× d, we denote with A⊤A its transpose. In case A is symmetric, λmin(A),λmax(A)∈ _ (A),\ _ (A)∈ denote its minimum and its maximum eigenvalue. We denote with Br(x)⊂dB_r(x)⊂ ^d the ball of radius r>0r>0 and center x∈dx∈ ^d. A.1 Lie and Clarke’s derivatives Definition 8 (Lie derivative). Let X:d→dX: ^d→ ^d be a locally Lipschitz continuous vector field and let ΦX:≥0→d _X: _≥ 0→ ^d denote the associated local flow. For a smooth function g:d→g: ^d→ , the Lie derivative of g along X is defined as ℒXg(x)=dt|t=0(g∘ΦXt)(x),L_Xg(x)= ddt|_t=0(g\ _X^t)(x), (20) and coincides with the directional derivative ℒXg(x)=∇Xg(x)=∇xg(x)⊤X(x).L_Xg(x)= _Xg(x)= _xg(x) X(x). (21) Lie derivatives [Lee, 2012, Chapters 3, 9] will be central for the characterization of the negative definiteness of the energy function along the trajectories generated by the EBM dynamics. We extend the notion of derivatives to nonsmooth functions by recalling Clarke’s generalized gradient and directional derivative [Clarke, 1975], which apply to any locally Lipschitz function. Definition 9 (Generalized gradient). Let f:d→f: ^d→ be a locally Lipschitz function. The generalized gradient of f at x∈dx∈ ^d is denoted as ∂f(x)∂ f(x) and is the convex hull of the set of limits limk→+∞∇f(x+hk), _k→+∞∇ f(x+h_k), (22) with hk→k→+∞0dh_k k→+∞ 0_d and f differentiable at x+hk∈dx+h_k∈ ^d for all k∈k . Definition 10 (Generalized directional derivative). Let f:d→f: ^d→ be a locally Lipschitz function and take v∈dv∈ ^d. Then the Clarke’s generalized directional derivative at x∈dx∈ ^d is defined as f∘(x,v) f (x;v) =limh→0dsupδ→0f(x+h+δv)−f(x+h)δ = _h→ 0_d _ → 0 f(x+h+ v)-f(x+h) (23) =maxξ∈∂f(x)ξ⊤v. = _ ∈∂ f(x) v. (24) In particular, combining (21) with (24), the Clarke Lie derivative of a locally Lipschitz function f along a vector v is given by ℒvf(x):=f∘(x,v)=maxξ∈∂f(x)ξ⊤v,L_vf(x):=f (x;v)= _ξ∈∂ f(x)ξ v, (25) extending the classical smooth definition to the nonsmooth setting. Appendix B Technical derivations B.1 Port-Hamiltonian energy balance A port-Hamiltonian model [Van der Schaft, 2007] is a dynamical system characterized by dynamics x˙ x =[J(x)−R(x)]∇ℋ+G(x)u+w =[J(x)-R(x)] +G(x)u+w (26) yp y_p =G(x)⊤∇ℋ =G(x) (27) where x∈dx∈ ^d is the state of the model, J(x)=−J(x)⊤∈d×dJ(x)=-J(x) ∈ ^d× d is the skew-symmetric interconnection matrix, 0⪯R(x)=R(x)⊤∈d×d0 R(x)=R(x) ∈ ^d× d is the symmetric positive-semidefinite dissipation matrix, G(x)∈d×kG(x)∈ ^d× k is the input port, u∈ku∈ ^k the external time varying input, w∈dw∈ ^d is an unknown but bounded disturbance, and yp∈ky_p∈ ^k the power-conjugate output port. Along (4), ℒfℋ _fH =∇ℋ⊤(J−R)∇ℋ+∇ℋ⊤(Gu+w) = (J-R) + (Gu+w) (28) =∇ℋ⊤J∇ℋ−∇ℋ⊤R∇ℋ+yp⊤u+∇ℋ⊤w. = J - R +y_p u+ w. (29) For every real vector q, skew-symmetry gives q⊤Jq=(q⊤Jq)⊤=q⊤J⊤q=−q⊤Jq,q Jq=(q Jq) =q J q=-q Jq, (30) so q⊤JθJq=0q J_ _Jq=0. Therefore, dissipativity reduces to ℒfℋ=−∇ℋ⊤R∇ℋ+yp⊤u+∇ℋ⊤w.L_fH=- R +y_p u+ w. (31) No convexity property of ℋH is required for dissipativity. Coercivity is used instead to make energy sublevel sets compact and thereby preclude escape to infinity along bounded-energy trajectories. Convergence to a particular minimum requires additional information on the invariant set where the dissipation vanishes and is not implied by coercivity alone. B.2 Derivation and coercivity of the hybrid Hamiltonian The original energy function for a L-layered modern Hopfield network has been introduced in [Krotov and Hopfield, 2020, Hoover et al., 2022]. We identify the visible state with g1=z∈dg_1=z∈ ^d and denote the hidden variables by gh∈Nhg_h∈ ^N_h, h=2,…,Lh=2,…,L. For each layer let ℱh:Nh→F_h: ^N_h→ be a proper convex differentiable potential and define the activation Ψh(gh)=∇ℱh(gh) _h(g_h)= _h(g_h). Adjacent layers are coupled by matrices Wh(h+1)∈Nh×Nh+1W_h(h+1)∈ ^N_h× N_h+1 with W(h+1)h=Wh(h+1)⊤W_(h+1)h=W_h(h+1) and each layer has bias term bh∈Nhb_h∈ ^N_h. The modern Hopfield energy is then ℰ(g1,…,gL)= (g_1,…,g_L)= −∑h=1L−1Ψh(gh)⊤Wh(h+1)Ψh+1(gh+1) - _h=1^L-1 _h(g_h) W_h(h+1) _h+1(g_h+1) +∑h=1L[gh⊤(Ψh(gh)−bh)−ℱh(gh)] + _h=1^L [g_h ( _h(g_h)-b_h )-F_h(g_h) ] (32) = = g1⊤Ψ1(g1)−ℱ1(g1)−g1⊤b1 \>\>g_1 _1(g_1)-F_1(g_1)-g_1 b_1 +∑h=2L[−Ψh−1(gh−1)⊤W(h−1)hΨh(gh)+gh⊤(Ψh(gh)−bh)−ℱh(gh)]. + _h=2^L [- _h-1(g_h-1) W_(h-1)h _h(g_h)+g_h ( _h(g_h)-b_h )-F_h(g_h) ]. (33) The recurrence of the dynamics, with all hidden layers taking contributions from both the layers above and below, makes training the full dynamical architecture a computational and temporal challenge. For this reason, in Betteti and Laurenti [2026] a hybrid modern Hopfield energy has been introduced, where only the visible layer maintains its dynamical characterization and all the hidden layers are expressed as static maps ϕh _h of the layer below. Specifically, by expressing x(h)=ϕh(x(h−1))x_(h)= _h(x_(h-1)) for h=2,…,Lh=2,…,L, both the energy and the autonomous energy-gradient contribution x˙∝−∇ℰ(x) x - (x) can be expressed in terms of the visible layer state. With g1=zg_1=z, the visible layer internal contribution to the Hamiltonian is z⊤Ψ1(z)−ℱ1(z) z _1(z)-F_1(z) =z⊤z−12‖z‖22 =z z- 12\|z\|_2^2 (34) =12‖z‖22. = 12\|z\|_2^2. (35) For every hidden layer h≥2h≥ 2, the feedforward constraint (7) yields gh⊤Ψh(gh) g_h _h(g_h) =Ψh−1(gh−1)⊤Wh(h−1)⊤Ψh(gh) = _h-1(g_h-1) W_h(h-1) _h(g_h) (36) =Ψh−1(gh−1)⊤W(h−1)hΨh(gh), = _h-1(g_h-1) W_(h-1)h _h(g_h), (37) which cancels exactly the coupling term between layers h−1h-1 and h in (33). Applying this cancellation for h=2,…,Lh=2,…,L leaves precisely ℰ(z,g2,…,gL) (z,g_2,…,g_L) ≡ℋθH(z) _ _H(z) =12‖z‖22−z⊤b1−∑h=2L(ℱh(gh(z))+gh(z)⊤bh) = 12\|z\|_2^2-z b_1- _h=2^L (F_h(g_h(z))+g_h(z) b_h ) (38) and the associated neural network action has autonomous contribution z˙∝−∇ℋθH(z) z - _ _H(z) (39) We next prove (43). Suppose the first hidden activation has bounded image: there exists M2<+∞M_2<+∞ such that ‖Ψ2(x)‖2≤M2\| _2(x)\|_2≤ M_2 for all x∈N2x∈ ^N_2. Since Ψ2=∇ℱ2 _2= _2, the mean-value inequality gives |ℱ2(z)−ℱ2(0)|≤M2‖z‖2,|F_2(z)-F_2(0)|≤ M_2\|z\|_2, (40) so ℱ2(W21z)≤|ℱ2(0)|+M2‖W21‖2‖z‖2.F_2(W_21z)≤|F_2(0)|+M_2\|W_21\|_2\|z\|_2. (41) Moreover, g3=W32Ψ2(g2)g_3=W_32 _2(g_2) belongs to a compact set independent of z. By continuity of Ψ3 _3, its image over this compact set is bounded; recursively, every gh(z)g_h(z) for h≥3h≥ 3 lies in a compact set. Continuity of ℱhF_h and compactness of the layer variables ghg_h then yields constants Ch<+∞C_h<+∞ such that ℱh(gh(z))+gh(z)⊤bh≤Ch,h=3,…,L.F_h(g_h(z))+g_h(z) b_h≤ C_h, h=3,…,L. (42) Substituting these estimates into (6) yields ℋθH(z)≥12‖z‖22−c1‖z‖2−c0H_ _H(z)≥ 12\|z\|_2^2-c_1\|z\|_2-c_0 (43) with c1=‖W21‖2(M2+b2)+b1,c0=|ℱ2(0)|+∑h=3LCh.c_1=\|W_21\|_2(M_2+b_2)+b_1, c_0=|F_2(0)|+ _h=3^LC_h. (44) The quadratic term dominates the affine lower bound as ‖z‖2→+∞\|z\|_2→+∞, proving coercivity. B.3 Barrier functions and invariance Consider a controlled system z˙=f(z,u,w) z=f(z,u,w) and a continuously differentiable function h:ℝd→ℝh:R^d . The superlevel set :=z:h(z)≥0C:=\z:h(z)≥ 0\ (45) is invariant if every solution initialized in C remains in C for all future times for which the solution exists. A standard way to certify this property is through a zeroing barrier function. Let α be an extended class-K function, i.e., a continuous strictly increasing function satisfying α(0)=0α(0)=0. For a fixed admissible input/disturbance pair, the inequality ℒfh(z,u,w)+α(h(z))≥0L_fh(z,u,w)+α(h(z))≥ 0 (46) prevents h from decreasing through zero. Under the usual solution-existence assumptions, if (46) holds on a neighborhood of C, then C is invariant. If it holds for every u in a prescribed input set and every w in a prescribed disturbance set, the same set is robustly invariant for those exogenous signals [Ames et al., 2017, Ames et al., 2019]. On a regular boundary point z∈∂z∈ , where h(z)=0h(z)=0 and ∇h(z)≠0∇ h(z)≠ 0, the condition reduces to the geometric inward-pointing requirement ℒfh(z,u,w)≥0.L_fh(z,u,w)≥ 0. (47) Equivalently, the vector field must lie in the tangent cone of C; this is the local form of Nagumo’s invariance condition. B.4 BF invariance, robust input sets, and boundary geometry We first distinguish three levels of certification that are used in the main text: a pointwise zeroing-BF condition, a boundary-only invariance condition, and a state-independent set of exogenous inputs that satisfies the boundary condition everywhere. In the following treatment, we will refer to a specific condition being robustly satisfied when it holds over the entire support W of the unknown disturbances. Let hϵ(z)=ϵ−ℋ(z)h_ε(z)=ε-H(z) and consider the perturbed port-Hamiltonian dynamics z˙=[J(z)−R(z)]∇ℋ(z)+G(z)u+w,w∈(z), z=[J(z)-R(z)] (z)+G(z)u+w, w (z), (48) where (z)⊂dW(z)⊂ ^d is nonempty and compact. Write dH:=∇ℋ⊤R∇ℋd_H:= R , and aH:=G⊤∇ℋa_H:=G . Then ℒfhϵ=dH−aH⊤u−∇ℋ⊤w.L_fh_ε=d_H-a_H u- w. (49) Therefore the zeroing-BF inequality ℒfhϵ+α(hϵ)≥0L_fh_ε+ (h_ε)≥ 0 for every w∈(z)w (z) is equivalent to aH(z)⊤u+σ(z)(∇ℋ(z))≤dH(z)+α(hϵ(z)).a_H(z) u+ _W(z)( (z))≤ d_H(z)+ (h_ε(z)). (50) Setting (z)=0W(z)=\0\ recovers the result without bounded perturbayions. If aH(z)=0a_H(z)=0, no division by ‖aH(z)‖∗\|a_H(z)\|_* is performed: the energy effect of the input is identically zero to first order and feasibility depends only on dissipation, the barrier slack, and the disturbance term. For any norm ∥⋅∥\|·\| on input space with dual norm ∥⋅∥∗\|·\|_*, Hölder’s inequality gives aH⊤u≤‖aH‖∗‖u‖.a_H u≤\|a_H\|_*\|u\|. (51) Thus (10) is sufficient in the nominal case. It is the largest origin-centered norm ball contained in the unconstrained half-space because sup‖u‖≤ρaH⊤u=ρ‖aH‖∗. _\|u\|≤ a_H u= \|a_H\|_*. (52) B.5 Exact robust boundary certificate Let z⋆z_ be an isolated local minimum and define ϵ,⋆:=Compz⋆z:hϵ,⋆(z)≥0,hϵ,⋆(z):=ϵ−V⋆(z),Γϵ,⋆:=∂ϵ,⋆.C_ε, :=Comp_z_ \z:h_ε, (z)≥ 0\, h_ε, (z):=ε-V_ (z), _ε, := _ε, . (53) The energy component is both connected (by definition) and compact by coercivity of ℋH, and the solutions of (48) exist uniquely for the admissible measurable inputs under consideration. On Γϵ,⋆ _ε, one has by definition that hϵ,⋆=0h_ε, =0, therefore the robust inward-pointing condition is aH(z)⊤u+σ(z)(∇ℋ(z))≤dH(z).a_H(z) u+ _W(z)( (z))≤ d_H(z). (54) Before proceeding with the first theorem, we provide a general notion from differential geometry, that is, that of regular hypersurface, which will be used in reference to the boundaries of our candidate safe sets ϵ,⋆C_ε, . Proposition 11 (Regular (d−1)(d-1)-dim hypersurface). Let f:d→f: ^d→ be a C1C^1 function, and for s∈s∈ define the level set Γs:=z∈d:f(z)=s _s:=\z∈ ^d:\ f(z)=s\ (55) Then we say that Γs _s is a regular (d−1)(d-1)-dim C1C^1-hypersurface if ∇f(z)≠0∀z∈Γs.∇ f(z)≠ 0 ∀ z∈ _s. (56) When the Hamiltonian ℋH has isolated critical points, the shells Γϵ,⋆ _ε, are regular almost everywhere ϵ>0ε>0, as there will be finitely many values of ϵε such that a saddle point or local maxima lie in Γϵ,⋆ _ε, . The existence of entire connected open subsets ⊂dV⊂ ^d where the regularity assumption is not satisfied, and hence such that ℋ(z)=cH(z)=c for all z∈z , is possible but unlikely. The existence of such regions would require the following two conditions to hold simultaneously ∇ℋ(z)=0,det∇2ℋ(z)=0∀z∈. (z)=0, ∇^2H(z)=0 ∀ z . (57) More specifically, if the Hamiltonian ℋH is constituted only by primitives ℱhh=1L\F_h\_h=1^L of real analytic activation functions - such as the softmaxsoftmax, tanhtanh, logisticlogistic, and softplussoftplus - then the existence of open subsets V with constant energy is impossible due to analytic continuation. If such V existed, it would imply that ℋ≡cH≡ c on all d, which contraddicts the coercivity of ℋH. Therefore, the condition in (57) can only be satisfied by neural networks characterized by less regular activation functions with broad plateaus - such as ReLUReLU - and even in such cases remains unlikely. In the following theorem, we specialize the definition of safe sets as defined by the energy BF to our port-Hamiltonian framework. Theorem 12 (Maximal uniform robust input set). Let Γϵ,⋆ _ε, be the regular boundary of the zero level set hϵ,⋆=0h_ε, =0 in the sense of Proposition 11, and define ¯ϵ,⋆rob:=∩⋂z∈Γϵ,⋆u:aH(z)⊤u+σ(z)(∇ℋ(z))≤dH(z). U_ε, ^rob:=U∩ _z∈ _ε, \u:a_H(z) u+ _W(z)( (z))≤ d_H(z) \. (58) Then (i) for any locally essentially bounded input satisfying u(t)∈¯ϵ,⋆robu(t)∈ U_ε, ^rob almost everywhere, ϵ,⋆C_ε, is robustly invariant for every measurable disturbance w(t)∈(z(t))w(t) (z(t)); (i) ¯ϵ,⋆rob U_ε, ^rob is maximal with respect to set inclusion among state-independent subsets of U whose every element satisfies (54) at every boundary state. Proof. (i) let z∈Γϵ,⋆z∈ _ε, , and observe that aH(z)⊤u+σ(z)(∇ℋ(z))≤dH(z)⇔ℒfhϵ,⋆(z,u,w)≥0a_H(z) u+ _W(z)( (z))≤ d_H(z) _fh_ε, (z,u,w)≥ 0 for every w∈(z)w (z). Intersecting these pointwise admissible half-spaces over the full boundary yields exactly ¯ϵ,⋆rob U_ε, ^rob. ℒfh≥0L_fh≥ 0 on Γϵ,⋆ _ε, implies that the neural network action is in the tangent cone of the sublevel component, so Nagumo’s tangent condition yields robust invariance. (i) Let ϵ,⋆−U_ε, ^- be such that ϵ,⋆−/¯ϵ,⋆rob≠∅U_ε, ^-/ U_ε, ^rob≠ . Then there exists u∈ϵ,⋆−u _ε, ^- such that aH(z)⊤u+σ(z)(∇ℋ(z))>dH(z)a_H(z) u+ _W(z)( (z))>d_H(z) (59) and ϵ,⋆−U_ε, ^- does not make ϵ,⋆C_ε, robustly invariant. ∎ It is now easy to generalize the previous definition of safety ensured input set to arbitrary geometries instead of normed balls. Corollary 13 (Certification of arbitrary compact input sets). Let 0⊆U_0 be compact. Then 0⊆¯ϵ,⋆robU_0 U_ε, ^rob if and only if σ0(aH(z))+σ(z)(∇ℋ(z))≤dH(z)∀z∈Γϵ,⋆ _U_0(a_H(z))+ _W(z)( (z))≤ d_H(z) ∀ z∈ _ε, (60) Proof. Let ζ(z,u)=aH⊤u+σ(z)(∇ℋ(z))−dH(z) (z,u)=a_H u+ _W(z)( (z))-d_H(z). The function ζ is continuous in both of its arguments. Since both Γϵ,⋆ _ε, and 0U_0 are compact, then by Weierstrass theorem ζ constrained to this sets admits maximum and minimum with respect to both of its arguments. The certification of arbitrary compact input sets then requires 0 0 <minz∈Γϵ,⋆maxu∈0ζ(z,u) < _z∈ _ε, _u _0 (z,u) =minz∈Γϵ,⋆[maxu∈0aH⊤u]+σ(z)(∇ℋ(z))−dH(z) = _z∈ _ε, [ _u _0a_H u ]+ _W(z)( (z))-d_H(z) =minz∈Γϵ,⋆σ0(aH(z))+σ(z)(∇ℋ(z))−dH(z) = _z∈ _ε, _U_0(a_H(z))+ _W(z)( (z))-d_H(z) (61) ∎ Define the residual boundary dissipation after worst-case disturbance injection, δ(z):=dH(z)−σ(z)(∇ℋ(z)). (z):=d_H(z)- _W(z)( (z)). (62) Assume δ(z)≥0 (z)≥ 0 on the boundary, so that zero input is itself robustly admissible. The exact radius of the largest origin-centered norm ball contained in (58) is ρϵ,⋆:=infz∈Γϵ,⋆aH(z)≠0δ(z)‖aH(z)‖∗, _ε, := _ subarraycz∈ _ε, \\ a_H(z)≠ 0 subarray (z)\|a_H(z)\|_*, (63) Boundary points with aH(z)=0a_H(z)=0 impose no input-radius restriction as long as δ(z)≥0 (z)≥ 0. Consequently, every ρ≤ρϵ,⋆ ≤ _ε, satisfies u:‖u‖≤ρ∩⊆¯ϵ,⋆rob\u:\|u\|≤ \ U_ε, ^rob. A frequently useful but more conservative separable bound follows from minΓϵ,⋆δ(z)maxΓϵ,⋆‖aH(z)‖∗≤infz∈Γϵ,⋆aH(z)≠0δ(z)‖aH(z)‖∗. _ _ε, (z) _ _ε, \|a_H(z)\|_*≤ _ subarraycz∈ _ε, \\ a_H(z)≠ 0 subarray (z)\|a_H(z)\|_*. (64) This inequality corrects the tempting but generally false replacement of the infimum of a ratio by the ratio of independent extrema. B.6 Regular energy shells and normal robustness decomposition The regularity condition for the shell has deeper geometric implications, since it allows to clearly outline the vectors resulting in inward and outward motion with respect to the set ϵ,⋆C_ε, . Proposition 14 (Regular energy boundary). Assume that ∇ℋ(z)≠0,z∈Γϵ,⋆. (z)≠ 0, z∈ _ε, . (65) Then every point of Γϵ,⋆ _ε, belongs locally to the regular level set z:hϵ,⋆=0\z:h_ε, =0\, which is a C1C^1 embedded hypersurface of dimension (d−1)(d-1), and we can characterize the tangent space and the tangent cone at Γϵ,⋆ _ε, as TzΓϵ,⋆ T_z _ε, =v:∇ℋ(z)⊤v=0, =\v: (z) v=0\, (66) Tϵ,⋆(z) T_C_ε, (z) =v:∇ℋ(z)⊤v≤0. =\v: (z) v≤ 0\. (67) In addition, if the boundary is compact, then κϵ,⋆:=infΓϵ,⋆‖∇ℋ‖2>0 _ε, := _ _ε, \| \|_2>0. Proof. Consider the BF hϵ,⋆(z)h_ε, (z), and its directional derivative Dhϵ,⋆(z)[v]=−∇ℋ(z)⊤vDh_ε, (z)[v]=- (z) v. Since hϵ,⋆h_ε, maps d to , the derivative is surjective at all points z∈Γϵ,⋆z∈ _ε, such that ∇ℋ(z)≠0 (z)≠ 0. Since by assumption ∇ℋ(z)≠0 (z)≠ 0, we have that hϵ,⋆=0h_ε, =0 is a submersion and Γϵ,⋆ _ε, is a C1C^1 (d−1)(d-1)-hypersurface. Because ϵ,⋆C_ε, is the superlevel set hϵ,⋆≥0h_ε, ≥ 0, the feasible first-order directions are defined by Dhϵ,⋆(z)[v]≥0Dh_ε, (z)[v]≥ 0,or equivalently ∇ℋ(z)⊤v≤0 (z) v≤ 0, giving (67). Finally, by Weierstrass theorem the continuous function ‖∇ℋ‖2\| \|_2 attains a positive maximum and minimum on the compact boundary, which must be strictly positive given the assumption ∇ℋ(z)≠0 (z)≠ 0 on Γϵ,⋆ _ε, . ∎ The characterization of tangent space and tangent cone of Γϵ,⋆ _ε, , and in particular the nonzero-gradient condition is essential for a first-order barrier test. If ∇ℋ=0 =0 at a boundary point, ℒfhϵ,⋆L_fh_ε, vanishes and Dhϵ,⋆(z)[v]≡0Dh_ε, (z)[v]≡ 0 for any vector v∈dv∈ ^d: therefore, the BF condition fails to distinguish inward from outward motion. Thereby, regularity makes the energy gradient a valid boundary normal for the BF argument. For z∈Γϵ,⋆z∈ _ε, define the outward energy normal η(z):=∇ℋ(z)‖∇ℋ(z)‖2. (z):= (z)\| (z)\|_2. (68) Positive homogeneity of support functions and the identities dH(z) d_H(z) =‖∇ℋ(z)‖22η(z)⊤R(z)η(z), =\| (z)\|_2^2 (z) R(z) (z), (69) ‖aH(z)‖∗ \|a_H(z)\|_* =‖∇ℋ(z)‖2‖G(z)⊤η(z)‖∗, =\| (z)\|_2\|G(z) (z)\|_*, (70) σ(z)(∇ℋ(z)) _W(z)( (z)) =‖∇ℋ(z)‖2σ(z)(η(z)) =\| (z)\|_2 _W(z)( (z)) (71) show that, whenever G(z)⊤η(z)≠0G(z) (z)≠ 0, δ(z)‖aH(z)‖∗=‖∇ℋ(z)‖2η(z)⊤R(z)η(z)−σ(z)(η(z))‖G(z)⊤η(z)‖∗. (z)\|a_H(z)\|_*= \| (z)\|_2 (z) R(z) (z)- _W(z)( (z))\|G(z) (z)\|_*. (72) Interestingly, (72) reveals that, with respect to the input port, only the directions that are normal to the level set consumes the energy margin, while tangent directions leave it unaffected. In particular, if ‖G(z)⊤η(z)‖∗=0\|G(z) (z)\|_*=0 for all z∈Γϵ,⋆z∈ _ε, , the input is tangent to the Hamiltonian shells over the entire boundary and the energy certificate imposes no magnitude restriction on u; only the disturbance margin must remain nonnegative. Let rϵ,⋆ r_ε, :=infΓϵ,⋆η⊤Rη, := _ _ε, R , (73) g¯ϵ,⋆⟂ g_ε, :=supΓϵ,⋆‖G⊤η‖∗, := _ _ε, \|G \|_*, (74) w¯ϵ,⋆⟂ w_ε, :=supΓϵ,⋆σ(z)(η). := _ _ε, _W(z)( ). (75) If g¯ϵ,⋆⟂>0 g_ε, >0, and rϵ,⋆κϵ,⋆≥w¯ϵ,⋆⟂r_ε, _ε, ≥ w_ε, , then (72) yields ρϵ,⋆≥rϵ,⋆κϵ,⋆−w¯ϵ,⋆⟂g¯ϵ,⋆⟂, _ε, ≥ r_ε, _ε, - w_ε, g_ε, , (76) which proves (72). For each attractor component one may define the maximal regular-shell energy ϵ⋆reg:=supϵ¯>0:∇ℋ(z)≠0∀z∈Γe,⋆,∀ϵ∈(0,ϵ¯). _ ^reg:= \ ε>0: (z)≠ 0 ∀ z∈ _e, ,\>∀ε∈(0, ε) \. (77) Every compact shell below this threshold is eligible for the regular-boundary certificate. When a shell approaches a critical energy containing a saddle or another stationary point, κϵ,⋆ _ε, can collapse toward zero, quantitatively revealing the loss of first-order robustness. After the saddle point, the ρ -safety certificate exit resumes, but the new certified region will contain the sets on "both" sides of the saddle. B.7 Bounds under local PL and persistent forcing Assume the component ϵ,⋆C_ε, lies inside a neighborhood ⋆N_ of a local minimum z⋆∈dz_ ∈ ^d on which the local Polyak-Łojasiewicz (PL) condition [Karimi et al., 2016] holds for the Hamiltonian ℋH. 12‖∇ℋ(z)‖22≥μ⋆(ℋ(z)−ℋ(z⋆))=μ⋆V⋆(z),dH(z)≥r⋆‖∇ℋ(z)‖22, 12\| (z)\|_2^2≥ _ (H(z)-H(z_ ) )= _ V_ (z), d_H(z)≥ r_ \| (z)\|_2^2, (78) for constants μ⋆,r⋆>0 _ ,r_ >0 where r⋆≤minz∈ϵ,⋆λmin(R(z))r_ ≤ _z _ε, _ (R(z)). The PL condition is associated to a region of convex growth around a minima, and since we are requiring it to be satisfied for all points in ⋆N_ , then necessarily we will have that ⋆⊂ϵ⋆reg,⋆N_ _ _ ^reg, , where ϵ⋆reg _ ^reg has been defined in (77). Indeed, it holds that limϵ→ϵ⋆regκϵ,⋆=0 _ε→ _ ^reg _ε, =0 (79) and ⋆⊃ϵ⋆reg,⋆N_ _ _ ^reg, would lead, for a point z∈Γϵ⋆reg,⋆z∈ _ _ ^reg, to the contraddiction 0≥μ⋆V⋆(z)⏟>00≥ _ V_ (z)_>0 (80) Let R(x)≻0R(x) 0 over the entire component ϵ,⋆C_ε, . Then the PL condition implies dH(z)≥2r⋆μ⋆V⋆(z).d_H(z)≥ 2r_ _ V_ (z). (81) For the zeroing-BF condition inside the component, (81) show that aH(z)⊤u≤2r⋆μ⋆V⋆(z)+α(ϵ−V⋆(z))a_H(z) u≤ 2r_ _ V_ (z)+ (ε-V_ (z) ) (82) is a sufficient conservative approximation of the exact admissible half-space. If α(s)=γs (s)= s, the right-hand side is affine in q=V⋆(z)∈[0,ϵ]q=V_ (z)∈[0,ε]. Considering the first derivative of l(q)=2r⋆μ⋆q+γ(ϵ−q)l(q)=2r_ _ q+ (ε-q), we have l˙(q)=2r⋆μ⋆−γ l(q)=2r_ _ - . Consequently, if 2r⋆μ⋆>γ2r_ _ > , then l is strictly increasing over [0,ϵ][0,ε] and its minimum on the interval is ϵγε . Conversely, if γ>2r⋆μ⋆ >2r_ _ , the l is strictly decreasing over [0,ϵ][0,ε] and its minimum is 2r⋆μ⋆ϵ2r_ _ ε. Therefore, the minimum of the r.h.s. can be expressed as βϵ=ϵmin2r⋆μ⋆,γ. _ε=ε \2r_ _ , \. (83) Consequently, with Aϵ:=supz∈ϵ,⋆‖aH(z)‖∗<∞A_ε:= _z _ε, \|a_H(z)\|_*<∞, the bound ‖u‖≤βϵAϵ\|u\|≤ _εA_ε (84) ensures the chosen zeroing-BF inequality throughout the component whenever Aϵ>0A_ε>0. This is stronger than boundary-only invariance and is correspondingly more conservative. On Γϵ,⋆ _ε, , PL yields ‖∇ℋ(z)‖2≥2μ⋆ϵ>0.\| (z)\|_2≥ 2 _ ε>0. (85) Hence the boundary is regular by Proposition 14 and κϵ,⋆≥2μ⋆ϵ _ε, ≥ 2 _ ε. Substitution in the normal robustness bound gives the robust PL corollary ρϵ,⋆≥rϵ,⋆2μ⋆ϵ−w¯ϵ,⋆⟂g¯ϵ,⋆⟂, _ε, ≥ r_ε, 2 _ ε- w_ε, g_ε, , (86) whenever the numerator is nonnegative. This is an invariance-only certificate; it need not satisfy the selected zeroing-BF inequality at every interior point. B.8 Input-to-energy robustness under persistent forcing In this subsection we are going to present an additional bound on the confinement of the port-Hamiltonian trajectories z(t)z(t) whenever they are in a component ϵ,⋆C_ε, and we have uniform bounds on the input u, the input-port G and the disturbances w. Proposition 15 (Input-to-energy tube). Let the uniform bounds ‖G(z)‖2≤g¯,‖u(t)‖2≤u¯,‖w(t)‖2≤w¯.\|G(z)\|_2≤ g, \|u(t)\|_2≤ u, \|w(t)\|_2≤ w. (87) hold in the component ϵ,⋆C_ε, and set c:=g¯u¯+w¯c:= g u+ w. If c≤r⋆2μ⋆ϵ,c≤ r_ 2 _ ε, (88) then ϵ,⋆C_ε, is robustly invariant and every trajectory initialized in the component satisfies (). In particular, lim supt→∞V⋆(z(t))≤c22r⋆2μ⋆. _t→∞V_ (z(t))≤ c^22r_ ^2 _ . (89) Proof. Along the trajectories of the latent port-Hamiltonian EBM model, defined by dynamics (48) with vector field f(z)=[J(z)−R(z)]∇ℋ(z)+G(z)u+wf(z)=[J(z)-R(z)] (z)+G(z)u+w, we have the following Lie derivative on the difference function V⋆V_ . ℒfV⋆(z,u) _fV_ (z,u) =ℒfℋ(z,u) =L_fH\>(z,u) =−∇ℋ(z)R(z)∇ℋ(z)+(G(z)u+w)⊤∇ℋ(z) =- (z)R(z) (z)+ (G(z)u+w ) (z) ≤−r⋆‖∇ℋ(z)‖22+c‖∇ℋ(z)‖2 ≤-r_ \| (z)\|_2^2+c\| (z)\|_2 ≤−r⋆2‖∇ℋ(z)‖22+c22r⋆ ≤- r_ 2\| (z)\|_2^2+ c^22r_ ≤−r⋆μ⋆V⋆(z)+c22r⋆ ≤-r_ _ V_ (z)+ c^22r_ (90) where in the fourth passage we have used Young’s inequality ab≤εa2/2+b2/(ε2)ab≤ a^2/2+b^2/( 2) with a=‖∇ℋ‖a=\| \|, b=cb=c, and ε=r⋆ =r_ and in the last passage the PL condition. Under Assumption 88 we actually have on the level set hϵ,⋆=0h_ε, =0 −ℒfhϵ,⋆(z,u) -L_fh_ε, \>(z,u) =ℒfV⋆(z,u) =L_fV_ \>(z,u) ≤−r⋆μ⋆ϵ+(r⋆2μ⋆ϵ)22r⋆=0 ≤-r_ _ ε+ (r_ 2 _ ε)^22r_ =0 (91) and consequently we also have that ϵ,⋆C_ε, is robustly invariant for all (u,w)(u,w) satisfying Assumption 87. Finally, applying the comparison lemma on (90) we obtain V⋆(t)≤e−r⋆μ⋆tV⋆(0)+c22r⋆2μ⋆(1−e−r⋆μ⋆t).V_ (t)≤ e^-r_ _ tV_ (0)+ c^22r_ ^2 _ (1-e^-r_ _ t ). (92) Taking the limit t→+∞t→+∞ we conclude. ∎ Notice finally that when c<r⋆2μ⋆ϵc<r_ 2 _ ε, we are actually stating that asymptotically the trajectories of the system will be bounded to level sets hϵ,⋆=sh_ε, =s with s≥ϵ−c2/(2r⋆μ⋆)s≥ε-c^2/(2r_ _ ). B.9 Activation-dependent guidance for local PL We begin the current subsection by proving the general result of Proposition 5, and the proceed to specialize the result to an example architecture and different classes of activation functions. Proof. Since ℋ∈Cloc1,1(d)H∈ C^1,1_loc( ^d), then the Hessian ∇2ℋ(z)∇^2H(z) is defined a.e. z∈dz∈ ^d, and at the points of non-differentiability we relay on Clarke’s generalized Jacobian z↦∂(∇ℋ)(z)z ∂( )(z) as defined in Definition 9. The set ∂(∇ℋ)(z)∂( )(z) is nonempty, compact, convex, and upper semi-continuous, and consequently Q(z)=Q(z)⊤∀Q(z)∈∂(∇ℋ)(z).Q(z)=Q(z) ∀ Q(z)∈∂( )(z). (93) Let m⋆=minQ(z)∈∂(∇ℋ)(z)λmin(Q(z))m_ = _Q(z)∈∂( )(z) _ (Q(z)), and since z⋆z_ is a local minimum, by compactness we have that Q(z)≻0Q(z) 0 for all Q(z)∈∂(∇ℋ)(z)Q(z)∈∂( )(z) and consequently m⋆>0m_ >0. By upper semicontinuity of the generalized Jacobian there exists δ>0 >0 and a neighbourhood δN_ of the local minimum z⋆z_ such that ‖Q(z)−Q(z⋆)‖2≤δ∀z∈δ\|Q(z)-Q(z_ )\|_2≤ ∀ z _ (94) Fix μ⋆∈(0,m⋆) _ ∈(0,m_ ) such that δ<m⋆−μ⋆ <m_ - _ . Then is holds that Q(z)≻μ⋆Id∀Q(z)∈∂(∇ℋ)(z),∀z∈δQ(z) _ I_d ∀ Q(z)∈∂( )(z),\>∀ z _ (95) By Taylor expansion of ℋH at z⋆z_ in δN_ we obtain ℋ(z⋆) (z_ ) ≥ℋ(z)+∇ℋ(z)⊤(z⋆−z)+12(z⋆−z)⊤Q(z)(z⋆−z)Q(z)∈∂(∇ℋ)(z) (z)+ (z) (z_ -z)+ 12(z_ -z) Q(z)(z_ -z) Q(z)∈∂( )(z) ≥ℋ(z)+∇ℋ(z)⊤(z⋆−z)+μ⋆2∥z⋆−z∥22 (z)+ (z) (z_ -z)+ _ 2\|z_ -z\|_2^2 (96) where in the last passage we have applied (95). Re-ordering the terms, it holds that ℋ(z)−ℋ(z⋆) (z)-H(z_ ) ≤−∇ℋ(z)⊤(z⋆−z)−μ⋆2∥z⋆−z∥22 ≤- (z) (z_ -z)- _ 2\|z_ -z\|_2^2 ≤∇ℋ(z)⊤(z−z⋆)−μ⋆2∥z−z⋆∥22 ≤ (z) (z-z_ )- _ 2\|z-z_ \|_2^2 (97) The r.h.s. is a function ar−μ⋆r2/2ar- _ r^2/2 that reaches a maximum in r>0r>0 for r=a/μ⋆r=a/ _ . Taking now a=∇ℋa= and r=z−z⋆r=z-z_ we obtain ℋ(z)−ℋ(z⋆)≤12μ⋆‖∇ℋ(z)‖22.H(z)-H(z_ )≤ 12 _ \| (z)\|_2^2. (98) ∎ For a one-hidden-layer hybrid Hamiltonian, ℋ(z)=12‖z‖22−ℱ2(W21z),∇ℋ(z)=z−W21⊤Ψ2(W21z).H(z)= 12\|z\|_2^2-F_2(W_21z), (z)=z-W_21 _2(W_21z). (99) If Ψ2 _2 is locally Lipschitz, the Clarke sum and chain rules give the sound inclusion ∂(∇ℋ)(z⋆)⊆I−W21⊤QW21:Q∈∂(Ψ2)(W21z⋆).∂( )(z_ ) \I-W_21 QW_21:Q∈∂( _2)(W_21z_ ) \. (100) Thus the tractable sufficient condition infQ∈∂(Ψ2)(W21z⋆)λmin(I−W21⊤QW21)>0 _Q∈∂( _2)(W_21z_ ) _ \! (I-W_21 QW_21 )>0 (101) implies Proposition 5. The inclusion may be strict for nonsmooth compositions, but positivity over the displayed outer set remains a sound certificate. For the one-hidden-layer Hamiltonian (99), the gradient is the affine term z minus the composition of the locally Lipschitz activation Ψ2 _2 with two linear maps. The Clarke generalized-Jacobian sum and chain rules yield (100). Therefore, requiring every matrix I−W21⊤QW21I-W_21 QW_21 with Q∈∂(Ψ2)(W21z⋆)Q∈∂( _2)(W_21z_ ) to be uniformly positive definite is a directly checkable sufficient condition for Proposition 5. The following specializations make this condition explicit. Softmax. Let ℱ2(s)=1βlog(∑i=1N2eβsi),p(s)=Ψ2(s)=softmax(βs),β>0.F_2(s)= 1β ( _i=1^N_2e^β s_i ), p(s)= _2(s)=softmax(β s), β>0. (102) Then ∂(Ψ2)(s)=DΨ2(s)∂( _2)(s)=\D _2(s)\ with DΨ2(s)=β[diag(p(s))−p(s)p(s)⊤],D _2(s)=β [diag(p(s))-p(s)p(s) ], (103) and a sufficient local PL condition at z⋆z_ is I−βW21⊤[diag(p⋆)−p⋆p⋆⊤]W21≻0,p⋆=p(W21z⋆).I-β W_21 [diag(p_ )-p_ p_ ]W_21 0, p_ =p(W_21z_ ). (104) Softmax is smooth and bounded, so it is also compatible with the bounded-first-hidden-layer condition used to establish coercivity. Hyperbolic tangent. For the componentwise activation Ψ2,i(si)=tanh(βsi) _2,i(s_i)= (β s_i), one may take ℱ2(s)=∑i1βlogcosh(βsi),F_2(s)= _i 1β (β s_i), (105) with ∂(Ψ2)(s)=βdiag(sech2(βsi)).∂( _2)(s)= \ \! (sech^2(β s_i) ) \. (106) The generalized-Hessian certificate reduces to I−βW21⊤diag(sech2(β(W21z⋆)i))W21≻0.I-β W_21 diag\! (sech^2(β(W_21z_ )_i) )W_21 0. (107) The activation is smooth and bounded and therefore satisfies both the local regularity requirement and the first-hidden-layer boundedness requirement. Logistic sigmoid. For Ψ2,i(si)=σ(βsi) _2,i(s_i)= (β s_i), with σ(r)=(1+e−r)−1 (r)=(1+e^-r)^-1, ∂(Ψ2)(s)=βdiag(σ(βsi)(1−σ(βsi))),∂( _2)(s)= \ \! ( (β s_i)(1- (β s_i)) ) \, (108) and local PL follows from I−βW21⊤diag(σi⋆(1−σi⋆))W21≻0,σi⋆=σ(β(W21z⋆)i).I-β W_21 diag\! ( _i (1- _i ) )W_21 0, _i = (β(W_21z_ )_i). (109) Sigmoid is again smooth and bounded. ReLU and activation kinks. For the componentwise ReLU activation Ψ2,i(si)=[si]+ _2,i(s_i)=[s_i]_+, a convex C1,1C^1,1 primitive is ℱ2,i(si)=12[si]+2F_2,i(s_i)= 12[s_i]_+^2. Its Clarke generalized Jacobian is ∂(Ψ2)(s)=diag(q):qi=1,si>0,qi=0,si<0,qi∈[0,1],si=0.∂( _2)(s)= \diag(q): array[]lq_i=1,&s_i>0,\\ q_i=0,&s_i<0,\\ q_i∈[0,1],&s_i=0 array \. (110) Consequently, local PL is certified even at a kink whenever I−W21⊤DW21≻0for every D∈∂(Ψ2)(W21z⋆).I-W_21 DW_21 0 every D∈∂( _2)(W_21z_ ). (111) Let DmaxD_ set the slope to one for every active or zero preactivation and to zero for every strictly inactive preactivation. Since 0⪯D⪯Dmax0 D D_ for every admissible D, condition (111) is equivalent to the single worst-case check I−W21⊤DmaxW21≻0.I-W_21 D_ W_21 0. (112) This removes the need to assume a fixed activation pattern around z⋆z_ . ReLU remains unbounded, however, so it is not compatible with the bounded-first-hidden-layer hypothesis used above to establish global coercivity when it is chosen as Ψ2 _2; the local Clarke certificate is still applicable when ReLU appears only in deeper layers or when coercivity is guaranteed by another architectural mechanism. Bounded piecewise-linear activations. The generalized certificate also covers nonsmooth activations that are compatible with the hybrid coercivity argument. For example, let Ψ2,i(si)=clip(si,−1,1) _2,i(s_i)=clip(s_i,-1,1) (hard-tanh ). A convex C1,1C^1,1 primitive is ℱ2,i(si)=12si2,|si|≤1,|si|−12,|si|>1.F_2,i(s_i)= cases 12s_i^2,&|s_i|≤ 1,\\ |s_i|- 12,&|s_i|>1. cases (113) The generalized slopes are one for |si|<1|s_i|<1, zero for |si|>1|s_i|>1, and any value in [0,1][0,1] at si=±1s_i=± 1. Thus the same worst-case diagonal test as (112), with DmaxD_ defined from the unsaturated and kink coordinates, certifies local PL. Unlike ReLU, hard-tanh is bounded and therefore can simultaneously satisfy the first-hidden-layer boundedness condition used for coercivity. Deeper hybrid networks. For deeper architectures, the activation-independent statement remains Proposition 5: it is sufficient that every element of the full projected Clarke generalized Hessian at the learned equilibrium be positive definite. Computing that set exactly may be conservative or expensive because generalized-Jacobian chain rules through the feedforward composition can yield set-valued inclusions. Any tractable outer approximation ⋆Q_ satisfying ∂(∇ℋ)(z⋆)⊆⋆∂( )(z_ ) _ (114) therefore provides a sound sufficient certificate whenever every M∈⋆M _ has a common positive spectral margin. The one-hidden-layer formulas above are the simplest exact or directly computable instances of this principle. Appendix C Experimental details and reproducibility This section records the experimental choices used in the released code and complements the compact discussion in Section 4. The benchmark-specific architectures are intentionally allowed to vary: the common object across tasks is the pH-EBM construction, i.e., a learned Hopfield storage function coupled to free port-Hamiltonian interconnection, dissipation, and input maps. Unless stated otherwise, continuous-time rollouts are integrated with fourth-order Runge–Kutta (RK4), and all reported certificates are evaluated after training without modifying the learned parameters. Table 2: Prediction accuracy and certified robustness across nonlinear dynamical systems. Comparative performance of the pH-EBM across benchmarks addressing nonlinear system identification, nonconvex dynamics, and higher-dimensional controlled systems. Results reported within square brackets, [..][..], correspond to evaluations on distinct test sets of the same benchmark, while entries separated by semicolons refer to variants of the reference architecture. Benchmark Metric ↓ Published portHNN-u pH-EBM Silverbox RMSE [0.289,0.334,0.257][0.289,0.334,0.257] [51.343,50.981,41.138][51.343,50.981,41.138] [0.395,0.545,0.327][0.395,0.545,0.327] CED RMSE [0.062,0.047][0.062,0.047] [0.241,0.366][0.241,0.366] [0.073,0.055][0.073,0.055] Duffing double-well RMSE - 0.2540.254 0.0480.048 33-link pendulum RMSE 0.162;0.1620.162;0.162 0.051;0.0510.051;0.051 0.0140.014 NanoDrone RMSE - [1.532,5.260,2.150,15.541][1.532,5.260,2.150,15.541] [1.368,3.550,1.343,11.648][1.368,3.550,1.343,11.648] NanoDrone (Melon) RMSE [4.051,12.581,6.059,35.5735][4.051,12.581,6.059,35.5735] [3.939,12.445,5.574,30.802][3.939,12.445,5.574,30.802] [11.860,28.014,7.169,36.597][11.860,28.014,7.169,36.597] C.1 Datasets, splits, and preprocessing Silverbox. We use the canonical SISO Silverbox benchmark through the nonlinear_benchmarks interface. Input and output are standardized using the mean and standard deviation of the training record only; the same affine transformation is then applied unchanged to every test record. The final model uses a four-dimensional latent state. Its initial state is inferred from a burn-in window before the free rollout begins. The released evaluation scores all three official test records—multisine, full arrow, and arrow without extrapolation - rather than selecting a single favorable trajectory. The RMSE score reported in the main is the average RMSE over the three test datasets: check Table 2 for the individual scores. The champion is trained on windows of 192 prediction steps with stride 16 and batch size 32. The benchmark-native output RMSE is reported in mV, together with the dimensionless NRMSE used internally for diagnostics. (a) Projected safe-set geometry and input-radius sweep. (b) Affine Hamiltonian slice and profiled energy envelope. Figure 4: Silverbox: learned energy geometry and certificate diagnostics. The upper composite shows the profiled two-dimensional energy view at the nominal shell: traced boundary states are colored by their pointwise admissible-input radius, the red cross marks the least tolerant sampled state, and the adjacent panel resolves the same radii along the shell. The lower sweep compares the sampled shell minimum with the regular-shell estimate as the relative energy budget increases; green diamonds mark detected saddle-energy shells. The second composite contrasts an affine Hamiltonian slice with the profiled envelope in both 2D and 3D. Their close agreement indicates that the selected plane captures the dominant single-well geometry over the explored region. These plots visualize the learned Hamiltonian and should not be interpreted as recovery of a unique physical Silverbox potential. Coupled Electric Drive (CED). CED contains two SISO realizations corresponding to different actuation amplitudes. Both training realizations are concatenated and standardized with statistics computed jointly from the combined training data, while the two test records remain separate and are rolled out independently, with no state carried between them. The benchmark-provided state-initialization window has length 10. The selected configuration uses a four-dimensional latent state, 64-step training windows, unit stride, and batch size 16. We report each test RMSE in ticks/s and retain the mean NRMSE only as an aggregate diagnostic. This separation is important because the two records probe distinct operating amplitudes rather than repeated measurements of the same trajectory. (a) Projected safe-set boundary and epsilon sweep. (b) Affine slice and profiled Hamiltonian view. Figure 5: CED: learned energy geometry and certificate diagnostics. The upper composite shows the profiled safe-set boundary, colored by pointwise admissible-input radius, together with the shell-wise radius distribution and the least tolerant sampled state. The energy sweep compares the sampled minimum radius with the regular-shell estimate from the nominal shell into higher-energy regions. The second composite compares the affine slice and profiled envelope in 2D and 3D. The smooth, nearly monotone radius trend and close slice/profile agreement are consistent with an effectively single-well learned geometry over the explored region; they do not establish convexity of an underlying physical energy. Duffing double well. The Duffing experiment is state observed and therefore differs deliberately from the output-only identification benchmarks. Data are generated from the physical dynamics q˙=p q=p, p˙=q−q3−rp+u p=q-q^3-rp+u with r=0.4r=0.4, sampling interval Δt=0.02 t=0.02 s, and 20 s trajectories. The frozen dataset contains 200 training, 40 validation, and 100 nominal-test trajectories. Initial conditions alternate between the two wells, with 20% of the trajectories drawn from a higher-energy regime. Training inputs combine piecewise-constant, pseudo-random binary, and multisine forcing. Four held-out forcing families - step, chirp, unseen-frequency sinusoid, and long constant forcing - are generated only for OOD evaluation. Each OOD family contains 25 trajectories in the supplied artifact. No state or input normalization is used: the learned state is directly z=(q,p)z=(q,p), the physical initial condition is provided to the model, and the state itself is the prediction target. This makes the learned Hamiltonian directly comparable with the analytic mechanical energy. Figure 6: Duffing Hamiltonian reconstruction in physical coordinates. Analytic (left) and learned pH-EBM (right) energies are shown directly in (q,p)(q,p), so no latent projection is involved. Both the contour and 3D views recover the characteristic double-well topology and the central saddle separating the two basins. The comparison supports qualitative recovery of the physical energy geometry; it does not imply unique or exact identification of the underlying Hamiltonian. n-link pendulum. The pendulum experiments import the arrays produced by the reference Deep Dissipative Dynamics data generator; the pH-EBM code does not regenerate an “equivalent” dataset. Inputs, partial observations, and full states are loaded from the corresponding .npy files. The primary mode is output-only: the observation consists of the first-joint angle/velocity pair, whereas the complete physical state is retained only for diagnostics. All supplied trajectories use Δt=0.05 t=0.05 s and duration 10 s. The pH-EBM uses an eight-dimensional latent state, an input-aware 10-sample encoder, and the same output-only convention for both the two- and three-link experiments. The final released configurations use 50-step windows, stride 50, batch size 8, and 500 epochs. The supplementary figures report both nominal reconstruction/OOD diagnostics and the learned projected energy geometry. The supplied release contains complete champion artifacts for n=2n=2 and n=3n=3. (a) n=2n=2: prediction, OOD tests, and certificate diagnostics. (b) n=2n=2: affine slice versus profiled envelope. Figure 7: 22-link pendulum: prediction, OOD forcing, and certificate geometry. The first composite reports the nominal pulse input and the two official OOD inputs, representative trajectory reconstructions, per-trajectory RMSE across pH-EBM and dissipative/naive DDM baselines, the empirical safety margin of the evaluated OOD trajectories, the admissible-input-radius sweep across energy shells, and the projected shell colored by pointwise radius. The radius decreases near several detected saddle-energy shells before increasing beyond the nominal operating energy. The second composite compares an affine Hamiltonian slice with the profiled envelope, showing the nonconvex structure that is hidden by any single trajectory projection. (a) n=3n=3: prediction, OOD tests, and certificate diagnostics. (b) n=3n=3: affine slice versus profiled envelope. Figure 8: 33-link pendulum: prediction, OOD forcing, and certificate geometry. The first composite mirrors the 22-link protocol: nominal and official OOD inputs, representative reconstructions, per-trajectory RMSE, empirical safety margins, and an energy-shell sweep of the sampled and regular-shell admissible radii. Detected saddle-energy shells coincide with pronounced changes in the radius curve. The projected energy field is shown for geometric context, but the shell boundary is intentionally omitted because too few rays met the reconstruction threshold; it should therefore not be read as a complete projected boundary certificate. The second composite compares affine-slice and profiled-envelope views of the learned nonconvex Hamiltonian. NanoDrone. The NanoDrone loader follows a frozen external dataset revision and validates all file names, timestamps, and quaternion norms before training. The training families are Square, Random, and Chirp, with four runs per family; runs 1-3 form the development training split and run 4 is used for development validation. The three Melon trajectories form the held-out extrapolation family and are explicitly excluded from scaler fitting and from the S3 model-selection protocol. Samples are spaced by Δt=0.01 t=0.01 s. The physical state is represented by 12 coordinates: position, linear velocity, the SO(3)SO(3) logarithm of attitude, and angular velocity. The four inputs are rotor angular speeds. State and input scalers are fitted only on the training trajectories, and the primary pH-EBM is trained in state-observed mode with an identity readout. The final S3 confirmation campaign evaluates horizons of 50, 100, 250, 500, and 1000 steps; the primary long-horizon comparison is 500 steps (5 s), ten times the 50-step (0.5 s) training/evaluation horizon emphasized by the benchmark. Four frozen seeds (11, 23, 47, 89) are used in the paired pH/black-box confirmation protocol, with model selection performed on the training-supported families only. (a) Chirp long-horizon prediction and learned safe envelope. (b) Square long-horizon prediction and learned safe envelope. Figure 9: NanoDrone: training-supported long-horizon Chirp and Square rollouts. For each excitation, the top panels report channel-averaged position, velocity, and angular-velocity errors; the 3D panel compares the nominal trajectory with pH-EBM and black-box rollouts inside the projected learned-energy envelope; the exit-timing panel tracks the most critical projected coordinate; and the final panel reports physical MAE over the extended horizon. The pH-EBM remains inside the displayed envelope for the shown 5s5\,s rollouts, whereas the black-box trajectory exits at the indicated time. The envelope is a visualization of the high-dimensional certificate, not a substitute for its full-state evaluation. (a) Held-out Melon long-horizon prediction. (b) Projected Hamiltonian and profiled envelope. Figure 10: NanoDrone: held-out Melon extrapolation and learned Hamiltonian geometry. The first composite evaluates the unseen Melon family over a 5s5\,s rollout, reporting channel-wise error, projected trajectory evolution relative to the learned-energy envelope, exit timing, and physical MAE. The black-box prediction is initially competitive but develops rapidly growing long-horizon error and leaves the displayed envelope, whereas the pH-EBM remains bounded in the shown rollout. The second composite compares an affine two-dimensional Hamiltonian slice with the profiled envelope obtained by optimizing over unshown coordinates, in both contour and 3D views, exposing the nonconvex geometry of the learned high-dimensional energy. Figure 11: NanoDrone projected certificate. Left: projection of the traced full-state energy shell, with boundary samples colored by their pointwise admissible-input radius and the least tolerant sampled state marked explicitly. Right: the same pointwise radii ordered along the projected shell, together with the sampled shell minimum and regular-shell estimate. The projection is diagnostic only: certification is evaluated in the full model state space, and the displayed contour is not asserted to contain every high-dimensional boundary extremum. C.2 Per-benchmark model configurations The selected models share the same pH-EBM construction but not a forced common width. This is intentional: state dimension, observation structure, and input dimension differ substantially across the tasks. In all code-supported configurations below, J is free and skew-symmetric by construction, R is learned through a Cholesky factor and is positive semidefinite (strictly positive definite when the learned factor has nonzero diagonal), and G is learned rather than tied to the scalar Hamiltonian parameterization. “Active parameters” count only terms used by the selected computational graph. Table 3: Per-benchmark pH-EBM configurations. Architectures and active parameter counts reconstructed from the supplied champion configurations and frozen provenance artifacts. Entries marked “not in release” correspond to rows already present in the draft for which no experiment code or champion configuration was supplied. Benchmark dzd_z Hamiltonian architecture Activation J R G Parameters Silverbox 4 →64128\!→\!64 softmax(22)→ (44) free skew LL⊤L free G∈ℝ4×1G ^4× 1 ≈9.1≈ 9.1k CED 4 →4896\!→\!48 softmax(22)→ (44) free skew LL⊤L free G∈ℝ4×1G ^4× 1 ≈5.25≈ 5.25k Duffing double-well 2 →3264\!→\!32 tanh(2)→ (2) (44) free skew LL⊤L free G∈ℝ2×1G ^2× 1 ≈2.28≈ 2.28k n-link pendulum 8 →3264\!→\!32 tanh(2)→ (2) (44) free skew LL⊤L free G∈ℝ8×1G ^8× 1 ≈2.99≈ 2.99k NanoDrone† 12 →64128\!→\!64 tanh→ (44) free pH field learned PSD unrestricted 4-port 43,738 C.3 Training, validation, and model selection Across the standard pH-EBM interfaces, optimization is performed on finite rollout windows with a Huber observation loss, Adam/AdamW updates, global gradient clipping, and RK4 integration. The first training epochs use a stop-gradient rollout warm-up to limit long-horizon gradient pathologies; subsequent epochs differentiate through the complete rollout. This warm-up changes the optimization path only and does not alter the model evaluated at test time. Benchmark-specific schedules. The Silverbox champion uses 800 epochs, learning rate 4.59×10−44.59× 10^-4, 192-step windows, stride 16, and batch size 32. CED uses a longer 8000-epoch schedule, learning rate 10−410^-4, 64-step windows, stride 1, and batch size 16. Duffing uses 1000 epochs, learning rate 1.5×10−31.5× 10^-3, 100-step windows, stride 25, and batch size 64 in the frozen champion. The two supplied n-link champions use 500 epochs, learning rate 3×10−43× 10^-4, 50-step windows, stride 50, and batch size 8. For NanoDrone, the final S3 confirmation protocol fine-tunes a frozen parent model for 120 epochs with fresh Adam, batch size 64, a five-epoch warm-up, and a learning-rate schedule peaking at 3×10−53× 10^-5 and ending at 10−710^-7. Validation is evaluated at a fixed list of epochs, and the S3 campaign is repeated for four seeds. The Melon family is inaccessible to this selection stage by protocol. Initial conditions and readouts. Silverbox, CED, and the primary n-link experiments are output-only systems. Their latent initial condition is inferred from a short observation window; CED and n-link use an input-aware encoder so the burn-in actuation is available when estimating z(0)z(0). Their measured output is distinct from the power-conjugate port output. Duffing and NanoDrone instead use observed-state training: the physical state initializes the rollout directly and no learned state decoder is required. This distinction is preserved deliberately because the corresponding datasets expose different information to the learner. C.4 Baselines and comparison protocol The baseline is chosen to answer the scientific question posed by each experiment rather than forcing one architecture across unrelated datasets. Silverbox and CED are compared with benchmark-native literature references under the same published test records and physical error units. Duffing uses the released portHNN-u implementation, including a version with free full-matrix J,R,GJ,R,G, so that the comparison isolates the effect of the Hamiltonian parameterization rather than fixing the pH operators to a scalar ansatz. The n-link comparison imports the exact upstream Deep Dissipative Dynamics trajectories and uses its reference model on those same arrays. NanoDrone uses the accompanying black-box neural model and a paired confirmation protocol: pH and black-box models see the same S3 families, horizons, validation starts, and seed set. The frozen protocol records 43,738 active parameters for the pH model and 18,508 for the reference black-box model; the comparison is therefore a benchmark-reference comparison rather than a parameter-matched capacity ablation. Missing certification for an unconstrained baseline is denoted as unavailable, not as a zero robustness radius. C.5 Prediction metrics and statistical reporting For scalar-output identification benchmarks we report the benchmark-native root-mean-square error RMSE=1N∑i=1N(y^i−yi)2,RMSE= 1N _i=1^N( y_i-y_i)^2, (115) and, for internal cross-dataset diagnostics, NRMSE=RMSE/std(y)NRMSE=RMSE/std(y). Silverbox RMSE is converted back to mV and CED is reported in ticks/s. Duffing and the state-observed diagnostics use the same definitions directly in physical state coordinates. NanoDrone is evaluated in physical units after undoing the training scalers. Its primary paired S3 estimand is the difference between pH-EBM and black-box physical trajectory MAE at 500 steps, with additional reporting at 50, 100, 250, and 1000 steps and channel-family diagnostics for position, velocity, orientation, and angular velocity. Quaternion observations are converted to SO(3)SO(3) logarithmic coordinates for the 12-state model; orientation error is additionally tracked geometrically in the evaluation code. For consistency with official validation script, Nanodrone errors are also computed as per-output channel RMSE. C.6 Operational energy-shell selection The certificate calculations are performed after fitting. We first search the learned Hamiltonian for low-energy wells by multi-start gradient descent. Half of the starts are drawn from states visited by the fitted model when such states are available, and the remaining starts are randomized around their empirical center and scale. Descents are restricted to the operational state ball used by the numerical model, and terminal points are clustered into candidate wells. Candidates with a large energy-gradient norm are rejected before they are used as stationary anchors. This avoids treating a low-energy but nonstationary optimizer endpoint as a physical attractor. For a chosen threshold ϵε, the full-state shell ℋ(z)=ϵH(z)=ε is traced from every discovered well whose energy lies below the threshold. Random unit directions are launched from the well; along each ray the algorithm expands an outer bracket until the first crossing of the target energy is found and then refines the crossing by bisection. A short damped correction in the local gradient-normal direction finally reduces the residual |ℋ(z)−ϵ||H(z)-ε|. Rays that do not reach the shell within the configured search radius are recorded as failures rather than silently projected to an unrelated point. Consequently, all extrema reported from this procedure are sampled shell estimates; a finite ray set by itself is not presented as an exhaustive proof of a continuous high-dimensional boundary. C.7 Numerical computation and verification of certificates Hamiltonian projection. A two-dimensional picture of a high-dimensional Hamiltonian is necessarily a reduction, so we use two complementary views and label them separately. Let c∈dc∈ ^d be the selected well, let U∈d×2U∈ ^d× 2 contain two orthonormal display directions, and let V span the orthogonal complement. The affine slice Hslice(ζ)=ℋ(c+Uζ)H_slice(ζ)=H(c+Uζ) (116) is an exact cross-section: all unshown coordinates are frozen. The profiled envelope H~(ζ)=minωℋ(c+Uζ+Vω) H(ζ)= _ωH(c+Uζ+Vω) (117) instead minimizes over the unshown coordinates. Up to numerical optimization error, the set ζ:H~(ζ)≤ϵ\ζ: H(ζ)≤ε\ is the planar projection of the full sublevel set z:ℋ(z)≤ϵ\z:H(z)≤ε\. The minimization is solved independently over the display grid with Adam and multiple warm starts, including the orthogonal coordinates of discovered wells. We retain the orthogonal-gradient norm as a stationarity diagnostic for the profile. The display plane is chosen from physically or dynamically informative directions. With multiple discovered wells, the first direction follows the inter-well geometry and the second is taken from the dominant orthogonal variation of the sampled boundary/trajectory cloud. With a single well, the default is PCA of the sampled safety boundary; ordinary trajectory PCA is used only when boundary information is unavailable. Duffing is already two dimensional and therefore requires no projection. The plane limits are derived from robust quantiles of the projected reference states, boundary samples, and discovered wells, with padding added only for readability. The supplementary slice/profile pairs in Figs. 4, 5, 7–8, and 9, 10 use this procedure. Input-radius computation. At every sampled shell state we evaluate the raw Hamiltonian gradient and the learned R and G matrices, then compute the pointwise boundary margin described in Section 3.1. For a chosen input norm, the sampled “exact” shell radius is the minimum of the pointwise ratio dH(z)−σ(z)(∇ℋ(z))∥G(z)⊤∇ℋ(z)∥∗ d_H(z)- _W(z)( (z))\|G(z) (z)\|_* (118) over the traced states for which the input couples to the energy gradient. In parallel, we compute the regular-boundary lower estimate from the sampled extrema of ‖∇H‖\|∇ H\|, normal dissipation, normal input gain, and disturbance support. The implementation records both values because they answer different questions: the first measures the least favorable sampled direction directly, while the second exposes the geometric factors entering the analytic lower bound. Both remain numerical shell estimates unless the continuous extrema are independently bounded. Epsilon sweeps. The epsilon sweeps in the supplementary figures do not rescale a single boundary. For every energy level, the full-state shell is traced again from scratch and both radii are recomputed. By default, 300 levels are sampled in relative energy coordinates ϵ−Hminε-H_ , from 0.100.10 to 5/35/3 times the nominal energy gap ϵnom−Hmin _nom-H_ . Thus the sweep covers the certified operating shell and extends by two thirds of its nominal gap to reveal what happens near and beyond critical levels. The boundary search radius is enlarged by a factor 1.5 during this sweep so that outer shells are not artificially truncated. For Duffing, the range crosses the saddle energy; after components merge, the shell is treated globally rather than assigning points to their original well. Failed levels are stored as missing values and are not interpolated. Every sweep artifact stores the absolute ϵε, HminH_ , relative energy, sampled radius, regular-boundary estimate, number of successful boundary points, failed-ray fraction, and component summaries so the plotted curves can be audited independently. Hamiltonian-normalized empirical certificate. The robustness radius ρϵ,⋆ _ε, introduced above is an exact certificate defined at a specific energy level ϵε and expressed in the native coordinates of the model input. As a consequence, its numerical value may vary under simple rescalings of either the energy function or the input variables, making direct comparisons across architectures difficult. To obtain a dimensionless quantity suitable for cross-model evaluation, we introduce an empirical normalization procedure based exclusively on validation data. Common energy level. Given a fixed validation rollout (zi,ui)i=1Nval,\(z_i,u_i)\_i=1^N_val, we define the reference energy level as the q-quantile of the validation energy distribution, ϵq:=Qq(ℋ(zi)i=1Nval),q=0.99. _q:=Q_q\! (\H(z_i)\_i=1^N_val ), q=0.99. (119) This choice ensures that all models are evaluated on a shell corresponding to the same validation-state occupancy rather than at an arbitrary absolute energy value. Input normalization. To remove the dependence on the physical units and scaling of the inputs, we normalize the input space using the empirical second moment Mu:=1Nval∑i=1Nvaluiui⊤=LuLu⊤,M_u:= 1N_val _i=1^N_valu_iu_i =L_uL_u , (120) and introduce normalized coordinates u=Luv.u=L_uv. When MuM_u is rank-deficient, LuL_u is restricted to the subspace effectively excited by the validation data. In these coordinates, distances are measured relative to the variability observed during validation rather than in the original input units. Normalized robustness margin. Using the runtime vector field fθrunf_θ^run, define the directional energy dissipation rate ℓH(z,u):=∇ℋ(z)⊤fθrun(z,u). _H(z,u):= (z) f_θ^run(z,u). (121) A negative value of ℓH _H corresponds to local energy dissipation, whereas ℓH≥0 _H≥ 0 identifies inputs capable of halting or reversing this decrease. We therefore define the normalized robustness margin at state z as rH(z):=inf‖v‖2:ℓH(z,Luv)≥0.r_H(z):= \\|v\|_2: _H(z,L_uv)≥ 0 \. (122) Equivalently, rH(z)r_H(z) is the smallest perturbation, measured in normalized input coordinates, capable of eliminating the local decrease of the energy function. By convention, rH(z)=0r_H(z)=0 whenever ℓH(z,0)≥0 _H(z,0)≥ 0. For input-affine dynamics, fθrun(z,u)=fθrun(z,0)+Bθ(z)u,f_θ^run(z,u)=f_θ^run(z,0)+B_θ(z)u, the above optimization admits the closed-form expression rH(z)=−∇ℋ(z)⊤fθrun(z,0)∥Lu⊤Bθ(z)⊤∇ℋ(z)∥2,r_H(z)= - (z) f_θ^run(z,0) \|L_u B_θ(z) (z) \|_2, (123) whenever the denominator is nonzero. For nonlinear forcing structures, rH(z)r_H(z) is obtained numerically by locating the first zero crossing of ℓH _H along the normalized input direction. Empirical normalized radius. Let Γ^q _q denote a numerical reconstruction of the energy shell Γq:=z:ℋ(z)=ϵq _q:=\z:H(z)= _q\ obtained from validation data. We then define ρ^q,norm:=minz∈Γ^qrH(z),q=0.99. ρ_q,norm:= _z∈ _qr_H(z), q=0.99. (124) The quantity ρ^q,norm _q,norm measures the smallest normalized input perturbation required to destroy local energy dissipation on the validation shell. By construction, it is invariant under any strictly increasing reparameterization of ℋH and uses a common, data-derived input geometry for all models. Since the shell itself is reconstructed numerically from finite validation trajectories, ρ^0.99,norm _0.99,norm should be interpreted as an empirical benchmarking metric for cross-model comparison. In contrast, the radius ρϵ,⋆ _ε, remains the exact continuous-state robustness certificate associated with a prescribed energy level. C.8 Additional benchmark results Figures 4-11 collect the supplementary outputs generated by the released experiment pipeline. In particular, the Silverbox and CED plates show that the same certificate machinery can be applied when the learned energy is effectively single-well (or almost); the Duffing figure validates recovery of a genuinely nonconvex physical energy; the n-link figures compare prediction/OOD behavior with high-dimensional projected geometry; and the NanoDrone figures expose the distinction between short-horizon fit and long-horizon safe recursive deployment. These visualizations are diagnostic companions to the quantitative table in the main paper and do not replace the full-state certificate computations from which the shell radii are obtained. C.9 Ablations and sensitivity analyses The released experiments prioritize ablations that alter the scientific certificate rather than exhaustive optimizer sweeps. Duffing compares the modern-Hopfield Hamiltonian with the portHNN-u reference while keeping the port-Hamiltonian interpretation explicit. The n-link family changes mechanical dimension under a common eight-dimensional latent pH-EBM template. NanoDrone pairs the structured model with the benchmark black-box model over identical S3 starts and horizons. Finally, the epsilon sweeps provide a post-training sensitivity analysis of the certificate itself: they reveal how the minimum shell gradient, normal dissipation, input exposure, and admissible-input radius change as the certified operating region expands toward critical energy levels.