Paper deep dive
PHAST: Port-Hamiltonian Architecture for Structured Temporal Dynamics Forecasting
Shubham Bhardwaj, Chandrajit Bajaj
Intelligence
Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 90%
Last extracted: 7/20/2026, 10:59:00 PM
Summary
The paper introduces PHAST (Port-Hamiltonian Architecture for Structured Temporal dynamics), a neural architecture designed to forecast dissipative dynamical systems from position-only observations. PHAST decomposes the Hamiltonian into potential, mass, and damping components, allowing for structured learning across KNOWN, PARTIAL, and UNKNOWN knowledge regimes. It guarantees passivity and long-horizon stability through a port-Hamiltonian framework and Strang splitting integration, achieving superior forecasting and parameter recovery compared to baselines like HNNs, S5, and LinOSS.
Entities (11)
Relation Signals (10)
PHAST → decomposes → Hamiltonian
confidence 95% · PHAST... decomposes the Hamiltonian into potential V(q), mass M(q), and damping D(q)...
PHAST → uses → Port-Hamiltonian Framework
confidence 95% · We introduce PHAST (Port-Hamiltonian Architecture for Structured Temporal dynamics)... The port-Hamiltonian framework makes the conservative-dissipative split explicit...
PHAST → uses → Strang Splitting
confidence 92% · uses efficient low-rank PSD/SPD parameterizations, and advances dynamics with Strang splitting.
Port-Hamiltonian Framework → guarantees → Passivity
confidence 90% · guaranteeing dH/dt<=0 when R>=0... the system is passive by construction
PHAST → supports → KNOWN Regime
confidence 90% · uses efficient low-rank PSD/SPD parameterizations... across three knowledge regimes (KNOWN, PARTIAL, UNKNOWN)
PHAST → supports → PARTIAL Regime
confidence 90% · uses efficient low-rank PSD/SPD parameterizations... across three knowledge regimes (KNOWN, PARTIAL, UNKNOWN)
PHAST → supports → UNKNOWN Regime
confidence 90% · uses efficient low-rank PSD/SPD parameterizations... across three knowledge regimes (KNOWN, PARTIAL, UNKNOWN)
PHAST → outperforms →
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:Real physical systems are dissipative -- a pendulum slows, a circuit loses charge to heat -- and forecasting their dynamics from partial observations is a central challenge in scientific machine learning. We address the \emph{position-only} (q-only) problem: given only generalized positions~$q_t$ at discrete times (momenta~$p_t$ latent), learn a structured model that (a)~produces stable long-horizon forecasts and (b)~recovers physically meaningful parameters when sufficient structure is provided. The port-Hamiltonian framework makes the conservative-dissipative split explicit via $\dot{x}=(J-R)\nabla H(x)$, guaranteeing $dH/dt\le 0$ when $R\succeq 0$. We introduce \textbf{PHAST} (Port-Hamiltonian Architecture for Structured Temporal dynamics), which decomposes the Hamiltonian into potential~$V(q)$, mass~$M(q)$, and damping~$D(q)$ across three knowledge regimes (KNOWN, PARTIAL, UNKNOWN), uses efficient low-rank PSD/SPD parameterizations, and advances dynamics with Strang splitting. Across thirteen q-only benchmarks spanning mechanical, electrical, molecular, thermal, gravitational, and ecological systems, PHAST achieves the best long-horizon forecasting among competitive baselines and enables physically meaningful parameter recovery when the regime provides sufficient anchors. We show that identification is fundamentally ill-posed without such anchors (gauge freedom), motivating a two-axis evaluation that separates forecasting stability from identifiability.
Tags
Links
- Source: https://arxiv.org/abs/2602.17998v1
- Canonical: https://arxiv.org/abs/2602.17998v1
Trouble viewing inline? Open PDF directly →
Full Text
184,862 characters extracted from source content.
Expand or collapse full text
PHAST: Port-Hamiltonian Architecture for Structured Temporal Dynamics Forecasting Shubham Bhardwaj Chandrajit Bajaj Abstract Real physical systems are dissipative—a pendulum slows, a circuit loses charge to heat—and forecasting their dynamics from partial observations is a central challenge in scientific machine learning. We address the position-only (q-only) problem: given only generalized positions qtq_t at discrete times (momenta ptp_t latent), learn a structured model that (a) produces stable long-horizon forecasts and (b) recovers physically meaningful parameters when sufficient structure is provided. The port-Hamiltonian framework makes the conservative–dissipative split explicit via x˙=(J−R)∇H(x) x=(J-R)∇ H(x), guaranteeing dH/dt≤0dH/dt≤ 0 when R⪰0R 0. We introduce PHAST (Port-Hamiltonian Architecture for Structured Temporal dynamics), which decomposes the Hamiltonian into potential V(q)V(q), mass M(q)M(q), and damping D(q)D(q) across three knowledge regimes (KNOWN, PARTIAL, UNKNOWN), uses efficient low-rank PSD/SPD parameterizations, and advances dynamics with Strang splitting. Across thirteen q-only benchmarks spanning mechanical, electrical, molecular, thermal, gravitational, and ecological systems, PHAST achieves the best long-horizon forecasting among competitive baselines and enables physically meaningful parameter recovery when the regime provides sufficient anchors. We show that identification is fundamentally ill-posed without such anchors (gauge freedom), motivating a two-axis evaluation that separates forecasting stability from identifiability. Port-Hamiltonian systems, physics-informed machine learning, dissipative dynamics, structure-preserving neural networks 1 Introduction Forecasting the future states of dynamical systems from observational data is a central problem in scientific machine learning. Physical systems span a wide spectrum of complexity: simple conservative systems (e.g., an undamped pendulum or a frictionless spring) preserve total energy and trace closed orbits in phase space; dissipative systems (e.g., a pendulum slowing due to air drag, an RLC circuit losing charge to resistive heating) lose energy over time and converge toward attractors; and chaotic systems (e.g., a double pendulum, turbulent fluid flows) exhibit sensitive dependence on initial conditions, making long-horizon prediction fundamentally difficult. A practical forecasting framework must handle all three regimes—and, crucially, most real-world systems are dissipative: a swinging pendulum eventually stops, a robot arm loses energy through joint friction, and electrical circuits dissipate energy as heat. Problem setting. We formalize the forecasting task in the position-only (q-only) setting: we observe only generalized positions qtq_t at discrete times; momenta ptp_t are latent and never measured. The goal is to learn a structured dynamical model from q-only data that (a) produces stable long-horizon open-loop forecasts and (b) when partial physics is available, recovers physically interpretable parameters (potential, mass, damping). Port-Hamiltonian framework. The port-Hamiltonian framework (Van Der Schaft & Jeltsema, 2014) provides a principled decomposition for dissipative systems: x˙=(J−R)∇H(x)+Gu, x=(J-R)∇ H(x)+Gu, (1) where x=(q,p)x=(q,p) denotes generalized coordinates and momenta. The term GuGu represents external actuation: u is a generalized force (e.g., torque at a joint), and the port matrix G selects which degrees of freedom are actuated. The conjugate port output yport=G⊤∇H(x)y^port=G ∇ H(x) is velocity-like for mechanical systems, and the product u⊤yportu y^port is the instantaneous power supplied to the system—the rate at which external work flows in through the port. Our forecasting benchmarks have no external actuation, so we set u≡0u≡ 0 throughout Sections 1–4; the forced form and its port structure become central in the Energy–Casimir control study (Sec. 3.5). In port-Hamiltonian form, dynamics decompose into energy geometry (via the mass tensor M(q)M(q)), environmental structure (via the potential V(q)V(q)), and dissipation (D(q)D(q)), enabling principled combinations of known physics and learned components (Sec. 1). For mechanical systems, the Hamiltonian takes the form H(q,p)=V(q)+T(q,p)H(q,p)=V(q)+T(q,p) with kinetic energy T(q,p)=12p⊤M(q)−1pT(q,p)= 12p M(q)^-1p, where M(q)≻0M(q) 0 is the generalized mass (inertia) tensor—it induces a Riemannian metric on configuration space under which free trajectories are geodesics—and we define the generalized velocity v:=∇pH(q,p)=M(q)−1pv:= _pH(q,p)=M(q)^-1p. In our experiments we use a separable constant-mass approximation M(q)≈M(q)≈ M for efficiency (Sec. 1), but the architecture supports the general configuration-dependent case (Appendix H.2). The operator J=−J⊤J=-J encodes conservative energy exchange via the canonical symplectic structure. Dissipation is introduced through a structured operator R(q)=(000D(q)),D(q)=D(q)⊤⪰0,R(q)= pmatrix0&0\\ 0&D(q) pmatrix, D(q)=D(q) 0, so that damping acts only on the momentum dynamics. This corresponds to standard mechanical dissipation, where friction affects generalized velocities but not positions. With the canonical symplectic structure and the block-diagonal dissipation, the dynamics decompose explicitly as q˙ q =∇pH(q,p), = _pH(q,p), (2) p˙ p =−∇qH(q,p)−D(q)∇pH(q,p). =- _qH(q,p)-D(q)\, _pH(q,p). (3) This structure yields an explicit energy balance: H˙(x)=∇H(x)⊤(J−R)∇H(x)=−∇H(x)⊤R∇H(x)≤0, H(x)=∇ H(x) (J-R)∇ H(x)=-∇ H(x) R∇ H(x)≤ 0, (4) since ∇H⊤J∇H=0∇ H J∇ H=0 by skew-symmetry. As a result, the system is passive by construction, guaranteeing non-increasing energy regardless of the specific parameterization of H(x)H(x). Oscillatory dynamics. In the conservative case (D=0D=0), the dynamics reduce to x˙=J∇H(x) x=J∇ H(x); linearizing around an elliptic equilibrium (e.g., a local minimum of the potential with M≻0M 0) yields purely imaginary eigenvalues (phase-space rotations), corresponding to oscillatory exchange between kinetic and potential energy. Adding D⪰0D 0 shifts eigenvalues into the closed left half-plane (e.g., −σ±iω-σ± iω), yielding damped oscillations; Appendix A.4 provides a concise linear-algebra statement and proof for the linearized model x˙=(J−R)Qx x=(J-R)Qx (Mehl et al., 2018). This work develops neural architectures that respect this structure (Appendix A.2). Motivation. We are motivated by three converging observations: (i) Physics-informed networks like HNNs (Greydanus et al., 2019) and LNNs (Cranmer et al., 2020) show that structure improves generalization, but they assume purely conservative systems (R=0R=0). (i) State-space models like S5 (Smith et al., 2023), LinOSS (Hasani et al., 2024), and D-LinOSS (Boyer et al., 2025) achieve efficiency but typically do not provide explicit physical guarantees—e.g., they do not expose an explicit port-Hamiltonian decomposition or a passivity certificate. (i) Structure-preserving networks like Volume-Preserving Transformers (Brantner et al., 2024) enforce geometric constraints such as volume preservation, i.e., their one-step map ϕΔt:xt↦xt+1 _ t:x_t x_t+1 has unit Jacobian determinant det(∂ϕΔt/∂xt)=1 \! (∂ _ t/∂ x_t )=1, which is naturally aligned with conservative, divergence-free dynamics. However, they do not provide an explicit dissipation mechanism or a certificate of passivity for lossy systems, and strictly volume-preserving maps cannot represent generic dissipative flows that contract phase-space volume. We propose PHAST (Port-Hamiltonian with Strang Splitting), which properly models both the conservative (J) and dissipative (R) components. The key insight is that the port-Hamiltonian structure provides passivity regardless of how the components are parameterized: whether the potential V(q)V(q), mass M(q)M(q), and damping D(q)D(q) are structured functions or neural networks, the continuous-time dynamics satisfy dH/dt≤0dH/dt≤ 0 by construction. PHAST addresses the q-only problem via a three-stage pipeline: a causal velocity observer reconstructs q˙ q from position history, a canonicalizer maps to phase state (q,p^)(q, p), and the port-Hamiltonian core integrates forward with Strang splitting. Given a short burn-in window q0:K−1q_0:K-1, PHAST rolls out open-loop; Section 3 details the components. Knowledge regimes and the Hamiltonian template. A central challenge in learning dynamics from data is that trajectories alone do not uniquely identify physical parameters: multiple combinations of mass, potential, and damping can generate identical motions (gauge freedom). PHAST addresses this by allowing problem-specific structure to be imposed on the Hamiltonian components, breaking gauge freedom when partial physics knowledge is available. We formalize this flexibility through three knowledge regimes: Known, Partial, and Unknown, which differ only in which components of the Hamiltonian are specified versus learned. Anchoring either the potential V(q)V(q) or the mass matrix M(q)M(q) is sufficient to render the remaining parameters identifiable, while leaving both free leads to an underdetermined problem. Table 1 illustrates representative systems supported by PHAST, spanning classical mechanics, robotics, molecular dynamics, electrical circuits, thermodynamics, and ecological modeling. Our experiments span all six domains in Table 1: the port-Hamiltonian formulation is not restricted to mechanics—any system whose dynamics admit a storage function and a conservative–dissipative split can be cast in this form. Regime Domain System Forecast q V(q)V(q): landscape M(q)M(q): geometry D(q)D(q): dissipation Ident.? Known Mechanics Pendulum Angle θ(t)θ(t) −gcosθ-g θ Scalar m Air drag d(θ)d(θ) Yes Mechanics Spring–mass Displacement x(t)x(t) 12kx2 12kx^2 Scalar m Viscous bx˙b x Yes Molecular Lennard–Jones Positions ri(t)r_i(t) 4ϵ[(σ/r)12−(σ/r)6]4ε[(σ/r)^12-(σ/r)^6] Diagonal mImI Solvent friction γIγ I Yes Electrical RLC circuit Charge qC(t)q_C(t) q2/2Cq^2/2C Inductance L Resistance R Yes Thermal Heat exchange Temperature T(t)T(t) 12cT2 12c\,T^2 Time constant τ Heat loss κ Yes Partial Mechanics Pendulum (unk. m) Angle θ(t)θ(t) −gcosθ-g θ (given) mImI (learn m) Learned (bounded) Yes Robotics Robot arm Joint angles θi(t) _i(t) Gravity regressor CAD inertia template Joint friction (learned) Yes Astrophysics N-body gravity Positions ri(t)r_i(t) −∑i<jGmimj/rij-\! _i<j\!Gm_im_j/r_ij miIm_iI (learn mim_i) Drag (learned) Yes† General Any partially known System-dep. q(t)q(t) Template +ΔV+ V m0I+UU⊤m_0I+U Learned (bounded) Partial Unknown Any Black-box dynamics System-dep. q(t)q(t) Neural network Neural SPD net Neural PSD net No Ecology Predator–prey Densities ni(t)n_i(t) Neural network Neural SPD net Neural PSD net No Table 1: Hamiltonian decomposition across knowledge regimes. Each system is described through the three components of the port-Hamiltonian template: the energy landscape V(q)V(q) (environmental forces and equilibria), the geometry M(q)M(q) (Riemannian metric governing the cost of motion), and the dissipation D(q)D(q) (energy loss structure). The regime determines which components are specified versus learned. Anchoring either V(q)V(q) or M(q)M(q) is sufficient to break gauge freedom and render the remaining parameters identifiable; leaving both unconstrained (Unknown) yields an underdetermined inverse problem even when forecasting is accurate. Our experiments cover all six domains; the port-Hamiltonian formulation applies to any system admitting a storage function and a conservative–dissipative split. †Identifiable up to an overall mass scale; fixing the gravitational constant G resolves the remaining one-parameter gauge. Interpretive structure of the Hamiltonian. Appendix B.5 provides the complete port-Hamiltonian mapping (state, potential, mass, damping, and identifiability analysis) for every system in Table 1, including derivations and schematic diagrams. At a high level, the table illustrates that diverse systems—mechanical, electrical, thermal, and ecological—admit a common decomposition into three components with distinct physical roles: 1. Potential V(q)V(q): encodes the environmental energy landscape—external forces, equilibrium structure, and constraints (e.g., gravity for a pendulum, capacitive energy for an RLC circuit, entropy potential for a thermal system). 2. Mass M(q)M(q): defines the intrinsic geometry of motion. The kinetic energy T=12q˙⊤M(q)q˙T= 12 q M(q) q induces a Riemannian metric on configuration space Q; free (unforced, undamped) trajectories are geodesics under this metric. M(q)M(q) determines how effort translates to motion and encodes inertial coupling between degrees of freedom. 3. Damping D(q)D(q): encodes dissipative structure—how and where energy is lost, governing contraction toward attractors. The three knowledge regimes differ only in which of these components are specified versus learned; the compositional template—and its passivity guarantee—remains unchanged (Fig. 1). Concrete instantiations. To make the mapping from physical system to (V,M,D)(V,M,D) template explicit, we walk through three representative entries from Table 1 (Appendix B.5 provides the full derivation for every row): • Simple pendulum (KNOWN). State q=θq=θ (angle), momentum p=mℓ2θ˙p=m ^2 θ (angular momentum). The potential V(θ)=−mgℓcosθV(θ)=-mg θ encodes gravity; mass M=mℓ2M=m ^2 is a known scalar inertia; the only unknown is position-dependent air drag D(θ)D(θ), which PHAST learns. Because both V and M are given, the inverse problem is fully identifiable. • Cart-pole (PARTIAL). State q=(cart position,θ)q=(cart position,\,θ), a two-DOF system. Gravity gives V(q)=mgℓcosθV(q)=mg θ (template provided). The mass tensor M(q)M(q) is a 2×22×2 configuration-dependent inertia coupling cart and pole (given from the Lagrangian or CAD model). Dissipation D(q)D(q) has separate friction channels for the cart rail and pivot joint—learned with bounded strength. • RLC circuit (KNOWN, non-mechanical). State q=qCq=q_C (charge on capacitor), momentum p=Lip=Li (flux linkage, where i is current and L is inductance). Capacitive energy gives V(q)=q2/2CV(q)=q^2/2C; inductance plays the role of mass (M=LM=L); resistance dissipates energy as heat (D=RD=R). The same (V,M,D)(V,M,D) template applies, illustrating that the port-Hamiltonian decomposition extends beyond mechanics. UNKNOWNV(q)V(q): learnedM(q)M(q): learned (SPD)D(q)D(q): learned (PSD) PARTIALV(q)=V¯+εV~V(q)= V+ [rgb]0,0,0.7 [named]pgfstrokecolorrgb0,0,0.7 VM(q)M(q): givenD(q)D(q): learned (bounded) KNOWNV(q)V(q): givenM(q)M(q): givenD(q)D(q): learned (PSD) Shared template (all regimes):• Port-Hamiltonian state x=(q,p)x=(q,p) with conservative–dissipative split• Dissipation acts on momentum only (D(q)⪰0D(q) 0); mass is SPD (M(q)≻0M(q) 0)• Structure-preserving time stepping (Strang-split symplectic core) ⟵ less prior knowledge more prior knowledge ⟶ Figure 1: PHAST unifies three knowledge regimes under one port-Hamiltonian template. Regimes differ only in what is given vs. learned for (V,M(q))(V,M(q)); dissipation D(q)D(q) is always learned. All regimes share the same continuous-time structure x˙=(J−R)∇H x=(J-R)∇ H and a structure-preserving discrete-time transition map ΦΔt _ t (Strang splitting with a symplectic core), ensuring definiteness constraints by construction. (A) Windy Pendulumwindmmℓ θ (θ)d(θ)ggd(θ)=d0+Δd|sinθ|d(θ)=d_0+ d\,|\! θ|(B) Windy Double Pendulumwindm1m_1θ1 _1d1(θ1)d_1( _1)m2m_2θ2 _2d2(θ2)d_2( _2)ggM(q)M(q)di(θi)=bi+Δbi|sinθi|d_i( _i)=b_i+ b_i\,|\! _i|(C) Windy Cart-PolewindMMxxmmℓ θ (θ)d(θ)dcd_cggM(q)M(q)D(q)=diag(dc,d0+Δd|sinθ|)D(q)=diag( [rgb]0.90234375,0.46875,0.1171875 [named]pgfstrokecolorrgb0.90234375,0.46875,0.1171875d_c,\; [rgb]0.70703125,0.1953125,0.1953125 [named]pgfstrokecolorrgb0.70703125,0.1953125,0.1953125d_0+ d\,|\! θ|) (D) Quantitative comparison W. Pendulum W. Double Pend. W. Cart-Pole Model MSE100 RD2R^2_D MSE100 RD2R^2_D MSE100 RD2R^2_D PHAST (ours) KNOWN 0.064 1.001.00 0.532 1.001.00 0.033 −-1.10 PARTIAL 0.0560.056 0.18 0.5260.526 −-1.09 0.0140.014 −0.24-0.24 UNKNOWN 0.737 ≪0 \!0 0.866 ≪0 \!0 0.243 ≪0 \!0 S5 0.954 — 2.910 — 5.145 — LinOSS 1.536 — 3.339 — 4.167 — VPT 3.227 — 2.464 — 2.103 — Figure 2: Benchmark suite overview. (A–C) Schematics of the three representative dissipative systems with position-dependent damping: Windy Pendulum, Windy Double Pendulum (configuration-dependent M(q)M(q)), and Windy Cart-Pole (ℝ×1R×S^1, per-DOF damping, M(q)M(q)). The Cart-Pole features two distinct damping mechanisms: constant viscous cart friction (dcd_c) and position-dependent angular wind damping (d(θ)d(θ)). (D) Wrapped-angle rollout MSE at H=100H=100 (lower== better) and damping identifiability RD2R^2_D (higher== better; bold== best per column). PHAST (KNOWN) achieves near-perfect identifiability (RD2≈1R^2_D≈ 1) on Pendulum and Double Pendulum; PHAST (PARTIAL) achieves the best forecasting across all three systems. All models use 3k–9k parameters (exact counts vary by state dimension; see Appendix). Baselines do not expose explicit damping fields (—). Appendix B describes the full suite of thirteen benchmarks, including conservative variants, the harmonic oscillator, and five non-mechanical systems (Tables 20–24). A critical but often overlooked issue emerges in practice: forecasting accuracy and parameter recovery are distinct objectives. A model may achieve good rollouts by using D(q)D(q) as a “stabilizer” that absorbs errors from V or M, rather than learning the true physical damping. More broadly, without additional structure the inverse problem can be non-identifiable (gauge freedom), even in conservative systems; Appendix H.8 provides a controlled illustration on the double pendulum. We therefore separate forecasting and identifiability in both the model design and the evaluation (Sec. 4.1.2), and empirically observe that stronger physical structure improves identifiability (Table 5). Our contributions are: 1. Unified architecture: A single port-Hamiltonian template spanning three knowledge regimes (KNOWN, PARTIAL, UNKNOWN) by swapping component instantiations for V, M, D. 2. Householder-style parameterizations: Low-rank outer-product expansions that guarantee D(q)⪰0D(q) 0 and M(q)≻0M(q) 0 by construction with efficient structured primitives: applying D(q)D(q) to a vector costs O(nr)O(nr), and computing v=M(q)−1pv=M(q)^-1p via Woodbury is O(nr2+r3)O(nr^2+r^3), where n is the number of degrees of freedom and r is the number of rank-1 terms (Appendix H.4). 3. Damping bounds for identifiability. We show that unconstrained dissipation can absorb model mismatch and obscure physical parameter recovery. PHAST therefore bounds the effective damping strength, improving identifiability in grey-box settings while preserving stability (Table 5; Appendix H) 4. Structure-preserving integration: Strang splitting that preserves the continuous-time dissipation structure to O(Δt2)O( t^2). 5. Two-axis evaluation: Explicit separation of forecasting metrics (rollout MSE) from identifiability metrics (damping R2R^2), revealing trade-offs invisible to single-axis evaluation. Across thirteen q-only benchmarks—including non-mechanical systems (RLC circuit, Lennard–Jones cluster, heat exchange, N-body gravity, predator–prey)—PHAST achieves the best long-horizon open-loop stability among competitive baselines (Tables 17–24); on Windy Pendulum, PHAST (KNOWN) recovers the true position-dependent damping with R2=0.996R^2=0.996 (Table 5). These results suggest that, on these q-only benchmarks, structure is a stronger driver of long-horizon stability than additional capacity. Code to reproduce our experiments will be released. 2 Related Work Learning dynamics from data spans energy-based methods, structure-preserving networks, and efficient sequence models; we position PHAST relative to each. Energy-based learning. Hamiltonian Neural Networks (HNNs) (Greydanus et al., 2019) learn Hθ(x)H_θ(x) and set x˙=J∇Hθ x=J∇ H_θ, while Lagrangian Neural Networks (LNNs) (Cranmer et al., 2020) learn Lθ(q,q˙)L_θ(q, q) and derive Euler–Lagrange dynamics. Both are conservative by construction (R=0R=0), and LNNs additionally require Hessian computations. Dissipative extensions. Dissipative SymODEN (Zhong et al., 2020) and DHNN (Sosanya & Greydanus, 2022) extend HNN/LNN with damping terms. However, dissipation is often unconstrained and can entangle with errors in V or M, producing reasonable rollouts but poor physical recovery. PHAST enforces D(q)⪰0D(q) 0 by construction, optionally bounds its spectrum, and evaluates forecasting vs. identifiability separately. Structure-preserving networks. SympNets (Jin et al., 2020) and Volume-Preserving Transformers (Brantner et al., 2024) enforce symplectic/volume-preserving updates aligned with conservative Hamiltonian flows. For dissipative systems, strict volume preservation cannot represent generic phase-space contraction and provides no passivity guarantee. Port-Hamiltonian learning. Prior work on learning port-Hamiltonian systems (Desai et al., 2021; Eidnes & Riemer-Sørensen, 2023) largely assumes known structure, while concurrent work (Anonymous, 2026) learns generalized Hamiltonian dynamics from phase-space data with probabilistic constraints. PHAST targets partially known/unknown structure and q-only learning while keeping the J–R decomposition explicit. Hamiltonian methods for optimal control. The PHAST architecture builds upon a sustained research trajectory investigating the integration of geometric, Hamiltonian, and stochastic structures into learning-based dynamics and control. Early work addressed the fundamental brittleness of Physics-Informed Neural Networks (PINNs) in the presence of noise, proposing Gaussian Process smoothing to recover robust learning when physical priors fail (Bajaj et al., 2023). Complementing this, the Motion Code framework established methods for disentangling complex, noisy time-series data into underlying stochastic processes (Bajaj & Nguyen, 2024). To bridge deep learning with optimal control, the NeuralPMP framework was introduced, leveraging the Pontryagin Maximum Principle to learn reduced Hamiltonian dynamics for continuous control tasks (Bajaj & Nguyen, 2021), a methodology subsequently adapted to navigate the complex energy landscapes of molecular dynamics (Bajaj et al., 2024). This paper’s focus on Hamiltonian structures was generalized to learn diverse dynamical classes—including conservative, dissipative, and port-Hamiltonian systems—from noisy trajectories by enforcing stability constraints (McLennan et al., 2025). Simultaneously, the theoretical underpinnings of Reinforcement Learning were reformulated through a differential dual perspective in Differential Policy Optimization (dfPO), enabling pointwise policy refinement (Nguyen & Bajaj, 2025a), which was further extended to stochastic settings using rough path theory (Nguyen & Bajaj, 2025b). Most recently, these geometric principles were applied to simultaneous navigation and mapping (SNAM), demonstrating how learned Hamiltonian energy landscapes can guide robotic navigation in unknown environments (Ellendula et al., 2025). Collectively, these papers motivate PHAST’s unified approach to modeling dissipative systems across varying regimes of prior knowledge. Efficient sequence models. State-space models such as S4/S5 and Mamba (Gu et al., 2022; Smith et al., 2023; Gu & Dao, 2023) provide linear-time recurrence; oscillatory variants like LinOSS and D-LinOSS (Hasani et al., 2024; Boyer et al., 2025) learn stable oscillation and dissipation time scales. Unlike PHAST, these models typically do not expose explicit energies or configuration-dependent damping fields, and they provide no port-Hamiltonian passivity guarantee. The gap. Table 2 summarizes the landscape and the key properties we target. Table 2: Comparison of dynamics learning approaches. PHAST combines dissipative modeling (R≠0R≠ 0), passivity-aware structure, spectral control, and efficient primitives for long-horizon rollouts. Passivity indicates a by-construction storage-function inequality. Method R≠0R≠ 0 Passivity Spectral Efficient HNN ✗ ✓ ✗ ✓ LNN ✗ ✓ ✗ ✗ DHNN / Diss. SymODEN ✓ ✗ ✗ ✓ SympNets / VPT ✗ ✗ ✗ ✓ SSMs (S5, LinOSS, D-LinOSS) — ✗ ✗ ✓ PHAST (ours) ✓ ✓ ✓ ✓ 3 Method: PHAST for Dissipative Dynamics This section describes how PHAST enforces the structural decomposition introduced above (Table 1, Fig. 1) and preserves it under time discretization. We begin with the observation model, data contract, and passivity guarantees (Sec. 3.1), followed by the parameterizations that enforce definiteness (Sec. 3.2), the structure-preserving integrator (Sec. 3.3), the training objective (Sec. 3.4), and port-based control (Sec. 3.5). 3.1 Problem Formulation: Partial Observability (q-only) As described in Section 1, PHAST operates in the q-only setting where only generalized positions qtq_t are observed and momenta ptp_t are latent. We now detail the three pipeline stages—velocity observer, canonicalizer, and port-Hamiltonian core—that map a short position history to long-horizon phase-space rollouts. Formally, we observe yt=qt,xt=(qt,pt),y_t=q_t, x_t=(q_t,p_t), (5) where ptp_t is the (generally unobserved) conjugate momentum in the underlying phase state. Velocity observer. We first compute a simple finite-difference estimate q˙tfd:=0,t=0,(qt−qt−1)/Δt,t≥1, q_t^fd:= cases0,&t=0,\\ (q_t-q_t-1)/ t,&t≥ 1, cases (6) and then apply a causal Temporal Convolutional Network (TCN) observer oϕo_φ that predicts an additive correction: q˙^t=q˙tfd+δt,δ0:t=oϕ([q0:t,q˙0:tfd]). q_t= q_t^fd+ _t, _0:t=o_φ([q_0:t, q^fd_0:t]). (7) Canonicalization. The estimated velocity is mapped to a phase state via a canonicalizer Cψ:(qt,q˙^t)↦(qt,p^t)C_ψ:(q_t, q_t) (q_t, p_t). For mechanical systems, the physical relation is p=M(q)q˙p=M(q) q. When M is known, we use the mass-consistent mapping p^=M(q)q˙ p=M(q) q. In q-only benchmarks with constant mass, we use the identity canonicalizer p^=q˙ p= q and treat p p as a velocity-like latent variable that yields an approximately Markov phase state. Rollout modes. Given a burn-in window q0:K−1q_0:K-1, PHAST supports two rollout modes. In takeover mode, an initial phase state is inferred at the end of the context and then integrated open-loop: x^K−1 x_K-1 =Cψ(qK−1,q˙^K−1), =C_ψ(q_K-1, q_K-1), x^t+1 x_t+1 =ΦΔt(x^t),q^t=Πq(x^t). = _ t( x_t), q_t= _q( x_t). (8) In autoregressive (self-conditioned) mode, the observer is re-applied to the generated position prefix at each step. For cross-model comparisons, we report autoregressive rollouts; takeover rollouts are additionally used as a diagnostic of the learned dynamics. When M(q)M(q) is available from physics (KNOWN regime), we use the mass-consistent canonicalizer p^=M(q)q˙ p=M(q) q; combined with a nonseparable Hamiltonian integrator, this yields substantial gains over the constant-mass approximation (e.g., ∼3× 3× lower H=100H=100 rollout error on Windy Cart-Pole; Appendix H.2, Table 26). Unforced dynamics (u≡0u≡ 0). All forecasting and filtering experiments in this work consider the autonomous system x˙=(J−R)∇H(x) x=(J-R)∇ H(x) (Eq. (1) with u=0u=0), which admits the energy balance dHdt=−∇H(x)⊤R∇H(x)≤0. dHdt=-∇ H(x) R∇ H(x)≤ 0. (9) Equation (9) is the passivity certificate: the learned Hamiltonian H serves as a storage function that is non-increasing along trajectories, ensuring bounded energy and stable long-horizon rollouts. Forced dynamics and ports (u≠0u≠ 0). The input u enters through a port matrix G, with conjugate output yport=G⊤∇H(x)y^port=G ∇ H(x). The supplied power is u⊤yportu y^port, yielding the balance dHdt=−∇H⊤R∇H+u⊤yport dHdt=-∇ H R∇ H+u y^port, which reduces to Eq. (9) when u=0u=0. Section 3.5 develops the forced form into a closed-loop Energy–Casimir controller. 3.2 Householder-Style Low-Rank Parameterizations To enforce definiteness constraints and keep inference fast, PHAST uses low-rank outer-product expansions (“Householder-style” additive updates; not Householder reflections). Damping (all regimes). We parameterize the damping matrix as D(q) D(q) =d0I+∑i=1rβi(q)ki(q)ki(q)⊤, =d_0I+ _i=1^r _i(q)\,k_i(q)k_i(q) , d0≥0,βi(q)≥0,‖ki(q)‖=1for i=1,…,r, d_0≥ 0,\ _i(q)≥ 0,\ \|k_i(q)\|=1\ for i=1,…,r, (10) where n is the number of degrees of freedom (so q∈ℝnq ^n) and r is the number of rank-1 terms. This guarantees D(q)⪰0D(q) 0 and admits O(nr)O(nr) matrix-vector products. Bounding damping strength (identifiability + stability). Optionally, we bound the total strength ∑i=1rβi(q)≤β¯⇒λmax(D(q))≤d0+β¯, _i=1^r _i(q)≤ β _ (D(q))≤ d_0+ β, (11) which (i) prevents D from acting as an “error sink” when V or M are misspecified, and (i) improves numerical stability of the explicit damping half-step (Algorithm 1, lines 2–3 and 10–11) when Δtλmax(D) t\, _ (D) is large. In our grey-box Windy benchmarks, we set β¯ β to a physics-calibrated magnitude (e.g., the known damping variation Δd d) in the PARTIAL regime (Appendix E). We enforce Eq. (11) by reparameterizing the strengths as βi(q)=βmaxσ(⋅) _i(q)= _ \,σ(·) with βmax=β¯/r _ = β/r, so that each term is nonnegative and uniformly bounded. Appendix Fig. 23 provides a schematic intuition for Eq. (11). Mass. PHAST supports a general configuration-dependent mass M(q)≻0M(q) 0. In our main experiments we use a constant-mass approximation M(q)≈M(q)≈ M for efficiency, enabling a fast leapfrog integrator core; the nonseparable M(q)M(q) variant uses an implicit-midpoint integrator and is studied in Appendix H.2. When learning M (UNKNOWN regime), we use an SPD diagonal-plus-low-rank parameterization and compute M−1pM^-1p via Woodbury; see Appendix Eq. 79 and Alg. 2. 3.3 Structure-Preserving Integration For dissipative systems we integrate with Strang splitting (Strang, 1968; Hairer et al., 2006): ΦΔt=ΦDΔt/2∘ΦHΔt∘ΦDΔt/2, _ t= _D t/2 _H t _D t/2, (12) where ΦHΔt _H t is a conservative (symplectic) update and ΦDΔt/2 _D t/2 is a half-step of damping. In our main experiments we use the leapfrog (Störmer–Verlet) conservative core; we additionally report ablations with implicit midpoint and Cayley–Neumann in settings where additional stability is helpful. Stacking and timestep. The PHAST transition ΦΔt _ t can optionally compose L integration substeps to produce a single predicted transition xt↦xt+1x_t x_t+1. All substeps reuse the same learned energy components (V, M, D) and differ only in the intermediate states they act on (and, if enabled, per-substep timesteps). In our experiments we set L=1L=1 for interpretability and efficiency, and initialize the internal timestep to the environment timestep Δt t; unless otherwise noted, this timestep is learnable (Appendix H reports L and timestep ablations). xt=(q,p)x_t=(q,p)V(q)V(q)M≻0M 0D(q)⪰0D(q) 0potentialmassdampingH=V+12p⊤M−1pH=V+ 12p M^-1pR=diag(0,D)R=diag(0,D)x˙=(J−R)∇H x=(J-R)∇ Hstructure-preserving ΦΔt _ txt+1x_t+1 Unroll for trajectory (L substeps per environment step):x0x_0substepsx1x_1substepsx2x_2substeps⋯·sH1≲H0H_1 H_0H2≲H1H_2 H_1H3≲H2H_3 H_2 Figure 3: PHAST core transition computation graph. Components V, M, and D assemble into Hamiltonian H and dissipation R; in the unforced case (u=0u=0), the port-Hamiltonian form implies dHdt≤0 dHdt≤ 0 in continuous time. PHAST advances the phase state using a structure-preserving discrete-time map ΦΔt _ t (Strang splitting) and optionally composes L integration substeps per environment step (Eq. (78)). The Hk+1≲HkH_k+1 H_k annotations are schematic; discrete-time energy monotonicity depends on timestep and integrator (Appendix A.2). Algorithm 1 PHAST step ΦΔt(q,p)→(q+,p+) _ t(q,p)→(q^+,p^+) (Strang splitting with leapfrog core; constant mass) 0: (q,p)∈ℝ2n(q,p) ^2n, time step Δt t 1: Half dissipation step 2: v←M−1pv← M^-1p 3: p←p−Δt2D(q)vp← p- t2\,D(q)\,v 4: Conservative step (leapfrog / Störmer–Verlet) 5: p1/2←p−Δt2∇V(q)p_1/2← p- t2\,∇ V(q) 6: v1/2←M−1p1/2v_1/2← M^-1p_1/2 7: q+←q+Δtv1/2q^+← q+ t\,v_1/2 8: p+←p1/2−Δt2∇V(q+)p^+← p_1/2- t2\,∇ V(q^+) 9: Half dissipation step 10: v+←M−1p+v^+← M^-1p^+ 11: p+←p+−Δt2D(q+)v+p^+← p^+- t2\,D(q^+)\,v^+ 12: return (q+,p+)(q^+,p^+) Discrete-time caveat. Continuous-time passivity in Eq. (9) does not automatically imply monotone energy decay at finite step size. With the explicit damping half-step, discrete energy monotonicity can fail if the system is stiff; bounding Eq. (11) (or using an implicit damping update) mitigates this in practice. 3.4 Training Objective PHAST supports end-to-end training with a weighted multi-term objective: ℒ=λdataℒdata+λpassℒpass+λenergyℒenergy+λrollℒroll.L= _dataL_data+ _passL_pass+ _energyL_energy+ _rollL_roll. (13) Loss breakdown. Let y0:T−1y_0:T-1 denote the observed sequence (in q-only, yt=qty_t=q_t), where T is the sequence length. Let err(y^,y)err( y,y) denote the appropriate per-step error on the manifold (Appendix F). For full-state training (when x=(q,p)x=(q,p) is available), we use teacher-forced one-step predictions x^t+1=ΦΔt(xt) x_t+1= _ t(x_t) and define Ht=H(xt)H_t=H(x_t) and H^t+1=H(x^t+1) H_t+1=H( x_t+1): ℒdata _data =1T−1∑t=0T−2err(y^t+1,yt+1), = 1T-1 _t=0^T-2err( y_t+1,y_t+1), (14) ℒpass _pass =1T−1∑t=0T−2max(0,H^t+1−Ht), = 1T-1 _t=0^T-2 (0, H_t+1-H_t), (15) ℒenergy _energy =1T−1∑t=0T−2|H^t+1−HtΔt+vt⊤D(qt)vt|, = 1T-1 _t=0^T-2 | H_t+1-H_t t+v_t D(q_t)v_t |, (16) vt:=M−1pt, v_t:=M^-1p_t, (17) ℒroll _roll =t0[1Hroll∑h=1Hrollerr(y~t0+h,yt0+h)]. =E_t_0 [ 1H_roll _h=1^H_rollerr( y_t_0+h,y_t_0+h) ]. (18) Here x~0=xt0 x_0=x_t_0, x~h+1=ΦΔt(x~h) x_h+1= _ t( x_h), and t0t_0 is sampled uniformly from valid start indices. In the q-only setting, energy diagnostics can be formed using finite-difference velocity estimates as in Appendix F. We always train with ℒdataL_data for next-step prediction; ℒpassL_pass penalizes discrete energy increases and ℒenergyL_energy matches a dissipation budget (when labels are available) or a label-free residual under Strang splitting. Optionally, ℒrollL_roll trains the model in the same open-loop mode used at evaluation, reducing teacher-forcing train–test mismatch. Unless otherwise noted, our main benchmark tables use λpass=λenergy=λroll=0 _pass= _energy= _roll=0 and report long-horizon open-loop stability and energy/damping diagnostics at evaluation time. The complete training pseudocode is given in Algorithm 3. 3.5 Energy-Casimir Control (forced dynamics mode) Why dynamic control. Static (memoryless) feedback can inject damping but cannot reshape the energy landscape to stabilize a non-natural equilibrium q⋆q . A dynamic controller that stores energy is needed: a port-Hamiltonian (pH) controller exchanges energy with the plant through power ports, adding a shaped potential that moves the minimum of the closed-loop energy to q⋆q (Fig. 7). Plant and controller as two pH systems. Let the plant be a controlled pH system x˙ x =(Jp(x)−Rp(x))∇Hp(x)+Gp(x)up, =(J_p(x)-R_p(x))∇ H_p(x)+G_p(x)\,u_p, (19) yp y_p =Gp(x)⊤∇Hp(x), =G_p(x) ∇ H_p(x), (20) with energy balance H˙p=−(∇Hp)⊤Rp(∇Hp)+yp⊤up H_p=-(∇ H_p) R_p(∇ H_p)+y_p u_p. Here upu_p is the control input (torque-like) and ypy_p is the conjugate port output (velocity-like); their product yp⊤upy_p u_p is the power supplied to the plant. The controller has a virtual state ξ with shaped storage Hc(ξ):=kc2(ξ−q⋆)2,H_c(ξ):= k_c2\,(ξ-q )^2, (21) acting as a virtual spring centered at the target q⋆q . Its port variables are ξ˙=uc,yc=∇Hc(ξ)=kc(ξ−q⋆). ξ=u_c, y_c=∇ H_c(ξ)=k_c(ξ-q ). (22) Power-preserving interconnection. We couple plant and controller through the port assignment up=−yc+v,uc=yp,u_p=-y_c+v, u_c=y_p, (23) where v is an auxiliary channel (set to zero for now). When v=0v=0 the cross-power cancels: yp⊤up+yc⊤uc=−yp⊤yc+yc⊤yp=0y_p u_p+y_c u_c=-y_p y_c+y_c y_p=0, so the coupling transfers energy between plant and controller but creates none. Casimir invariant. Define (q,ξ):=q−ξC(q,ξ):=q-ξ. Because uc=ypu_c=y_p and, for collocated mechanical systems, yp=q˙y_p= q, we obtain ˙=q˙−ξ˙=yp−uc=0 C= q- ξ=y_p-u_c=0. If ξ(0)=q(0)ξ(0)=q(0), then ξ(t)=q(t)ξ(t)=q(t) for all t, so Hc(ξ)H_c(ξ) acts as a shaped potential on the plant configuration: the controller’s virtual spring tracks the physical position forever. C is a Casimir function—a quantity conserved by the interconnection structure (i.e. ∇∈kerJcl⊤ ∈ J_cl ), independently of the energy functions. Closed-loop energy decrease. The closed-loop storage Hcl:=Hp+HcH_cl:=H_p+H_c satisfies H˙cl=−(∇Hp)⊤Rp(∇Hp)+yp⊤v. H_cl=-(∇ H_p) R_p(∇ H_p)+y_p v. (24) With damping injection v=−dinjypv=-d_inj\,y_p (dinj>0d_inj>0), this yields H˙cl≤−dinj‖yp‖2≤0 H_cl≤-d_inj\|y_p\|^2≤ 0, guaranteeing asymptotic convergence to the minimum of HclH_cl at q⋆q . Forced dynamics mode. In the experiments of Appendix 5 we additionally inject an exogenous signal vextv_ext (for excitation or tracking), setting v(t)=vext(t)−dinjy^(t),v(t)=v_ext(t)-d_inj\, y(t), (25) where y y is the (possibly estimated) port output. Substituting (23)–(25) into the power balance gives yp⊤up+yc⊤uc=yp⊤vext−dinjyp⊤y^.y_p u_p+y_c u_c=y_p v_ext-d_inj\,y_p y. (26) In the ideal case y^≡yp y≡ y_p, the closed-loop storage satisfies the forced passivity bound H˙cl≤yp⊤vext−dinj‖yp‖2, H_cl≤ y_p v_ext-d_inj\|y_p\|^2, (27) i.e. the closed loop is passive from vextv_ext to ypy_p with strict dissipation. The Casimir invariance is compatible with forced mode: vextv_ext enters only through upu_p and does not alter the constraint uc=ypu_c=y_p. Q-only limitation. In q-only settings ypy_p is not observed directly and must be replaced by an estimate y^≈yp y≈ y_p. The dissipative term then becomes −dinjyp⊤y^-d_inj\,y_p y, which can lose definiteness when y^≠yp y≠ y_p: estimation error may inject or dissipate energy. This is the central partial-observability bottleneck for port-based control. Discrete-time implementation. Let qtmeasq_t^meas denote q-only measurements and y^t y_t the online port estimate. The controller is implemented as ut u_t =−yc(ξt)+vext,t−dinjy^t, =-y_c( _t)+v_ext,t-d_inj\, y_t, (28) ξt+1 _t+1 =ξt+Δt(y^t+kξ(qtmeas−ξt)), = _t+ t ( y_t+k_ξ(q_t^meas- _t) ), (29) where kξ≥0k_ξ≥ 0 is an optional predictor–corrector gain that corrects ξ drift under partial observability: when y^=yp y=y_p and kξ=0k_ξ=0, the update reduces to a forward-Euler discretization of ξ˙=yp ξ=y_p. All variants share the same (kc,dinj,kξ)(k_c,d_inj,k_ξ) and forcing signal vext,tv_ext,t; they differ only in how y^t y_t is produced online. Port-estimation variants. We compare the same family of estimators as in the unforced setting (oracle/full-state, finite differences, fixed-lag MAP smoothing, learned observer, etc.) and evaluate stability and control-effort degradation as a function of port-estimation error (Appendix 5). 4 Experiments We evaluate PHAST on q-only systems spanning conservative and dissipative dynamics across six physical domains. Our evaluation separates two regimes with distinct goals: open-loop forecasting (autonomous rollouts) to assess long-horizon stability and physical identifiability, and closed-loop control (feedback) to assess whether the learned port-Hamiltonian structure remains useful under stabilization and to isolate the role of online port estimation under partial observability. 4.1 Open-loop forecasting (q-only) 4.1.1 Experimental setup Environments. We evaluate on thirteen q-only benchmarks spanning conservative and dissipative dynamics across six physical domains. Mechanical systems (eight benchmarks): single pendulum (conservative, constant damping, and position-dependent “windy” damping), a windy Cart-Pole on ℝ×1R×S^1, harmonic oscillator (conservative and damped), and double pendulum (conservative and damped). Non-mechanical systems (five benchmarks covering all three knowledge regimes in Table 1): series RLC circuit (KNOWN), Lennard–Jones 3-particle cluster (KNOWN), coupled heat exchange (KNOWN), 3-body gravitational system (PARTIAL), and Lotka–Volterra predator–prey (UNKNOWN). In all environments only configurations qtq_t are observed; the simulator evolves trajectories from random initial conditions (q0,p0)(q_0,p_0). All models are trained as next-step predictors that map a short history of observed configurations to q^t+1 q_t+1; PHAST internally uses the same history to infer a latent phase state via the FD+TCN observer (Sec. 3.1), while sequence baselines must implicitly learn any required state estimation. All forecasting benchmarks are unforced (u=0u=0 in Eq. (1)). For windy settings we use a position-dependent diagonal damping of the form d(q)=d0+Δd|sinq|d(q)=d_0+ d\,| q| (e.g., d0=0.3d_0=0.3, Δd=0.5 d=0.5), applied to the relevant angular coordinate. Some benchmarks (Cart-Pole and Double Pendulum) have configuration-dependent inertia M(q)M(q) in the simulator; unless otherwise stated, we use a separable constant-mass approximation for efficiency (Sec. 1) and study the nonseparable M(q)M(q) variant in Appendix H.2. Appendix B provides detailed environment descriptions and Hamiltonians. Data protocol. Trajectories have length T=200T=200. We use Δt=0.05 t=0.05 for single-pendulum environments, Δt=0.02 t=0.02 for Cart-Pole and Oscillator, and Δt=0.01 t=0.01 for Double Pendulum. For PHAST, we initialize the internal integrator step size to the environment timestep Δt t and use a single integration substep per observation step (L=1L=1); unless otherwise stated, the internal timestep is learnable. Dataset sizes are Ntrain/Nval/Ntest=1000/200/200N_train/N_val/N_test=1000/200/200. Unless otherwise noted, tables report mean ± std over 5 random model seeds using a fixed dataset shared across models, and all runs are executed on CPU for reproducibility. All models are trained for 50 epochs; complete hyperparameters are provided in Appendix E. Training and model selection. All methods use the same optimizer and schedule: AdamW (weight decay 10−510^-5) with a cosine learning-rate schedule. We select the model parameters with lowest validation MSE (evaluated every 10 epochs) and report test metrics for that selection. Baselines. We compare against GRU (Cho et al., 2014), S5 (Smith et al., 2023), LinOSS (Hasani et al., 2024), D-LinOSS (Boyer et al., 2025), a Transformer (Vaswani et al., 2017), and a Volume-Preserving Transformer (VPT) (Brantner et al., 2024). We report trainable parameter counts for all methods; parameter matching is treated as an experimental control rather than an assumption. All baselines are trained and evaluated as causal seq2seq predictors that map a q-only history to the next configuration (q0:t↦q^t+1q_0:t q_t+1); they do not use an explicit observer/canonicalizer pipeline. 4.1.2 Evaluation protocol: two axes Rollout MSE → (damping R2R^2↑ )IdealPARTIALKNOWNUNKNOWND-LinOSSS5GRUmore structure Figure 4: Two-axis evaluation (open-loop, conceptual). Forecasting accuracy (low rollout MSE) and physical identifiability (high damping R2R^2) are distinct objectives; model/regularizer choices induce trade-offs. Forecasting accuracy and physical parameter recovery are distinct objectives. A model may achieve good rollouts by using D(q)D(q) as a “stabilizer” rather than learning true physical damping. Conversely, recovering true D(q)D(q) does not guarantee good open-loop rollouts if V(q)V(q) or M have errors. We therefore evaluate open-loop performance along two axes (Fig. 4): Axis 1: Forecasting stability. For angular coordinates we report one-step wrapped-angle MSE and long-horizon open-loop rollout error at horizon H=100H=100 with burn-in context K=10K=10. For Euclidean q-only environments we report standard MSE and rollout MSE at H=100H=100. For mixed manifolds (Cart-Pole) we report a mixed-manifold rollout MSE at H=100H=100 that averages translation MSE and wrapped-angle MSE. All models are evaluated in open-loop (autoregressive) mode; PHAST additionally supports takeover rollouts (burn-in then integrate; Appendix C.3) as a diagnostic of the learned dynamics. Axis 2: Identifiability. When a model exposes an explicit damping field D(q)D(q), we report damping recovery via R2R^2 and MAE. We also report an energy-consistency diagnostic on open-loop rollouts (energy-budget residual at H=100H=100), computed from finite-difference velocity estimates and the benchmark’s simulator damping law Denv(⋅)D_env(·) for expected dissipation (Appendix F). Precise metric definitions (including rollout protocol and the discrete energy-budget residual) are provided in Appendix F. Regime-specific expectations. In KNOWN/PARTIAL regimes, physics-calibrated damping bounds (e.g., ∑iβi(q)≤β¯≈Δd _i _i(q)≤ β≈ d) can improve identifiability without destabilizing rollouts (Table 6; Appendix Table 31). In UNKNOWN, damping bounds can reveal a forecasting–identifiability trade-off unless additional anchors break gauge freedoms (Appendix Table 32). Open-loop error sources & diagnostics. We report rollout error (stability), damping R2R^2 (identifiability), energy-budget residual, and discrete-time passivity violations. Appendix F.1 provides the q-only computation graph and diagnostic interpretations (autoregressive vs. takeover, step-size/substepping, damping bounds, and time-scale ambiguity). 4.1.3 Main open-loop results Main qualitative results. We present qualitative results across three benchmarks that probe complementary aspects of dissipative dynamics learning (Fig. 5). The Windy Pendulum (1 DOF, 1S^1) is the simplest: constant mass, a single position-dependent damping coefficient d(θ)d(θ), and regular oscillatory trajectories. It tests whether a model can learn basic dissipative Hamiltonian dynamics from positions alone. The Windy Double Pendulum (2 DOF, 2T^2) escalates difficulty sharply: the dynamics are chaotic, the mass matrix M(q)M(q) couples the two links through cos(θ1−θ2) ( _1- _2), and each joint has independent wind damping. Small prediction errors grow exponentially, making this the hardest forecasting challenge in our suite. The Windy Cart-Pole (2 DOF, ℝ×1R×S^1) represents a different axis of difficulty: a mixed manifold (translation × rotation), configuration-dependent inertia M(q)M(q), and two qualitatively different damping mechanisms — constant viscous cart friction and position-dependent angular wind damping on the pole. Full forecasting tables across all thirteen q-only benchmarks are reported in Appendix G. Table 4 provides a compact summary across the full benchmark suite (see also Fig. 2 in Sec. 1 for a visual overview). Because the Double Pendulum and Cart-Pole have nonseparable Hamiltonians (M(q)≠constM(q) ), we use the implicit-midpoint integrator for all PHAST curves in the multi-environment figures. Appendix H.2 further discusses the nonseparable mass variant. Figure 5: Open-loop rollouts and phase-space portraits across three dissipative benchmarks (q-only). Top row: a single test trajectory is teacher-forced through a short burn-in window (grey region, vertical dashed line), then predicted open-loop for H=100H=100 steps. Column 1 (Windy Pendulum, θ): the simplest system — a single angle with position-dependent damping. PHAST (KNOWN/PARTIAL) tracks the decaying oscillation almost exactly; all three baselines diverge within a few periods. Columns 2–3 (Windy Double Pendulum, θ1 _1 and θ2 _2): a chaotic 2-DOF system with coupled configuration-dependent inertia M(q)M(q). This is the hardest benchmark: small errors grow exponentially, yet PHAST maintains trajectory coherence on both joints far longer than baselines. Showing both angles reveals that the model captures inter-joint coupling, not just marginal statistics. Column 4 (Windy Cart-Pole, θ): a mixed-manifold system (ℝ×1R×S^1) with two qualitatively different damping mechanisms (constant cart friction ++ angular wind damping). PHAST (PARTIAL) achieves the tightest tracking here. Bottom row: canonical phase-space portraits (θ,p)(θ,p) with p=M(q)q˙p=M(q) q during the open-loop segment. Momentum is latent in the q-only setting: for PHAST, q˙ q comes from the learned FD+TCN observer (Sec. 3); for baselines, q˙ q is approximated by finite differences. The closed orbits visible for PHAST confirm that the model has learned physically consistent Hamiltonian structure, whereas baseline phase portraits collapse or spiral outward. Table 3: Q-only open-loop forecasting (Cart-Pole). Mean ± std of mixed-manifold rollout MSE at horizon H=100H=100 over 5 seeds. The mixed metric averages translation MSE and wrapped-angle MSE (Appendix F). Lower is better. Model Params Cart-Pole (windy) PHAST (ours) PHAST (KNOWN) 3,589 0.063±0.0190.063± 0.019 PHAST (PARTIAL) 14,283 0.083±0.0220.083± 0.022 PHAST (UNKNOWN) 14,290 0.109±0.0220.109± 0.022 S5 (best baseline) 17,218 0.431±0.0770.431± 0.077 Table 3 shows that PHAST achieves strong long-horizon stability on Windy Cart-Pole. Figure 6 shows the learned damping profiles and energy traces across all three systems, illustrating the forecasting–identifiability trade-off: Table 5 highlights the two-axis nature of the problem: in the Windy setting, the KNOWN regime recovers a physically meaningful damping field, PARTIAL becomes identifiable once we impose a physics-calibrated magnitude bound β¯ β on D(q)D(q), and UNKNOWN still exhibits a forecasting–identifiability trade-off without additional anchors (Appendix F). Appendix H makes this trade-off explicit: Table 6 compares bounded vs. unbounded damping in the PARTIAL regime, showing that without the bound the model can achieve accurate rollouts while learning a physically meaningless damping field (RD2≪0R^2_D 0). Takeover rollouts (burn-in then integrate) can improve with better identifiability, but remain sensitive to single-shot momentum inference at the context boundary; we therefore treat takeover as a diagnostic rather than the primary q-only comparison metric. Table 4: Suite summary across thirteen q-only benchmarks. For each benchmark, we report the best PHAST regime and the best baseline (mean ± std over 5 seeds). The first eight rows use rollout MSE at horizon H=100H=100; the last five report next-step MSE (Tables 20–24). Gains are per-row ratios and are not directly comparable across environments. Benchmark Best PHAST Best baseline Gain Mechanical systems (rollout MSE at H=100H=100) Pendulum (cons) 0.680±0.0430.680± 0.043 (PARTIAL) 2.320±0.2242.320± 0.224 (Transformer) 3.4×3.4× Pendulum (damped) 0.017±0.0050.017± 0.005 (KNOWN) 0.450±0.2410.450± 0.241 (D-LinOSS) 26.5×26.5× Pendulum (windy) 0.092±0.0140.092± 0.014 (PARTIAL) 0.435±0.2390.435± 0.239 (D-LinOSS) 4.7×4.7× Cart-Pole (windy) 0.063±0.0190.063± 0.019 (KNOWN) 0.431±0.0770.431± 0.077 (S5) 6.8×6.8× Oscillator (cons) 0.0010±0.00020.0010± 0.0002 (KNOWN/PARTIAL) 1.087±0.2991.087± 0.299 (Transformer) 1.1×1031.1× 10^3 Oscillator (damped) 0.0011±0.00030.0011± 0.0003 (PARTIAL) 0.926±0.2540.926± 0.254 (Transformer) 8.4×1028.4× 10^2 Double pendulum (cons) 0.402±0.0470.402± 0.047 (PARTIAL) 0.618±0.0280.618± 0.028 (S5) 1.5×1.5× Double pendulum (damped) 0.320±0.0320.320± 0.032 (PARTIAL) 0.630±0.0310.630± 0.031 (S5) 2.0×2.0× Non-mechanical systems (next-step MSE) RLC circuit (KNOWN) 2.63×−2.63× 10^-5 (UNKNOWN) 4.81×10−44.81× 10^-4 (Transformer) 18×18× LJ-3 cluster (KNOWN) 4.59×−4.59× 10^-10 (PARTIAL) 2.05×10−42.05× 10^-4 (S5) 4.5×1054.5× 10^5 Heat exchange (KNOWN) 2.42×−2.42× 10^-6 (KNOWN) 4.46×10−44.46× 10^-4 (LinOSS) 1.8×1021.8× 10^2 N-body gravity (PARTIAL) 4.27×−4.27× 10^-8 (PARTIAL) 1.83×10−31.83× 10^-3 (Transformer) 4.3×1044.3× 10^4 Predator–prey (UNKNOWN) 0.01990.0199 (UNKNOWN) 0.1790.179 (Transformer) 9.0×9.0× Table 5: Windy Pendulum (q-only) identifiability and energy consistency. Mean ± std over 5 seeds (same setup as Table 17). Baselines do not expose a damping field, so damping recovery metrics are not applicable. Model WrapMSEθroll(H=100)WrapMSE_θ^roll(H=100) ↓ RD2R^2_D ↑ MAEDMAE_D ↓ EbudResroll(H=100)EbudRes^roll(H=100) ↓ Params PHAST (KNOWN) 0.106±0.0200.106± 0.020 0.996±0.0040.996± 0.004 0.007±0.0040.007± 0.004 1.56±0.021.56± 0.02 3,364 PHAST (PARTIAL) 0.092±0.0140.092± 0.014 0.654±0.0770.654± 0.077 0.063±0.0110.063± 0.011 1.50±0.041.50± 0.04 13,736 PHAST (UNKNOWN) 0.298±0.0480.298± 0.048 −96.546±15.897-96.546± 15.897 1.343±0.1151.343± 0.115 2.58±0.192.58± 0.19 13,738 GRU 1.796±0.6251.796± 0.625 — — 12.34±6.5812.34± 6.58 37,889 S5 0.600±0.0470.600± 0.047 — — 11.71±2.9111.71± 2.91 17,089 LinOSS 1.458±0.3241.458± 0.324 — — 56.14±37.4156.14± 37.41 17,089 D-LinOSS 0.435±0.2390.435± 0.239 — — 6.86±1.146.86± 1.14 33,793 Transformer 0.824±0.1340.824± 0.134 — — 7.71±0.637.71± 0.63 100,161 VPT 2.218±0.1352.218± 0.135 — — 40.48±8.0440.48± 8.04 16,833 Table 6: PARTIAL regime: effect of bounding total damping strength (Windy Pendulum, q-only). Mean ± std over 5 seeds on CPU (same data protocol as Sec. 4.1.1). Without the bound, damping becomes an “error sink” and is non-identifiable (RD2≪0R^2_D 0) even when rollouts are accurate; bounding ∑i=1rβi(q) _i=1^r _i(q) recovers a physically meaningful damping field. Damping cap WrapMSEθ(100) ↓ RD2R^2_D ↑ MAED ↓ EbudRes(100) ↓ Unbounded 0.059±0.0070.059± 0.007 −90.15±6.14-90.15± 6.14 1.326±0.0441.326± 0.044 1.50±0.031.50± 0.03 ∑iβi(q)≤0.5 _i _i(q)≤ 0.5 0.071±0.0040.071± 0.004 0.703±0.0720.703± 0.072 0.061±0.0100.061± 0.010 1.58±0.041.58± 0.04 Figure 6: Damping identifiability and energy consistency across environments (q-only). Top row: learned damping field d(θ)d(θ) vs. ground truth (black). Pendulum: PHAST (KNOWN) recovers the sinusoidal profile d(θ)=d0+Δd|sinθ|d(θ)=d_0+ d\,|\! θ| near-exactly (RD2≈1R^2_D≈1); PARTIAL captures the shape but with a magnitude offset, illustrating the forecasting–identifiability trade-off. Double Pendulum: solid and dashed curves show per-joint damping (d1d_1, d2d_2). Despite chaotic dynamics and a coupled M(q)M(q), KNOWN recovers both joints’ damping profiles to high accuracy. Cart-Pole: solid and dashed curves distinguish constant viscous cart friction (dcd_c, flat) from angular wind damping on the pole (d(θ)d(θ), sinusoidal) — two qualitatively different dissipation mechanisms that PHAST disentangles from q-only data. UNKNOWN shows poor identifiability across all systems (RD2≪0R^2_D 0; Table 4). Baselines do not expose explicit damping fields (marked N/A). Bottom row: total energy H(t)H(t) during open-loop rollouts. In a dissipative system, energy must monotonically decrease; this is a necessary physical consistency check. PHAST (PARTIAL) most closely tracks the ground-truth energy decay across all three systems. Baselines either diverge (energy blowup) or collapse to zero, both indicating physically inconsistent dynamics. For PHAST, energy is computed using the observer’s learned velocity; for baselines, via finite differences. Additional studies (open-loop). Appendix H reports additional open-loop ablations (integrator choices, observer capacity, and damping constraints, including β¯ β sweeps) and a full-state comparison. 4.2 Closed-loop control (Energy–Casimir stabilization) Demarcation from open-loop. Open-loop forecasting evaluates autonomous rollouts (u=0u=0) and physical identifiability of learned fields. Closed-loop evaluation instead probes whether the learned pH structure supports passivity-based stabilization under feedback. In q-only settings, closed-loop performance depends critically on the quality of the online port/velocity estimate used for damping injection. Scope (closed-loop). Throughout this subsection we keep the controller structure and gains fixed and vary only the port/velocity estimate y y used for damping injection and (optionally) for the controller-state update. Unless stated otherwise, the plant is the true pendulum integrated by RK4. Setup and controller. We evaluate stabilization on a single-degree-of-freedom pendulum plant with state x=(q,p)x=(q,p) and q-only measurements qtmeas=qt+ϵtq_t^meas=q_t+ _t. The Energy–Casimir-style controller maintains an internal state ξ and uses a port estimate y^t≈yp(t) y_t≈ y_p(t) (which reduces to y^t=q˙^t y_t= q_t for the pendulum): ut u_t =−kc(ξt−q⋆)−dinjy^t, =-k_c( _t-q )\;-\;d_inj\, y_t, (30) ξt+1 _t+1 =ξt+Δt(y^t+kξ(qtmeas−ξt)), = _t+ t ( y_t+k_ξ\,(q_t^meas- _t) ), (31) where the kξk_ξ term is a predictor–corrector measurement correction that reduces ξ drift under noisy q-only feedback. When kξ=0k_ξ=0 and y^=yp y=y_p, the update reduces to a discrete approximation of ξ˙=yp ξ=y_p. Closed-loop protocol and metrics. We evaluate N rollouts from random initial conditions over horizon TctlT_ctl and report: (i) success rate (fraction of trials that converge to the target), (i) final wrapped-angle error, (i) control effort ∑t=0Tctl−1‖ut‖22 _t=0^T_ctl-1\|u_t\|_2^2, and (iv) port-estimation error t‖y^t−yp(t)‖E_t\| y_t-y_p(t)\| when ypy_p is available in simulation. Precise thresholds, horizons, and noise settings follow Appendix 5. Port-estimation variants. We compare five Energy–Casimir variants that share the same shaping gain kck_c, damping-injection gain dinjd_inj, and drift-correction gain kξk_ξ, differing only in how y^t y_t is obtained online: (i) oracle/full-state, (i) finite differences, (i) fixed-lag MAP smoothing (label-free), (iv) FD+TCN observer trained offline with noise augmentation, and (v) PHAST-trained observer trained only through next-step prediction. Appendix 5 reports additional control ablations (noise mismatch, near-stable regime, and model-based velocity from learned H H). Closed-loop results (summary). Sec. 5 shows that PHAST models support passivity-based stabilization: using the learned PHAST Hamiltonian to compute the port output y^=∂H^/∂p y=∂ H/∂ p achieves 100% stabilization success with lower control effort than oracle velocities (245 vs. 263), indicating the learned energy landscape is accurate enough for feedback. Under q-only sensing, the bottleneck shifts to port/velocity estimation quality rather than model accuracy; noise-aware observers reduce control effort by ∼9× 9× compared to finite differences. 5 Energy-Casimir Control Energy–Casimir control is a passivity-based approach for stabilizing port-Hamiltonian plants by interconnecting the plant with a dynamic controller and shaping a closed-loop storage function (Hamiltonian plus Casimir). Although PHAST is evaluated primarily as an open-loop forecaster in the main paper, port-Hamiltonian models are often used in feedback settings. We include this study to (i) show a minimal closed-loop wiring compatible with the port variables in Sec. 3.1, and (i) isolate a practical bottleneck in q-only control: the controller requires the velocity-like port output yp=yporty_p=y^port (for torque actuation, yp=q˙y_p= q), which must be estimated online from noisy position measurements. Closed-loop scope. We evaluate Energy–Casimir stabilization under q-only feedback. Across all control experiments we keep the controller structure and gains (kc,dinj,kξ)(k_c,d_inj,k_ξ) fixed and vary only the online port estimate y^≈yp y≈ y_p used for damping injection (and, when stated, in the ξ update). Unless stated otherwise, the plant is the true pendulum integrated by RK4. What is a Casimir function? A Casimir function (z)C(z) is a structural invariant of a port-Hamiltonian system: it satisfies ∇∈kerJcl⊤∩kerRcl ∈ J_cl ∩ R_cl, so it is conserved by the interconnection structure regardless of the energy function. For the mechanical plant–controller pair, (q,ξ)=q−ξC(q,ξ)=q-ξ acts as a “configuration lock”: the interconnection forces the controller state ξ to track the plant position q at all times (Sec. 3.5). Plant (pH)x˙=(Jp−Rp)∇Hp+Gpup x=(J_p-R_p)∇ H_p+G_p\,u_pyp=Gp⊤∇Hpy_p=G_p ∇ H_p(velocity-like; yp=q˙y_p= q for collocated)upu_pypy_p Controller (pH)Hc(ξ)=kc2(ξ−q⋆)2H_c(ξ)= k_c2(ξ-q )^2ξ˙=uc ξ=u_cyc=∇Hc(ξ)y_c=∇ H_c(ξ)(q-only: ξ˙=uc+kξ(q~−ξ) ξ=u_c+k_ξ( q-ξ))ucu_cycy_cq⋆q uc=ypu_c=y_pup=−yc+vu_p=-y_c+vDamping injectionv=vext−dinjy^v=v_ext-d_inj\, yq-only: y^≈yp y≈ y_pvvCasimir: =q−ξC=q- when uc=yp=q˙u_c=y_p= qPower balanceyp⊤up+yc⊤uc=yp⊤vy_p u_p+y_c u_c=y_p v(=0=0 when v=0v=0: no energy created) Figure 7: Energy–Casimir control as port-Hamiltonian interconnection. The plant (blue, left) and controller (green, right) are pH systems coupled through power ports: uc=ypu_c=y_p feeds plant velocity to the controller, up=−yc+vu_p=-y_c+v returns the shaped restoring force plus an auxiliary channel v. When v=0v=0 the cross-power cancels (yp⊤up+yc⊤uc=0y_p u_p+y_c u_c=0), and the Casimir =q−ξC=q-ξ is invariant. Damping injection v=−dinjy^v=-d_inj y (red, dashed) ensures H˙cl≤0 H_cl≤ 0; in q-only settings y^≈yp y≈ y_p is estimated, so estimation error can break exact power balance. Having established the theoretical framework in Sec. 3.5, we now evaluate it experimentally. All experiments use the forced-dynamics mode of the Energy–Casimir controller (Eqs. 28–29) on a single-degree-of-freedom pendulum. For the pendulum with collocated torque actuation, yp=q˙y_p= q, so the port estimate reduces to y^=q˙ y= q. MAP smoother observer (label-free). For the MAP smoother baseline, we estimate q˙^t q_t by MAP smoothing over a fixed-lag window of length w (with w≥2w≥ 2) of q-only measurements qt−w+1:tmeasq_t-w+1:t^meas. We unwrap the angles within the window (integrating wrapped increments) and solve the convex quadratic problem qsmooth q^smooth ∈argminq∑s=0w−11σ2‖qs−qsmeas‖22 ∈ _q\; _s=0^w-1 1σ^2\, \|q_s-q^meas_s \|_2^2 +∑s=0w−31Δt4σa2‖qs+2−2qs+1+qs‖22, +\; _s=0^w-3 1 t^4 _a^2\, \|q_s+2-2q_s+1+q_s \|_2^2, (32) which corresponds to a Gaussian measurement model and a Gaussian prior on discrete acceleration. We then take a backward difference on the smoothed sequence: q˙^t=qw−1smooth−qw−2smoothΔt. q_t= q^smooth_w-1-q^smooth_w-2 t. (33) This baseline requires only qmeasq^meas and the assumed noise scales (σ,σa)(σ, _a); it does not use q˙ q supervision. Unless otherwise stated, we fix σa=10.0 _a=10.0 for all MAP results. Main results. Table 7 summarizes the overall results under (i) no measurement noise and (i) moderate q-only measurement noise σ=0.01σ=0.01. Under noise, finite differences incur large velocity error and increased control effort; the learned FD+TCN observer reduces velocity error by ∼ 6–7× while preserving convergence. The PHAST-trained observer also preserves convergence at σ=0.01σ=0.01, but its velocity estimate remains noisy (closer to finite differences), which increases damping-injection effort near the target. A simple MAP smoother provides a competitive label-free baseline, reducing velocity error and final error without requiring q˙ q supervision. This emphasizes that closed-loop performance hinges on the quality of the port output estimate y=q˙y= q. Table 7: Energy–Casimir control on the true pendulum (q-only sensing). Mean over 100100 trials (5 regimes × 20). All methods use the same predictor–corrector ξ update (kξ=5k_ξ=5). σ Method Velocity estimate Success ↑ Final error ↓ Effort ↓ |q˙^−q˙|| q- q| ↓ 0.00 Oracle oracle 1.00 1.31×10−41.31× 10^-4 262.9 0.000 0.00 Finite differences FD 1.00 1.21×10−41.21× 10^-4 268.6 0.007 0.00 MAP smoother MAP smoother 1.00 1.16×10−41.16× 10^-4 268.6 0.007 0.00 FD+TCN observer FD+TCN 1.00 1.03×10−31.03× 10^-3 270.7 0.008 0.00 PHAST-trained observer PHAST observer 1.00 5.71×10−35.71× 10^-3 272.4 0.019 0.01 Oracle oracle 1.00 5.49×10−35.49× 10^-3 262.9 0.000 0.01 Finite differences FD 1.00 1.67×10−21.67× 10^-2 342.2 1.117 0.01 MAP smoother MAP smoother 1.00 1.17×10−21.17× 10^-2 324.8 0.133 0.01 FD+TCN observer FD+TCN 1.00 2.39×10−22.39× 10^-2 335.5 0.172 0.01 PHAST-trained observer PHAST observer 1.00 1.93×10−21.93× 10^-2 327.7 0.979 Model-based velocity from a learned PHAST model (full-state). To isolate the effect of model mismatch from partial observability, we also evaluate a controller that uses a learned PHAST Hamiltonian to form the port output estimate y^=∂H^/∂p y=∂ H/∂ p from the observed state (q,p)(q,p). Table 8 shows that this model-based velocity estimate is accurate enough to preserve stability, supporting the use of PHAST-learned dynamics for control when full state is available. Table 8: Energy–Casimir with PHAST-based velocity (q,p observed). Mean over 100100 trials (5 regimes × 20). Both methods use the same predictor–corrector ξ update (kξ=5k_ξ=5). σ Method Velocity estimate Success ↑ Final error ↓ Effort ↓ |q˙^−q˙|| q- q| ↓ 0.00 Oracle oracle 1.00 1.31×10−41.31× 10^-4 262.9 0.000 0.00 PHAST model-based ∂H^/∂p∂ H/∂ p 1.00 5.46×10−55.46× 10^-5 245.0 0.011 0.01 Oracle oracle 1.00 5.49×10−35.49× 10^-3 262.9 0.000 0.01 PHAST model-based ∂H^/∂p∂ H/∂ p 1.00 5.32×10−35.32× 10^-3 245.1 0.011 Near-target efficiency under noisy q-only feedback. The strongest effect of learning a noise-aware observer is visible near the target, where finite differences can inject high-frequency damping torques. Table 9 reports the near-stable regime (2020 trials) under σ=0.01σ=0.01. Table 9: Near-stable regime (σ=0.01σ=0.01). Learning a noise-aware q-only observer reduces control effort near the target by ∼ 9× vs. finite differences. Method Final error ↓ Effort ↓ |q˙^−q˙|| q- q| ↓ Oracle 0.0041 6.3 0.000 Finite differences 0.0170 80.1 1.117 MAP smoother 0.0133 10.3 0.095 FD+TCN observer 0.0137 8.6 0.130 PHAST-trained observer 0.0188 54.9 0.911 Observer-noise match ablation (limiting factor). To isolate the key failure mode under noisy q-only feedback, we train the FD+TCN observer under different noise distributions and evaluate at σ=0.01σ=0.01. Table 10 shows that noise-matching (or a narrow range anchored to deployment noise) is the dominant factor: training the observer at σ=0σ=0 and deploying at σ=0.01σ=0.01 recovers neither accurate velocities nor low-effort behavior. Table 10: Observer training noise vs. evaluation noise (σeval=0.01 _eval=0.01), near-stable regime. Q-only FD+TCN observer used inside the Energy–Casimir controller. Observer train noise Final error ↓ Effort ↓ |q˙^−q˙|| q- q| ↓ σtrain=0.00 _train=0.00 (mismatch) 0.0168 82.4 1.134 σtrain=0.01 _train=0.01 (matched) 0.0120 8.5 0.132 σtrain∼[0,0.01] _train [0,0.01] (anchored range) 0.0137 8.6 0.130 Takeaway. These experiments support a simple conclusion: in the q-only setting, closed-loop Energy–Casimir behavior is primarily limited by the velocity estimate used to form the port output y=q˙y= q. Noise-aware observer training is therefore a prerequisite for meaningful closed-loop validation under partial observability. High-noise stress test (σ=0.05σ=0.05). We additionally stress-test q-only feedback at higher measurement noise (σ=0.05σ=0.05). The oracle-velocity controller remains stable at 100% success, indicating that the controller structure itself is not the limiting factor. In contrast, q-only finite differences fail completely, while a noise-conditioned FD+TCN observer recovers partial stability. This setting remains challenging and is best viewed as a supplementary robustness probe rather than a main result. Table 11: High-noise stress test (σ=0.05σ=0.05). Mean over 100100 trials. For the FD+TCN controller we use a noise-conditioned observer trained with σtrain∼[0,0.05] _train [0,0.05]. Method Velocity estimate Success ↑ Final error ↓ Effort ↓ |q˙^−q˙|| q- q| ↓ Oracle oracle 1.00 2.75×10−22.75× 10^-2 263.6 0.000 Finite differences FD 0.00 8.36×10−28.36× 10^-2 2109.0 5.582 FD+TCN observer FD+TCN (noise) 0.26 7.70×10−27.70× 10^-2 396.6 0.659 Figure 8: Success rate breakdown at σ=0.05σ=0.05 for q-only FD+TCN control. Each bar reports success over 20 trials for the given initial-condition regime, varying only the observer training noise distribution. 6 Conclusion We presented PHAST, a port-Hamiltonian framework for learning dissipative dynamics that unifies three knowledge regimes—KNOWN, PARTIAL, and UNKNOWN—within a single architecture. The key technical contributions are low-rank parameterizations that guarantee D(q)⪰0D(q) 0 and M(q)≻0M(q) 0 by construction (Sec. 3.2), and optional damping-strength constraints that can prevent learned dissipation from absorbing model mismatch. Our main empirical finding is that forecasting and identifiability are distinct objectives that require explicit two-axis evaluation. Without damping bounds, models can achieve reasonable rollouts while learning physically meaningless parameters (R2<0R^2<0). On Windy Pendulum (q-only), PHAST (PARTIAL) improves 100-step rollout θ-wrap MSE from 0.435 (best baseline) to 0.092, while PHAST (KNOWN) recovers the true damping with R2=0.996R^2=0.996. This supports the broader takeaway that, on these q-only benchmarks, structural priors can matter more than additional capacity for long-horizon stability. Limitations. Our evaluation is primarily in the q-only setting, so performance depends on inferring a momentum-like latent from a short context window and can degrade under severe sensor noise or partial observability; scaling PHAST to high-DOF robotics and continuum systems remains open. In PARTIAL/UNKNOWN regimes, physical recovery and discrete-time stability are not guaranteed without calibrated anchors (e.g., damping bounds) and appropriate step sizes (Appendix H; Appendix A). Future work. Promising directions include: (i) extending PHAST to more complex dynamical systems such as multi-body and continuum systems; (i) learning damping bounds from data rather than requiring physics calibration; (i) integrating with online estimation and closed-loop planning (e.g., model-predictive control); (iv) extending our Energy–Casimir results (Sec. 5) to stabilization at non-natural equilibria and multi-DOF systems; and (v) test-time adaptation of physics components on the burn-in window, leveraging PHAST’s modular regime structure for instance-specific tuning of (V,M,D)(V,M,D) via meta-learned initializations. Acknowledgements This research was supported in part by the Peter O’Donnell Foundation, the Jim Holland–Backcountry Foundation, and in part by a grant from the Army Research Office accomplished under Cooperative Agreement Number W911NF-19-2-0333. Impact Statement This paper presents work whose goal is to advance the field of Machine Learning by developing physics-informed neural architectures that respect fundamental physical principles such as energy dissipation. Our framework enables more accurate and interpretable modeling of physical systems—including mechanical, electrical, molecular, thermal, gravitational, and ecological dynamics—with potential applications in robotics, simulation, scientific computing, and beyond. There are many potential societal consequences of our work, none which we feel must be specifically highlighted here. References Anonymous (2026) Anonymous. Learning generalized Hamiltonian dynamics with stability from noisy trajectory data. Concurrent Submission to ICML, 2026. Uploaded as Supplementary Material. Bajaj & Nguyen (2021) Bajaj, C. and Nguyen, M. Physics-informed neural networks via stochastic hamiltonian dynamics learning. arXiv preprint arXiv:2111.08108, 2021. Bajaj & Nguyen (2024) Bajaj, C. and Nguyen, M. Motion code: Robust time series classification and forecasting via sparse variational multi-stochastic processes learning. arXiv preprint arXiv:2402.14081, 2024. Bajaj et al. (2023) Bajaj, C., McLennan, L., Andeen, T., and Roy, A. Recipes for when physics fails: recovering robust learning of physics informed neural networks. Machine Learning: Science and Technology, 4(1):015013, 2023. Bajaj et al. (2024) Bajaj, C., Nguyen, M., and Li, C. Reinforcement learning for molecular dynamics optimization: A stochastic pontryagin maximum principle approach. arXiv preprint arXiv:2212.03320, 2024. Binney & Tremaine (2008) Binney, J. and Tremaine, S. Galactic Dynamics. Princeton University Press, 2nd edition, 2008. Boyer et al. (2025) Boyer, J., Rusch, T. K., and Rus, D. Learning to dissipate energy in oscillatory state-space models. arXiv preprint arXiv:2505.12171, 2025. Brantner et al. (2024) Brantner, B., de Romemont, G., Kraus, M., and Li, Z. Volume-preserving transformers for learning time series data with structure. arXiv preprint arXiv:2312.11166, 2024. Cho et al. (2014) Cho, K., Van Merriënboer, B., Gulcehre, C., Bahdanau, D., Bougares, F., Schwenk, H., and Bengio, Y. Learning phrase representations using rnn encoder-decoder for statistical machine translation. arXiv preprint arXiv:1406.1078, 2014. Cranmer et al. (2020) Cranmer, M., Greydanus, S., Hoyer, S., Battaglia, P., Spergel, D., and Ho, S. Lagrangian neural networks. arXiv preprint arXiv:2003.04630, 2020. Desai et al. (2021) Desai, S. A., Mattheakis, M., Sondak, D., Protopapas, P., and Roberts, S. J. Port-hamiltonian neural networks for learning explicit time-dependent dynamical systems. Physical Review E, 104(3):034312, 2021. Eidnes & Riemer-Sørensen (2023) Eidnes, S. and Riemer-Sørensen, S. Pseudo-Hamiltonian neural networks for learning partial differential equations. Journal of Computational Physics, 479:112007, 2023. Ellendula et al. (2025) Ellendula, A. S., Wang, Y., Nguyen, M., and Bajaj, C. Grl-snam: Geometric reinforcement learning with path differential hamiltonians for simultaneous navigation and mapping in unknown environments. arXiv preprint arXiv:2601.00116, 2025. Greydanus et al. (2019) Greydanus, S., Dzamba, M., and Yosinski, J. Hamiltonian neural networks. Advances in neural information processing systems, 32, 2019. Gu & Dao (2023) Gu, A. and Dao, T. Mamba: Linear-time sequence modeling with selective state spaces. arXiv preprint arXiv:2312.00752, 2023. Gu et al. (2022) Gu, A., Goel, K., and Ré, C. Efficiently modeling long sequences with structured state spaces. arXiv preprint arXiv:2111.00396, 2022. Hairer et al. (2006) Hairer, E., Lubich, C., and Wanner, G. Geometric Numerical Integration: Structure-Preserving Algorithms for Ordinary Differential Equations. Springer, 2006. Hasani et al. (2024) Hasani, R., Wieser, E., Lechner, M., Kim, T. A., Amini, A., Rus, D., and Grosu, R. Linear oscillatory state space models. arXiv preprint arXiv:2410.03943, 2024. Horowitz & Hill (2015) Horowitz, P. and Hill, W. The Art of Electronics. Cambridge University Press, 3rd edition, 2015. Incropera et al. (2007) Incropera, F. P., DeWitt, D. P., Bergman, T. L., and Lavine, A. S. Fundamentals of Heat and Mass Transfer. John Wiley & Sons, 6th edition, 2007. Jin et al. (2020) Jin, P., Zhang, Z., Zhu, A., Tang, Y., and Karniadakis, G. E. Sympnets: Intrinsic structure-preserving symplectic networks for identifying hamiltonian systems. Neural Networks, 132:166–179, 2020. Jones (1924) Jones, J. E. On the determination of molecular fields.—I. From the equation of state of a gas. Proceedings of the Royal Society of London. Series A, 106(738):463–477, 1924. McLennan et al. (2025) McLennan, L., Wang, Y., Farell, R., Nguyen, M., and Bajaj, C. Learning generalized hamiltonian dynamics with stability from noisy trajectory data. arXiv preprint arXiv:2509.07280, 2025. Mehl et al. (2018) Mehl, C., Mehrmann, V., and Wojtylak, M. Linear algebra properties of dissipative Hamiltonian descriptor systems. SIAM Journal on Matrix Analysis and Applications, 39(3):1489–1519, 2018. Murray (2002) Murray, J. D. Mathematical Biology: I. An Introduction. Springer, 3rd edition, 2002. Nguyen & Bajaj (2025a) Nguyen, M. and Bajaj, C. A differential and pointwise control approach to reinforcement learning. arXiv preprint arXiv:2404.15617, 2025a. Nguyen & Bajaj (2025b) Nguyen, M. and Bajaj, C. Stochastic differential policy optimization: A rough path approach to reinforcement learning. In Theory of AI for Scientific Computing Workshop, 2025b. Smith et al. (2023) Smith, J. T., Warrington, A., and Linderman, S. W. Simplified state space layers for sequence modeling. arXiv preprint arXiv:2208.04933, 2023. Sosanya & Greydanus (2022) Sosanya, A. and Greydanus, S. Dissipative hamiltonian neural networks: Learning dissipative and conservative dynamics separately. In International Conference on Learning Representations Workshop on Physics for Machine Learning, 2022. Strang (1968) Strang, G. On the construction and comparison of difference schemes. SIAM Journal on Numerical Analysis, 5(3):506–517, 1968. Strogatz (2018) Strogatz, S. H. Nonlinear Dynamics and Chaos: With Applications to Physics, Biology, Chemistry, and Engineering. CRC Press, 2nd edition, 2018. Van Der Schaft & Jeltsema (2014) Van Der Schaft, A. and Jeltsema, D. Port-Hamiltonian systems theory: An introductory overview. Now Publishers Inc, 2014. Vaswani et al. (2017) Vaswani, A., Shazeer, N., Parmar, N., Uszkoreit, J., Jones, L., Gomez, A. N., Kaiser, Ł., and Polosukhin, I. Attention is all you need. Advances in neural information processing systems, 30, 2017. Zhong et al. (2020) Zhong, Y. D., Dey, B., and Chakraborty, A. Dissipative SymODEN: Encoding Hamiltonian dynamics with dissipation and control into deep learning. arXiv preprint arXiv:2002.08860, 2020. Appendix A Mathematical Details Organization. Appendix A provides supplementary mathematical details (passivity, discrete-time energy under splitting, and Woodbury identities). Appendix B describes the benchmark systems and their Hamiltonians/damping laws. Appendix C and Appendix D detail architectural components and parameterizations. Appendix E summarizes losses and hyperparameters. Appendix F defines evaluation metrics and q-only diagnostics. Sec. G and Appendix H collect additional tables and ablations. Finally, Appendix I lists notation, and Sec. 5 reports a separate closed-loop Energy–Casimir control study. A.1 Passivity Proof Theorem A.1 (Energy Balance and Passivity). Consider port-Hamiltonian dynamics with input x˙=(J−R)∇H(x)+Gu,x=(q,p), x=(J-R)∇ H(x)+Gu, x=(q,p), (34) and define the conjugate port output yport:=G⊤∇H(x)y^port:=G ∇ H(x) (Sec. 3.1; not to be confused with the q-only observation yt=qty_t=q_t). Then the energy balance is dHdt=−(∇H)⊤R(∇H)+yport⊤u≤yport⊤u. dHdt=-(∇ H) R(∇ H)+y^port u≤y^port u. (35) For mechanical systems with R=diag(0,D(q))R=diag(0,D(q)) and D(q)⪰0D(q) 0, this reduces to dHdt=−v⊤D(q)v+yport⊤u≤yport⊤u,v:=∂H∂p(=M−1pfor our separable Hamiltonian). dHdt=-v D(q)\,v+y^port u≤y^port u, v:= ∂ H∂ p\;(=M^-1p\ for our separable Hamiltonian). (36) In particular, the unforced system (u=0u=0) is passive with nonincreasing energy. Proof. Expanding the time derivative of the Hamiltonian: dHdt dHdt =∇H⊤x˙=∇H⊤(J−R)∇H+∇H⊤Gu =∇ H x=∇ H (J-R)∇ H+∇ H Gu (37) =∇H⊤J∇H⏟=0 (skew-sym)−∇H⊤R∇H = ∇ H J∇ H_=0 (skew-sym)-∇ H R∇ H (38) +(G⊤∇H)⊤u⏟=yport⊤u. + (G ∇ H) u_=y^port u. (39) With the block structure R=diag(0,D(q))R=diag(0,D(q)) and ∇pH=v _pH=v: dHdt=−v⊤D(q)v+yport⊤u≤yport⊤u, dHdt=-v D(q)\,v+y^port u≤y^port u, (40) since D(q)⪰0D(q) 0 by construction. ∎ A.2 Energy Budget Under Strang Splitting Proposition A.2. Let ΦHΔt _H t denote the (exact) flow of the conservative subsystem and ΦDΔt _D t the (exact) flow of the dissipative subsystem with q held fixed. Under Strang splitting ΦΔt=ΦDΔt/2∘ΦHΔt∘ΦDΔt/2 _ t= _D t/2 _H t _D t/2 with sufficiently small Δt t: H(q+,p+)−H(q,p)=−∫0Δtv(t)⊤D(q(t))v(t)t=−Δt⋅v⊤D(q)v+O(Δt2),H(q^+,p^+)-H(q,p)=- _0 tv(t) D(q(t))\,v(t)\,dt=- t· v D(q)\,v+O( t^2), (41) where v(t)=M−1p(t)v(t)=M^-1p(t) and the O(Δt2)O( t^2) term reflects variation of (q(t),v(t))(q(t),v(t)) over the step. Remark (discrete-time dissipation). Our dissipative half-step uses an explicit (Euler) update on the momentum dynamics p˙=−D(q)v p=-D(q)\,v: p←p−Δt2D(q)v,v=M−1p.p← p- t2\,D(q)\,v, v=M^-1p. (42) Since q is fixed in this half-step, only the kinetic energy changes. Writing h=Δt/2h= t/2, a direct expansion gives H(q,p+)−H(q,p)=−hv⊤D(q)v+h22v⊤D(q)M−1D(q)v,H(q,p^+)-H(q,p)=-h\,v D(q)\,v+ h^22\,v D(q)\,M^-1\,D(q)\,v, (43) so discrete-time energy monotonicity is step-size dependent. In particular, a sufficient condition for H(q,p+)≤H(q,p)H(q,p^+)≤ H(q,p) for all v is s:=Δt2λmax(D(q))λmin(M)≤2,s:= t2\, _ (D(q)) _ (M)≤ 2, (44) which can fail in stiff regimes. Bounding ∑i=1rβi(q)≤β¯ _i=1^r _i(q)≤ β (where r is the number of rank-1 terms; Eq. 11) controls λmax(D) _ (D) and improves stability in practice. Separately, our conservative map ΦHΔt _H t is implemented with symplectic integrators (leapfrog / implicit midpoint), so the true Hamiltonian need not be exactly preserved even when the damping half-step is stable; we therefore treat passivity violations and energy-budget residuals as empirical diagnostics of the net discrete-time behavior. A.3 Woodbury Identity for Mass Inverse For the low-rank mass M=Λ+UU⊤M= +U where Λ=diag(d) =diag(d) and U∈ℝn×rU ^n× r: M−1=Λ−1−Λ−1U(Ir+U⊤Λ−1U)−1U⊤Λ−1.M^-1= ^-1- ^-1U\,(I_r+U ^-1U)^-1\,U ^-1. (45) The r×r× r matrix (Ir+U⊤Λ−1U)(I_r+U ^-1U) is inverted once; total cost is O(nr2)O(nr^2) instead of O(n3)O(n^3). More precisely, the dominant terms are O(nr2+r3)O(nr^2+r^3) (with r≪nr n in all our experiments). For the log-determinant (needed if extending to configuration-dependent/Riemannian Hamiltonians): log|M|=log|Λ|+log|Ir+U⊤Λ−1U|, |M|= | |+ |I_r+U ^-1U|, (46) where the first term is O(n)O(n) and the second is O(r3)O(r^3). A.4 Spectral Properties of the Linearized System Consider the linearized port-Hamiltonian dynamics x˙=(J−R)Qx x=(J-R)Q\,x, where Q=Q⊤≻0Q=Q 0 is the Hessian of H at an equilibrium and R⪰0R 0 is PSD. Since J is skew-symmetric, we have A:=(J−R)Q.A:=(J-R)\,Q. (47) The eigenvalues of A lie in the closed left half-plane, i.e., Re(λi)≤0Re( _i)≤ 0 for all i. Proof. Let (λ,z)(λ,z) be an eigenpair of A, so (J−R)Qz=λz(J-R)Qz=λ z. Define w=Q1/2zw=Q^1/2z (well-defined since Q≻0Q 0), so Q1/2(J−R)Q1/2w=λwQ^1/2(J-R)Q^1/2w=λ w. Then Re(λ)=Re(w∗Q1/2(J−R)Q1/2w)/(w∗w)Re(λ)=Re(w^*Q^1/2(J-R)Q^1/2w)/(w^*w). Since w∗Q1/2JQ1/2w^*Q^1/2JQ^1/2w is purely imaginary (by skew-symmetry), Re(λ)=−w∗Q1/2RQ1/2w/(w∗w)≤0Re(λ)=-w^*Q^1/2RQ^1/2w/(w^*w)≤ 0 since R⪰0R 0. ∎ Conservative case (R=0R=0): All eigenvalues are purely imaginary, corresponding to oscillatory energy exchange (phase-space rotations). Dissipative case (R≻0R 0): Eigenvalues shift into the open left half-plane (e.g., −σ±iω-σ± iω with σ>0σ>0), yielding damped oscillations. See Mehl et al. (2018) for a comprehensive treatment. Appendix B Environment Descriptions This section provides detailed descriptions of the systems used in our q-only benchmarks, including their state spaces, Hamiltonians, and damping models. We cover four mechanical systems (Secs. B.1–B.4) and five non-mechanical systems (Sec. B.5). B.1 Single Pendulum m,Im,ICoMθℓc _cℓ (θ)θ˙d(θ) θxxyy Figure 9: Single pendulum. A point mass m at distance ℓ from a fixed pivot, with moment of inertia I about the pivot. The angle θ is measured from the downward vertical. Position-dependent damping d(θ)d(θ) acts at the pivot. Configuration space. The single pendulum has one degree of freedom: the angle θ∈1θ ^1 measured from the downward vertical (so θ=0θ=0 corresponds to the stable equilibrium with the pendulum hanging straight down). State space. The phase-space state is x=(θ,p)∈1×ℝx=(θ,p) ^1×R, where p=Iθ˙p=I θ is the angular momentum and I is the moment of inertia about the pivot. Kinematics. The Cartesian position of the center of mass is: =[ℓcsinθ−ℓccosθ],r= bmatrix _c θ\\ - _c θ bmatrix, (48) where ℓc _c is the distance from the pivot to the center of mass. Kinetic energy. For rotation about the fixed pivot: Ekin=12Iθ˙2=p22I,E_kin= 12I θ^2= p^22I, (49) where I=mℓc2I=m _c^2 for a point mass (or includes rotational inertia for an extended body). Potential energy. Taking the lowest point (θ=0θ=0) as the zero reference: V(θ)=mgℓc(1−cosθ).V(θ)=mg _c(1- θ). (50) Hamiltonian. The total energy is: H(θ,p)=p22I+mgℓc(1−cosθ).H(θ,p)= p^22I+mg _c(1- θ). (51) Equations of motion. The conservative dynamics are: θ˙ θ =∂H∂p=pI, = ∂ H∂ p= pI, (52) p˙ p =−∂H∂θ=−mgℓcsinθ. =- ∂ H∂θ=-mg _c θ. (53) Damping model. We consider three damping variants: • Conservative: d(θ)=0d(θ)=0. • Constant damping: d(θ)=γd(θ)=γ (constant viscous friction, γ=0.5γ=0.5). • Windy (position-dependent): d(θ)=d0+Δd|sinθ|d(θ)=d_0+ d\,| θ|, where d0=0.3d_0=0.3 and Δd=0.5 d=0.5. Windy damping: physical interpretation. The “windy” damping model d(θ)=d0+Δd|sinθ|d(θ)=d_0+ d\,| θ| captures position-dependent air resistance: when the pendulum is horizontal (θ=±π/2θ=±π/2), it presents maximum cross-sectional area to the airflow, yielding maximum damping dmax=d0+Δd=0.8d_ =d_0+ d=0.8. When vertical (θ=0θ=0 or π), the cross-section is minimal, yielding dmin=d0=0.3d_ =d_0=0.3. This is the primary benchmark for testing PHAST’s ability to learn configuration-dependent dissipation D(q)D(q). θ (θ)d(θ)−π-π−π2- π20π2 π2π 0=0.3d_0=0.3d0+Δd=0.8d_0+ d=0.8d(θ)=d0+Δd|sinθ|d(θ)=d_0+ d\,| θ|min dragmax drag Figure 10: Windy damping profile. The position-dependent damping coefficient d(θ)=d0+Δd|sinθ|d(θ)=d_0+ d\,| θ| varies between d0=0.3d_0=0.3 (vertical) and d0+Δd=0.8d_0+ d=0.8 (horizontal), modeling air resistance that depends on the pendulum’s cross-sectional area. Port-Hamiltonian dynamics with windy damping. The dissipative equations of motion are: θ˙ θ =pI, = pI, (54) p˙ p =−mgℓcsinθ−d(θ)⋅pI. =-mg _c θ-d(θ)· pI. (55) The energy dissipation rate is: dHdt=−d(θ)⋅v2=−(d0+Δd|sinθ|)⋅(pI)2≤0, dHdt=-d(θ)· v^2=- (d_0+ d\,| θ| )· ( pI )^2≤ 0, (56) which is always non-positive (passive) since d(θ)≥d0>0d(θ)≥ d_0>0. Simulator parameters. In our experiments: m=1m=1, ℓc=ℓ=1 _c= =1, g=9.81g=9.81, I=mℓc2=1I=m _c^2=1. B.2 Cart-Pole mcm_cxxmpm_pθℓ (θ)q˙d(θ) qxxyy Figure 11: Cart-pole system. A cart of mass mcm_c moves on a horizontal rail while a pole of mass mpm_p (concentrated at distance ℓ from the pivot) rotates about the cart. The angle θ is measured from the downward vertical (hanging configuration: θ=0θ=0). Windy damping acts isotropically on q˙ q via a scalar coefficient d(θ)d(θ). Configuration space. The cart-pole has two degrees of freedom: q=(x,θ)∈ℝ×1q=(x,θ) ×S^1, where x is the cart position and θ is the pole angle measured from the downward vertical (stable at θ=0θ=0). State space. The phase-space state is x=(q,p)=(x,θ,px,pθ)∈ℝ×1×ℝ2x=(q,p)=(x,θ,p_x,p_θ) ×S^1×R^2. Kinematics. The positions of the cart and pole center of mass are: c _c =[x0], = bmatrixx\\ 0 bmatrix, p _p =[x+ℓsinθ−ℓcosθ]. = bmatrixx+ θ\\ - θ bmatrix. (57) Kinetic energy. The velocities are ˙c=(x˙,0)⊤ r_c=( x,0) and: ˙p=[x˙+ℓcosθ⋅θ˙ℓsinθ⋅θ˙]. r_p= bmatrix x+ θ· θ\\ θ· θ bmatrix. (58) The total kinetic energy is: Ekin E_kin =12mcx˙2+12mp‖˙p‖2 = 12m_c x^2+ 12m_p\| r_p\|^2 (59) =12(mc+mp)x˙2+mpℓcosθ⋅x˙θ˙+12mpℓ2θ˙2, = 12(m_c+m_p) x^2+m_p θ· x θ+ 12m_p ^2 θ^2, (60) which can be written as Ekin=12q˙⊤M(q)q˙E_kin= 12 q M(q) q with the configuration-dependent mass matrix: M(θ)=[mc+mpmpℓcosθmpℓcosθmpℓ2].M(θ)= bmatrixm_c+m_p&m_p θ\\ m_p θ&m_p ^2 bmatrix. (61) Potential energy. Taking the hanging configuration (θ=0θ=0) as the zero reference: V(θ)=mpgℓ(1−cosθ).V(θ)=m_pg (1- θ). (62) Hamiltonian. With generalized momenta p=M(θ)q˙p=M(θ) q: H(q,p)=12p⊤M(θ)−1p+mpgℓ(1−cosθ).H(q,p)= 12p M(θ)^-1p+m_pg (1- θ). (63) Separable approximation. In our PHAST experiments, we use a separable Hamiltonian with constant M (Sec. 1), which is an approximation. The true Cart-Pole has configuration-dependent inertia as shown in Eq. (61). Damping model. We use windy (position-dependent) damping on the pole angle: d(θ)=d0+Δd|sinθ|,d(θ)=d_0+ d\,| θ|, (64) with d0=0.3d_0=0.3 and Δd=0.5 d=0.5. The damping acts on both the cart and pole velocities as a scalar coefficient: p˙=−∇qH−d(θ)⋅q˙. p=- _qH-d(θ)· q. (65) This models increased air resistance when the pole is horizontal. The energy dissipation rate is: dHdt=−d(θ)⋅‖q˙‖2≤0. dHdt=-d(θ)·\| q\|^2≤ 0. (66) Simulator parameters. In our experiments: mc=1m_c=1, mp=1m_p=1, ℓ=1 =1, g=9.81g=9.81, Δt=0.02 t=0.02. B.3 Harmonic Oscillator ω 1m_1q1q_1γω 2m_2q2q_2γ =2n=2 independentharmonic oscillatorswith unit mass Figure 12: Harmonic oscillator. Two independent harmonic oscillators with natural frequency ω and unit mass. Coordinates q1,q2q_1,q_2 measure displacements from equilibrium. Optional viscous damping γ acts on each oscillator. Configuration space. The oscillator benchmark uses n degrees of freedom (we use n=2n=2): q=(q1,…,qn)∈ℝnq=(q_1,…,q_n) ^n, where qiq_i is the displacement of oscillator i from equilibrium. State space. The phase-space state is x=(q,p)∈ℝ2nx=(q,p) ^2n, where p=q˙p= q (unit mass). Kinetic energy. With unit mass (M=IM=I): Ekin=12‖p‖2=12∑i=1npi2.E_kin= 12\|p\|^2= 12 _i=1^np_i^2. (67) Potential energy. Simple harmonic potential with natural frequency ω: V(q)=12ω2‖q‖2=12ω2∑i=1nqi2.V(q)= 12ω^2\|q\|^2= 12ω^2 _i=1^nq_i^2. (68) Hamiltonian. The total energy is: H(q,p)=12‖p‖2+12ω2‖q‖2.H(q,p)= 12\|p\|^2+ 12ω^2\|q\|^2. (69) Equations of motion. q˙ q =p, =p, (70) p˙ p =−ω2q−γp, =-ω^2q-γ p, (71) where γ is the viscous damping coefficient. Damping model. • Conservative: γ=0γ=0 (energy preserved). • Damped: γ=0.1γ=0.1 (constant viscous damping). Simulator parameters. In our experiments: n=2n=2 (degrees of freedom), ω=1ω=1 (natural frequency), Δt=0.02 t=0.02. B.4 Double Pendulum m1m_1ℓ1 _1m2m_2ℓ2 _2θ1 _1θ2 _2ggxxyybbbb Figure 13: Double pendulum (point-mass model). Two point masses m1,m2m_1,m_2 at the ends of massless rods of lengths ℓ1,ℓ2 _1, _2. The angles θ1 _1 and θ2 _2 are both measured from the downward vertical (absolute angles). Optional viscous damping b acts at both joints. Configuration space. The double pendulum has two degrees of freedom: q=(θ1,θ2)∈1×1q=( _1, _2) ^1×S^1, where θ1 _1 is the angle of link 1 from the downward vertical, and θ2 _2 is the angle of link 2 from the downward vertical (both absolute angles). State space. The phase-space state is x=(q,p)=(θ1,θ2,p1,p2)∈(1)2×ℝ2x=(q,p)=( _1, _2,p_1,p_2)∈(S^1)^2×R^2, where p=M(q)q˙p=M(q) q are the canonical momenta. Kinematics. Using shorthand si=sinθis_i= _i, ci=cosθic_i= _i: 1 _1 =[ℓ1s1−ℓ1c1], = bmatrix _1s_1\\ - _1c_1 bmatrix, 2 _2 =[ℓ1s1+ℓ2s2−ℓ1c1−ℓ2c2]. = bmatrix _1s_1+ _2s_2\\ - _1c_1- _2c_2 bmatrix. (72) Kinetic energy. For point masses at the ends of massless rods: Ekin=12m1‖˙1‖2+12m2‖˙2‖2=12q˙⊤M(q)q˙,E_kin= 12m_1\| r_1\|^2+ 12m_2\| r_2\|^2= 12 q M(q) q, (73) where the configuration-dependent mass matrix is: M(q)=[(m1+m2)ℓ12m2ℓ1ℓ2cos(θ1−θ2)m2ℓ1ℓ2cos(θ1−θ2)m2ℓ22].M(q)= bmatrix(m_1+m_2) _1^2&m_2 _1 _2 ( _1- _2)\\ m_2 _1 _2 ( _1- _2)&m_2 _2^2 bmatrix. (74) Note that M depends on the angle difference (θ1−θ2)( _1- _2). Potential energy. Taking the pivot level as the zero reference: V(θ1,θ2)=−(m1+m2)gℓ1cosθ1−m2gℓ2cosθ2.V( _1, _2)=-(m_1+m_2)g _1 _1-m_2g _2 _2. (75) Hamiltonian. With generalized momenta p=M(q)q˙p=M(q) q: H(q,p)=12p⊤M(q)−1p+V(θ1,θ2).H(q,p)= 12p M(q)^-1p+V( _1, _2). (76) Equations of motion. The Euler–Lagrange equations yield (in manipulator form): M(q)q¨+C(q,q˙)q˙+G(q)=−bq˙,M(q) q+C(q, q) q+G(q)=-b q, (77) where C(q,q˙)C(q, q) contains Coriolis/centrifugal terms involving sin(θ1−θ2) ( _1- _2), G(q)=∇qV(q)G(q)= _qV(q) is the gravity vector, and b is the viscous damping coefficient applied equally to both joints. Separable approximation. In our PHAST experiments, we use a separable Hamiltonian with constant M (Sec. 1). This is an approximation; the true double pendulum has configuration-dependent inertia as shown in Eq. (74). Damping model. • Conservative: b=0b=0. • Damped: b=0.2b=0.2 (viscous damping on both joints). Simulator parameters. In our experiments: m1=m2=1m_1=m_2=1, ℓ1=ℓ2=1 _1= _2=1, g=9.81g=9.81, Δt=0.01 t=0.01. B.5 Detailed Port-Hamiltonian Mappings for Table 1 Section 1 walks through three representative entries (simple pendulum, cart-pole, RLC circuit). Here we provide the complete (V,M,D)(V,M,D) decomposition for every system in Table 1, organized by knowledge regime. KNOWN regime. In these systems, the potential V and mass M are fully specified by the physics; only the dissipation D (or a small number of scalar parameters) must be learned. This makes the inverse problem well-posed: recovering D from trajectory data has a unique solution. kkbbmmxx Figure 14: Spring–mass–damper. A mass m connected to a wall by a spring (stiffness k) and dashpot (damping b) in parallel. Displacement x is measured from equilibrium. Spring–mass (KNOWN, mechanical). The damped harmonic oscillator is the simplest non-trivial port-Hamiltonian system and serves as the “hello world” for our benchmark: if a method cannot recover the damping coefficient of a linear spring, it will struggle with anything more complex. Despite its simplicity, the system already exhibits the conservative–dissipative split that defines the port-Hamiltonian framework: the spring stores energy, the dashpot removes it, and the mass mediates how quickly energy converts between potential and kinetic forms. • State: q=xq=x (displacement from equilibrium), p=mx˙p=m x (linear momentum). • Potential: V(x)=12kx2V(x)= 12kx^2 — a quadratic well whose gradient −kx-kx is the familiar restoring force. • Mass: M=mM=m — a scalar constant. In the port-Hamiltonian picture, m controls the “exchange rate” between momentum and velocity: x˙=p/m x=p/m. • Damping: D=bD=b — viscous friction (Fd=−bx˙F_d=-b x). Energy dissipation is H˙=−bx˙2≤0 H=-b x^2≤ 0, the defining passivity inequality. • Hamiltonian: H(x,p)=12kx2+p22mH(x,p)= 12kx^2+ p^22m. • Identifiability: With V and m both given, the only unknown is b—recoverable from the decay rate of oscillations. 11r1r_122r2r_233r3r_3r12r_12r13r_13r23r_23D=γI2N [rgb]1,0.35,0.35 [named]pgfstrokecolorrgb1,0.35,0.35D=γ I_2Nisotropic dampingM=diag(miI2) [rgb]0,0.5,0 [named]pgfstrokecolorrgb0,0.5,0M=diag(m_iI_2)atomic masses (known)V(r)V(r)rrrepulsionr−12r^-12attraction−r−6-r^-6reqr_eqσϵ (q)=∑i<j4ϵ[(σrij)12−(σrij)6] [rgb]0,0,0.6 [named]pgfstrokecolorrgb0,0,0.6V(q)= _i<j4ε\! [ ( σr_ij )^\!12- ( σr_ij )^\!6 ]Assumptions: 3 particles in 2D ⋅· softened potential (ϵsoft=0.1 _soft=0.1)⋅· isotropic damping ⋅· ϵ=1,σ=1,γ=0.1ε=1,\,σ=1,\,γ=0.1 Figure 15: Lennard–Jones cluster (Jones, 1924). N=3N=3 particles interact pairwise via the LJ potential. Inset: the repulsive core (r−12r^-12, red shading) dominates at r<σr<σ; the attractive well (−r−6-r^-6, blue shading) has depth ϵε near the equilibrium distance req=21/6σr_eq=2^1/6σ. Dissipation is modeled by isotropic velocity-proportional damping (D=γID=γ I). Lennard–Jones cluster (KNOWN, molecular dynamics). Molecular dynamics is a natural stress test for structure-preserving learners: the Lennard–Jones potential has a steep repulsive wall that amplifies integration errors, yet energy must remain bounded for the cluster to stay intact. A port-Hamiltonian integrator respects this energy budget by construction, making it a principled alternative to hand-tuned thermostats. We simulate N=3N=3 particles in 2D interacting through the classical LJ pair potential (Jones, 1924), with a velocity-proportional drag that models coupling to an implicit heat bath. • State: q=(r1,…,rN)∈ℝ2Nq=(r_1,…,r_N) ^2N (particle positions in 2D), p=(m1r˙1,…,mNr˙N)p=(m_1 r_1,…,m_N r_N) (momenta). • Potential: V(q)=∑i<j4ϵ[(σ/rij)12−(σ/rij)6] [rgb]0,0,0.6 [named]pgfstrokecolorrgb0,0,0.6V(q)= _i<j4ε [(σ/r_ij)^12-(σ/r_ij)^6 ], where rij=‖ri−rj‖r_ij=\|r_i-r_j\|, ϵε is the well depth, and σ is the particle diameter. The r−12r^-12 term encodes Pauli repulsion (atoms cannot overlap); the r−6r^-6 term encodes van der Waals attraction (induced-dipole interactions at longer range). Together they create a potential well at req=21/6σr_eq=2^1/6σ that traps the cluster. • Mass: M=diag(m1I2,…,mNI2) [rgb]0,0.5,0 [named]pgfstrokecolorrgb0,0.5,0M=diag(m_1I_2,…,m_NI_2) — each particle carries a known scalar mass; the block-diagonal structure reflects the absence of cross-particle inertial coupling. • Damping: D=γI2N [rgb]1,0.35,0.35 [named]pgfstrokecolorrgb1,0.35,0.35D=γ I_2N — isotropic friction that damps every degree of freedom equally. Physically, γ controls how quickly kinetic energy drains to the heat bath; in the conservative limit γ→0γ→0 the cluster oscillates forever. • Identifiability: All physics (V, M) is given, so the inverse problem reduces to recovering a single scalar γ. RRdamping D=RD=RLLmass M=LM=LCCpotentialqCq_Ci=q˙Ci= q_C Port-Hamiltonian mapping State: q=qCq=q_C (charge), p=Lip=Li (flux) V(q)=qC22CV(q)= q_C^22C (capacitive energy) M=LM=L (inductance == “mass”) D=RD=R (resistance == damping) H(q,p)=qC22C+p22L=qC22C+Li22H(q,p)= [rgb]0,0,0.6 [named]pgfstrokecolorrgb0,0,0.6 q_C^22C+ [rgb]0,0.55,0 [named]pgfstrokecolorrgb0,0.55,0 p^22L= [rgb]0,0,0.6 [named]pgfstrokecolorrgb0,0,0.6 q_C^22C+ [rgb]0,0.55,0 [named]pgfstrokecolorrgb0,0.55,0 Li^22Assumptions: ideal lumped elements ⋅· no mutual inductance ⋅· L=1,C=1,R=0.5L=1,\,C=1,\,R=0.5 Figure 16: Series RLC circuit (Horowitz & Hill, 2015). Charge qCq_C on the capacitor is the generalized coordinate; flux linkage p=Lip=Li is the conjugate momentum. Color coding shows the port-Hamiltonian mapping: capacitance C → potential V, inductance L → mass M, resistance R → damping D. The electrical energy H=qC2/2C+Li2/2H=q_C^2/2C+Li^2/2 is the storage function. Series RLC circuit (KNOWN, electrical). The series RLC circuit is perhaps the cleanest non-mechanical port-Hamiltonian system: every textbook component maps one-to-one onto the (V,M,D)(V,M,D) triple, making it an ideal sanity check for structure-preserving learning in a domain where the analogy is exact rather than approximate (Horowitz & Hill, 2015). • State: q=qCq=q_C (charge on the capacitor), p=Lip=Li (magnetic flux linkage through the inductor, conjugate to qCq_C). Current i=q˙Ci= q_C is the generalized velocity. • Potential: V(q)=qC2/2C [rgb]0,0,0.6 [named]pgfstrokecolorrgb0,0,0.6V(q)=q_C^2/2C — electrostatic energy stored in the capacitor. The restoring “force” −∂V/∂qC=−qC/C-∂ V/∂ q_C=-q_C/C is simply the capacitor voltage, driving current back around the loop. • Mass: M=L [rgb]0,0.55,0 [named]pgfstrokecolorrgb0,0.55,0M=L — inductance plays the role of inertia. Just as a heavy mass resists changes in velocity, a large inductance resists changes in current (Lenz’s law). • Damping: D=R [rgb]1,0.35,0.35 [named]pgfstrokecolorrgb1,0.35,0.35D=R — resistance dissipates electrical energy as Joule heat, exactly analogous to viscous friction converting kinetic energy to thermal energy. • Identifiability: With V and L given, the resistance R is uniquely recoverable from the oscillation decay envelope. T1T_1heat capacity c1c_1V1=12c1T12 [rgb]0,0,0.6 [named]pgfstrokecolorrgb0,0,0.6V_1= 12c_1T_1^2M1=τ1 [rgb]0,0.55,0 [named]pgfstrokecolorrgb0,0.55,0M_1= _1T2T_2heat capacity c2c_2V2=12c2T22 [rgb]0,0,0.6 [named]pgfstrokecolorrgb0,0,0.6V_2= 12c_2T_2^2M2=τ2 [rgb]0,0.55,0 [named]pgfstrokecolorrgb0,0.55,0M_2= _2Q˙12 Q_12D=κ [rgb]1,0.35,0.35 [named]pgfstrokecolorrgb1,0.35,0.35D=κ: Q˙=κ(T1−T2) Q=κ(T_1-T_2)loss κlossT1 _lossT_1loss κlossT2 _lossT_2H=12c1T12+12c2T22+p122τ1+p222τ2H= [rgb]0,0,0.6 [named]pgfstrokecolorrgb0,0,0.6 12c_1T_1^2+ 12c_2T_2^2+ [rgb]0,0.55,0 [named]pgfstrokecolorrgb0,0.55,0 p_1^22 _1+ p_2^22 _2 where q=(T1,T2)q=(T_1,T_2), pi=τiT˙ip_i= _i T_iAssumptions: lumped-parameter model (spatially uniform TiT_i) ⋅· two coupled thermal masses⋅· c1=c2=1c_1=c_2=1, coupling κ=0.5κ=0.5, loss κloss=0.1 _loss=0.1 ⋅· linear Newton cooling Figure 17: Lumped heat exchange (Incropera et al., 2007). Two thermal masses with temperatures T1T_1 and T2T_2 are the generalized coordinates (ndof=2n_dof=2). Stored thermal energy Vi=12ciTi2V_i= 12c_iT_i^2 is the potential for each mass; the thermal time constants τi _i play the role of inertia; and the coupling coefficient κ provides mutual heat exchange while individual losses κloss _loss dissipate energy to the environment. This is a lumped-parameter model: temperature is assumed spatially uniform within each thermal mass (valid at low Biot number). Heat exchange (KNOWN, thermal). This environment tests whether a port-Hamiltonian learner can handle a system where “position” is temperature and “momentum” has no mechanical meaning at all (Incropera et al., 2007). Two thermal masses exchange heat with each other through a coupling coefficient κ and independently lose heat to the environment at rate κloss _loss. The key analogy: temperature differences create a thermodynamic “force” (like a compressed spring), while thermal time constants resist instantaneous temperature change (like inertia resists acceleration). The lumped-parameter assumption—each mass has a single, spatially uniform temperature—is valid when internal conduction is fast relative to surface transfer (low Biot number). • State: q=(T1,T2)q=(T_1,T_2) (temperatures of the two masses), p=(τ1T˙1,τ2T˙2)p=( _1 T_1, _2 T_2) (thermal “momenta”). The Legendre structure ∂H/∂pi=pi/τi∂ H/∂ p_i=p_i/ _i gives a thermal “velocity”—the rate of temperature change. • Potential: V(T1,T2)=12c1T12+12c2T22+12κ(T1−T2)2 [rgb]0,0,0.6 [named]pgfstrokecolorrgb0,0,0.6V(T_1,T_2)= 12c_1\,T_1^2+ 12c_2\,T_2^2+ 12κ(T_1-T_2)^2 — the first two terms store thermal energy in each mass individually; the coupling term 12κ(T1−T2)2 12κ(T_1-T_2)^2 acts like a spring connecting the two temperatures, driving heat from hot to cold whenever T1≠T2T_1≠ T_2. • Mass: M=diag(τ1,τ2)=I [rgb]0,0.55,0 [named]pgfstrokecolorrgb0,0.55,0M=diag( _1, _2)=I — thermal time constants play the role of inertia: a large τ means the mass responds sluggishly to energy input, just as a heavy object accelerates slowly under force. Set to unity in our experiments. • Damping: D=κlossI [rgb]1,0.35,0.35 [named]pgfstrokecolorrgb1,0.35,0.35D= _lossI — heat loss to the environment. This is distinct from the inter-body coupling κ: coupling redistributes energy between masses, while κloss _loss removes it from the system entirely. • Identifiability: With V and M given, κloss _loss is uniquely recoverable—the only unknown is how fast energy leaves the system. PARTIAL regime. Here the functional form of V or M is known (e.g., from first principles or CAD), but some parameters must be learned from data. Anchoring part of the model breaks gauge freedom and makes the remaining unknowns identifiable. Pendulum with unknown mass (PARTIAL, mechanical). This variant of the simple pendulum (Fig. 9) tests a core prediction of our identifiability theory: when the potential V is anchored by physics, the mass–damping gauge freedom is broken and both become recoverable. The setup is deliberately minimal—a single unknown scalar m—so that success or failure is unambiguous. • State: q=θq=θ (angle from vertical), p=mℓ2θ˙p=m ^2 θ (angular momentum). • Potential: V(θ)=−gcosθV(θ)=-g θ — given from gravity. Because mgℓmg is absorbed into a single effective parameter g, the potential shape is fully specified independent of the unknown mass. • Mass: M=mℓ2M=m ^2 — a single positive scalar to be learned. This is the simplest possible PARTIAL scenario: one unknown inertia parameter. • Damping: Learned with bounded strength. • Identifiability: Anchoring V is the key. Without it, any rescaling V→αV→α V, M→αM→α M produces identical trajectories (gauge freedom). Fixing V breaks this symmetry, making both m and D uniquely recoverable. θ1 _1m1m_1ℓ1 _1θ2 _2m2m_2end effectorggfriction Figure 18: Robot arm (2-DOF manipulator). Joint angles θ1,θ2 _1, _2 are generalized coordinates. Link inertias and gravity torques are known from CAD; joint friction (red) is learned from trajectory data. Robot arm (PARTIAL, robotics). Industrial manipulators are a compelling PARTIAL-regime test case: CAD models provide excellent link inertias and gravity torques, but joint friction—the dominant source of trajectory error in practice—must be calibrated from motion data. This is exactly the scenario our framework targets: anchor the well-understood physics (V, M) and let the learner focus on the residual dissipation. The robot arm is also the first system in our benchmark with a configuration-dependent mass matrix, where off-diagonal entries couple different joints through Coriolis and centrifugal effects. • State: q=(θ1,…,θn)q=( _1,…, _n) (joint angles), p=M(q)q˙p=M(q) q (generalized momenta). • Potential: V(q)V(q) is a gravity regressor—a linear combination of known trigonometric basis functions weighted by link-mass–length products: V(q)=∑imigℓcm,isin(⋅)V(q)= _im_ig\, _cm,i (·). The functional form comes from rigid-body kinematics; parameters may be refined from data. • Mass: M(q)M(q) is the n×n× n configuration-dependent inertia tensor, computed from CAD link inertias via the recursive Newton–Euler equations. Unlike the scalar-mass systems above, M(q)M(q) varies with configuration and has non-zero off-diagonal entries that encode how motion at one joint creates torques at another. • Damping: D(q)D(q) encodes joint friction—typically diagonal, with separate Coulomb + viscous friction per joint, learned from trajectory data. • Identifiability: With both V and M anchored by CAD, the remaining friction parameters are fully identifiable. m1m_1m2m_2m3m_3F12=−Gm1m2/r122F_12=-Gm_1m_2/r_12^2D=γI [rgb]1,0.35,0.35 [named]pgfstrokecolorrgb1,0.35,0.35D=γ Idragr1r_1r2r_2r3r_3V(q)=−∑i<jGmimjrij2+ϵsoft2 [rgb]0,0,0.6 [named]pgfstrokecolorrgb0,0,0.6V(q)=- _i<j Gm_im_j r_ij^2+ _soft^2 Port-Hamiltonian mapping V: gravitational potential (known form) M=diag(miI)M=diag(m_iI): masses (learned) D=γID=γ I: drag (learned, bounded) Regime: Partial (mim_i learned) Assumptions: 3 bodies in 2D ⋅· softened potential (ϵsoft=0.1 _soft=0.1, avoids singularity at rij→0r_ij→0) ⋅· G=1G=1, γ=0.05γ=0.05 Figure 19: N-body gravitational system (Binney & Tremaine, 2008). N=3N=3 point masses interact pairwise via Newtonian gravity. The gravitational potential V is known in functional form but individual masses mim_i are learned (PARTIAL regime). A softened potential rij2+ϵsoft2r_ij^2+ _soft^2 replaces rij2r_ij^2 to avoid the singularity at close approach. Drag D=γID=γ I models residual dissipation (e.g., tidal friction or numerical damping); it is learned with tight bounds since point-mass gravitational systems are nearly conservative. N-body gravity (PARTIAL, astrophysics). Gravitational N-body problems are the prototypical example of a system where the functional form of the physics is known (Newton’s law) but the parameters—individual masses—must be inferred from observed trajectories (Binney & Tremaine, 2008). This makes N-body an ideal PARTIAL-regime benchmark: the learner is given the gravitational template and must fill in the mass values, a task that is meaningful in astrophysics (e.g., inferring stellar masses from orbital data). We use Plummer softening to regularize the r→0r→0 singularity, a standard technique in computational astrophysics that replaces the point-particle divergence with a smooth core. • State: q=(r1,…,rN)∈ℝ2Nq=(r_1,…,r_N) ^2N (positions in 2D), p=(m1r˙1,…,mNr˙N)p=(m_1 r_1,…,m_N r_N). • Potential: V(q)=−∑i<jGmimj/rij2+ϵsoft2 [rgb]0,0,0.6 [named]pgfstrokecolorrgb0,0,0.6V(q)=- _i<jGm_im_j/ r_ij^2+ _soft^2 — the gravitational potential is known in functional form, but the mass products mimjm_im_j are unknown parameters to be learned. This is what makes the system PARTIAL rather than KNOWN: the shape of V is given, but its strength is not. • Mass: M=diag(m1I2,…,mNI2) [rgb]0,0.55,0 [named]pgfstrokecolorrgb0,0.55,0M=diag(m_1I_2,…,m_NI_2) — individual masses appear both here and in V, coupling the two components. Each mim_i is learned as a positive scalar. • Damping: D=γI [rgb]1,0.35,0.35 [named]pgfstrokecolorrgb1,0.35,0.35D=γ I — a small drag modeling tidal dissipation or dynamical friction. For isolated point masses this is nearly zero; we learn it with tight bounds to test whether PHAST can recover a weak signal. • Identifiability: The masses mim_i appear as products mimjm_im_j in V and individually in M, so they are jointly identifiable up to an overall scale. Fixing G resolves this gauge freedom. General partially known system (PARTIAL, template). This row in Table 1 represents the generic PARTIAL recipe: • Potential: V(q)=V¯(q)+εV~(q)V(q)= V(q)+ \, V(q), where V¯ V is a known physics template and V~ V is a neural residual bounded by ε . This allows the model to correct for unmodeled effects (e.g., spring nonlinearities, thermal expansion) while staying close to the physics prior. • Mass: M=m0I+UU⊤M=m_0I+U — a diagonal-plus-low-rank SPD parameterization, where m0>0m_0>0 is a known base inertia and UU⊤U captures learned corrections. • Damping: Learned with bounded strength (∑iβi≤β¯ _i _i≤ β, Eq. 11). • Identifiability: Partial — depends on how much of V and M is anchored. UNKNOWN regime. When no physics is available, all three components (V,M,D)(V,M,D) are learned from data. The port-Hamiltonian structure still guarantees passivity (dH/dt≤0dH/dt≤ 0), which stabilizes long-horizon rollouts, but the recovered parameters are no longer uniquely identifiable (see the gauge-freedom discussion in Appendix H.8). Black-box dynamics (UNKNOWN, general). When no physics is available, all three components are learned: • Potential: V(q)=fθ(q)V(q)=f_θ(q), a neural network with no structural constraint beyond smoothness. • Mass: M(q)M(q) is parameterized as an SPD neural network (diagonal-plus-low-rank, Eq. 79). • Damping: D(q)D(q) is parameterized as a PSD neural network (Eq. 10). • Identifiability: No — the inverse problem is underdetermined. Multiple (V,M,D)(V,M,D) triples can generate identical trajectories (gauge freedom, Appendix H.8). Forecasting may still be accurate, but recovered parameters lack unique physical meaning. n1n_1preyn2n_2predator+βn1n2+β n_1n_2−βn1n2-β n_1n_2growth: αn1(1−n1/K)α n_1(1-n_1/K)decay: −δn2-δ n_2ecological dissipation D(q)⪰0 [rgb]1,0.35,0.35 [named]pgfstrokecolorrgb1,0.35,0.35D(q) 0mortality μn1μ n_1mortality μn2μ n_2 Lotka–Volterra ODE n˙1=αn1(1−n1K)−βn1n2−μn1 n_1=α n_1\! (1- n_1K )-β n_1n_2-μ n_1 n˙2=δn1n2−γn2−μn2 n_2=δ n_1n_2-γ n_2-μ n_2 pH mapping (all learned) q=n1q=n_1 (prey), p=n2p=n_2 (predator) V(q)=fθ(n1)V(q)=f_θ(n_1) neural potential M≻0M 0 neural SPD scalar D⪰0D 0 neural PSD dissipation Conservative limit HLV=δn1−γlnn1H_LV=δ n_1-γ n_1 +βn2−αlnn2 H_LV=+β n_2-α n_2 PH: dH/dt≤0dH/dt≤ 0 always Assumptions: 2 species, ndof=1n_dof=1 (q=n1q=n_1, p=n2p=n_2) ⋅· K=100K=100 ⋅· α=1α=1, β=0.1β=0.1, γ=0.4γ=0.4, δ=0.1δ=0.1, μ=0.01μ=0.01 ⋅· UNKNOWN regime Figure 20: Predator–prey (Lotka–Volterra) (Murray, 2002; Strogatz, 2018). Prey density n1n_1 is the generalized coordinate q; predator density n2n_2 serves as the conjugate momentum p (ndof=1n_dof=1). All port-Hamiltonian components are learned (UNKNOWN regime): V(q)V(q) is a neural potential generalizing the classical Lotka–Volterra conserved quantity, M is a learned SPD “ecological inertia,” and D is a learned PSD dissipation encoding mortality, disease, and resource depletion. The port-Hamiltonian structure guarantees bounded dynamics even without parameter recovery. Predator–prey (UNKNOWN, ecology). This is the most challenging benchmark in our suite: an ecological system with no known physics, where the port-Hamiltonian framework must be justified purely by its structural benefits rather than by any mechanical analogy (Murray, 2002; Strogatz, 2018). The classical Lotka–Volterra model admits a conserved quantity H=δn1−γlnn1+βn2−αlnn2H=δ n_1-γ n_1+β n_2-α n_2 in the conservative limit, hinting that a Hamiltonian viewpoint is natural. Our simulator breaks conservation with a carrying-capacity term αn1(1−n1/K)α n_1(1-n_1/K) and background mortality μ, making the dynamics dissipative—exactly the setting port-Hamiltonian structure is designed for. The system has ndof=1n_dof=1: prey density n1n_1 is the generalized coordinate and predator density n2n_2 serves as the conjugate momentum. This pairing has no interpretation as mass × velocity, but it is not arbitrary: it lets the skew-symmetric J matrix couple prey growth to predator response (and vice versa), while the positive-semidefinite D captures all irreversible losses. The guarantee dH/dt≤0dH/dt≤ 0 then enforces bounded population dynamics—populations cannot diverge—even though no physics is provided. • State: q=n1q=n_1 (prey density), p=n2p=n_2 (predator density). The (q,p)(q,p) pairing is non-standard but structurally motivated: it lets the pH framework enforce that predator–prey oscillations remain bounded. • Potential: V(q)=fθ(n1) [rgb]0,0,0.6 [named]pgfstrokecolorrgb0,0,0.6V(q)=f_θ(n_1) — a neural network of prey density alone, generalizing the prey-dependent part of the classical conserved quantity. The predator-dependent terms are absorbed into the Hamiltonian’s kinetic part p2/2Mp^2/2M. • Mass: M(q) [rgb]0,0.55,0 [named]pgfstrokecolorrgb0,0.55,0M(q) — a learned SPD scalar that controls how sensitively predator density (the “momentum”) responds to changes in the ecological potential. Larger M means more sluggish predator response, analogous to heavier inertia. • Damping: D(q) [rgb]1,0.35,0.35 [named]pgfstrokecolorrgb1,0.35,0.35D(q) — a learned PSD scalar encoding irreversible ecological losses: disease, emigration, starvation, and resource depletion. • Identifiability: No — the same gauge freedom as any fully black-box system applies. But identifiability is not the goal here; the value of pH structure is stability: long-horizon rollouts remain physical (bounded populations) where unconstrained models diverge. Appendix C Architecture Details C.1 Terminology: Depth vs. Substeps vs. Blocks The word layer is overloaded in the literature; we disambiguate: • Observer depth (FD+TCN layers) refers to the number of layers in the causal observer oϕo_φ in Eq. (7). • PHAST substeps L refer to numerical substepping of the same port-Hamiltonian update ΦΔt _ t within one environment step. All substeps reuse the same learned modules (V,M,D)(V,M,D) (parameters shared); only the intermediate states differ: xt+1=ΦΔt/L∘⋯∘ΦΔt/L⏟Lsubsteps(xt).x_t+1= _ t/L ·s _ t/L_L\ substeps(x_t). (78) • PHAST blocks (optional) refer to a sequence-model wrapper (projections + residual) around the PHAST transition; our q-only experiments unroll ΦΔt _ t directly and do not stack PHAST blocks. C.2 Velocity Computation PHAST supports a general configuration-dependent mass M(q)≻0M(q) 0; in our main experiments we use a constant-mass approximation M(q)≈M(q)≈ M for efficiency. In the UNKNOWN regime, PHAST learns this (constant) mass matrix M≻0M 0 with a diagonal-plus-low-rank parameterization: M=diag(d)+∑i=1rαikiki⊤,dj>0(j=1,…,n),αi≥0.M=diag(d)+ _i=1^r _i\,k_ik_i , d_j>0\ (j=1,…,n),\ _i≥ 0. (79) Equivalently, letting Λ=diag(d) =diag(d) and U=[α1k1,…,αrkr]U=[ _1k_1,…, _rk_r], we have M=Λ+UU⊤M= +U . We compute v=M−1pv=M^-1p via the Woodbury identity (Alg. 2). The full Strang-splitting update that uses this primitive is given in Algorithm 1 (main paper). Algorithm 2 Compute v=M−1pv=M^-1p via Woodbury 0: p∈ℝnp ^n, diagonal d∈ℝ>0nd _>0^n defining Λ=diag(d) =diag(d), and U∈ℝn×rU ^n× r such that M=Λ+UU⊤M= +U 1: z←p⊘dz← p d z=Λ−1pz= ^-1p; O(n)O(n) element-wise division 2: Z←U⊘dZ← U d Z=Λ−1UZ= ^-1U; O(nr)O(nr) column-wise division 3: A←Ir+U⊤ZA← I_r+U Z O(nr2)O(nr^2) 4: c←A−1(U⊤z)c← A^-1(U z) O(r3)O(r^3) 5: v←z−Zcv← z-Zc O(nr)O(nr) 6: return v C.3 Q-Only Pipeline When only positions qtq_t are observed (common in vision-based robotics), the Markov state is xt=(qt,pt)x_t=(q_t,p_t) but ptp_t is unobserved. PHAST uses a causal observer and a canonicalizer to infer a single initial canonical state from a burn-in context, then predicts open-loop with the same PHAST core transition. Observer and takeover rollouts. We use the finite-difference + TCN observer in Eq. (7) and infer a single initial state at the end of the burn-in window, then roll out open-loop via Eq. (8). This yields a parameter-efficient q-only pipeline that is forecasting-aligned (no measurements after burn-in). 1. Canonicalizer Cψ:(q,q˙^)↦(q,p^)C_ψ:(q, q) (q, p). In our q-only experiments we use the identity canonicalizer, so p^=q˙ p= q. When M is known (and not identity), one could instead map p^=Mq˙ p=M\, q. 2. PHAST core transition Φθ:(q,p^)→(q+,p^+) _θ:(q, p)→(q^+, p^+). Same as full-state; project to q+q^+ for output. Interpretation of p p (q-only). We use p p to denote the inferred conjugate coordinate for notational consistency with the port-Hamiltonian template. With identity canonicalization, p p is a velocity-like proxy used to form an approximately Markov phase state from q-only observations (it coincides with generalized momentum only when M=IM=I up to scaling). Table 12: What is learned vs. shared in the q-only PHAST pipeline. PHAST core transition substeps L reuse the same physics components (V,M,D)(V,M,D) and differ only through intermediate states (and, if enabled, per-substep timesteps). Regime choices determine which physics components are learned (Fig. 1). Component Parameters Shared across t? Shared across substeps L? Observer oϕo_φ (FD+TCN) ϕφ ✓ — Canonicalizer CψC_ψ ψ ✓ — Physics (V,M,D)(V,M,D) θ ✓ ✓ Substep timesteps δtss=1L\δ t_s\_s=1^L τ ✓ per-substep Rollout modes. • Takeover (forecasting): infer a single initial state from burn-in, then pure integration (no measurements after burn-in). • Self-conditioned (forecasting): re-infer velocity at each step from predicted positions (errors compound; more observer-dependent). • Predict–correct (filtering): online state estimation with measurement updates (requires measurements at every step). • Feedback control: online stabilization/tracking with inputs utu_t (requires measurements; see Sec. 5). Appendix D Potential Energy Parameterizations D.1 Structured Potentials (KNOWN Regime) Cosine (pendulum-like). V(q)=∑i=1nai(1−cosqi),∇V(q)=[a1sinq1,…,ansinqn]⊤.V(q)= _i=1^na_i\,(1- q_i), ∇ V(q)=[a_1 q_1,…,a_n q_n] . (80) Quadratic (spring-like). V(q)=12q⊤Ksq,∇V(q)=Ksq.V(q)= 12\,q K_s\,q, ∇ V(q)=K_s\,q. (81) where Ks⪰0K_s 0 is a stiffness matrix. D.2 Neural Potentials (UNKNOWN Regime) MLP potential. Vη(q)=fη(q)∈ℝ,∇Vη(q)computed by automatic differentiation.V_η(q)=f_η(q) , ∇ V_η(q)\ computed by automatic differentiation. (82) In our experiments, fηf_η is a 3-layer MLP with SiLU activations. Periodic neural potential. For angular coordinates, use Fourier features: ϕ(qi)=[sin(qi),cos(qi),sin(2qi),cos(2qi),…],Vη(q)=fη(ϕ(q)).φ(q_i)=[ (q_i), (q_i), (2q_i), (2q_i),…], V_η(q)=f_η(φ(q)). (83) D.3 Hybrid Potentials (PARTIAL Regime) V(q)=V¯(q)+εV~(q),V(q)= V(q)+ \, V(q), (84) where V¯(q) V(q) is structured (given), V~(q) V(q) is neural (learned), and ε=softplus(ρε) =softplus( _ ) is a learnable scale initialized small. Appendix E Training Details E.1 Loss Function Breakdown Let y0:T−1y_0:T-1 denote the observed sequence (in q-only, yt=qty_t=q_t), where T is the sequence length. Let err(y^,y)err( y,y) denote the appropriate per-step error on the manifold (Appendix F). For full-state training (when x=(q,p)x=(q,p) is available), we use teacher-forced one-step predictions x^t+1=ΦΔt(xt) x_t+1= _ t(x_t) and define Ht=H(xt)H_t=H(x_t) and H^t+1=H(x^t+1) H_t+1=H( x_t+1). We use the following optional losses: ℒdata _data =1T−1∑t=0T−2err(y^t+1,yt+1), = 1T-1 _t=0^T-2err( y_t+1,y_t+1), (85) ℒpass _pass =1T−1∑t=0T−2max(0,H^t+1−Ht), = 1T-1 _t=0^T-2 (0, H_t+1-H_t), (86) ℒenergy _energy =1T−1∑t=0T−2|H^t+1−HtΔt+vt⊤D(qt)vt|,vt:=M−1pt, = 1T-1 _t=0^T-2 | H_t+1-H_t t+v_t D(q_t)v_t |, v_t:=M^-1p_t, (87) ℒroll _roll =t0[1Hroll∑h=1Hrollerr(y~t0+h,yt0+h)],x~0=xt0,x~h+1=ΦΔt(x~h),t0∼Unif0,…,T−1−Hroll, =E_t_0 [ 1H_roll _h=1^H_rollerr( y_t_0+h,y_t_0+h) ], x_0=x_t_0,\ x_h+1= _ t( x_h),\ t_0 \0,…,T-1-H_roll\, (88) where ℒrollL_roll evaluates open-loop rollouts (no teacher forcing). In the q-only setting, energy diagnostics can be formed using finite-difference velocity estimates as in Appendix F. E.2 Optional Rollout-Aligned Training (Pseudocode) Algorithm 3 Optional rollout-aligned training (full-state; uses Alg. 1 for ΦΔt _ t) 0: Batch of trajectories x0:T−1(b)b=1B\x^(b)_0:T-1\_b=1^B, rollout horizon HrollH_roll, weights λ∙ _ 1: ℒdata,ℒpass,ℒenergy,ℒroll←0L_data,L_pass,L_energy,L_roll← 0 2: for b=1b=1 to B do 3: Teacher-forced next-step losses 4: for t=0t=0 to T−2T-2 do 5: x^t+1←ΦΔt(xt(b)) x_t+1← _ t(x^(b)_t) Alg. 1 6: (qt,pt)←xt(b)(q_t,p_t)← x^(b)_t unpack xt(b)x^(b)_t 7: vt←M−1ptv_t← M^-1p_t 8: ℒdata+=∥x^t+1−xt+1(b)∥2L_data +=\| x_t+1-x^(b)_t+1\|^2 9: ℒpass+=max(0,H(x^t+1)−H(xt(b)))L_pass += (0,H( x_t+1)-H(x^(b)_t)) 10: ℒenergy+=|H(x^t+1)−H(xt(b))Δt+vt⊤D(qt)vt|L_energy += | H( x_t+1)-H(x^(b)_t) t+v_t D(q_t)v_t | 11: end for 12: Open-loop rollout loss 13: Choose start index t0∈0,…,T−1−Hrollt_0∈\0,…,T-1-H_roll\ 14: x~0←xt0(b) x_0← x^(b)_t_0 15: for h=1h=1 to HrollH_roll do 16: x~h←ΦΔt(x~h−1) x_h← _ t( x_h-1) Alg. 1 17: ℒroll+=∥x~h−xt0+h(b)∥2L_roll +=\| x_h-x^(b)_t_0+h\|^2 18: end for 19: end for 20: ℒ←λdataℒdata+λpassℒpass+λenergyℒenergy+λrollℒrollL← _dataL_data+ _passL_pass+ _energyL_energy+ _roll\,L_roll 21: Update parameters with ∇ℒ E.3 Hyperparameters Table 13: Hyperparameters for the q-only experiments (Sec. 4). Hyperparameter KNOWN PARTIAL UNKNOWN Damping terms r 2 2 2 Damping bound β¯ β (Windy) unbounded Δd(=0.5) d\,(=0.5) unbounded d0d_0 policy (Windy) learned fixed learned Potential hidden dim — 64 64 Potential layers — 3 3 Mass terms (if learned) — — min(4,n) (4,n) Observer FD+TCN FD+TCN FD+TCN Observer hidden dim 32 32 32 Observer layers 2 2 2 Canonicalizer identity identity identity PHAST core transition trainable yes yes yes Optimizer AdamW AdamW AdamW Learning rate 10−310^-3 10−310^-3 10−310^-3 Batch size 64 64 64 Epochs 50 50 50 Rollout horizons (eval) 10,50,100\10,50,100\ 10,50,100\10,50,100\ 10,50,100\10,50,100\ Burn-in context K 10 10 10 λdata _data 1.0 1.0 1.0 λpass _pass 0.0 0.0 0.0 λenergy _energy 0.0 0.0 0.0 λroll _roll 0.0 0.0 0.0 E.4 Baseline Hyperparameters Table 14: Baseline architecture hyperparameters (q-only experiments). All baselines are trained as causal seq2seq predictors mapping q0:t↦q^t+1q_0:t q_t+1 with shared global settings hidden dimension d=64d=64 and depth 22 (unless otherwise noted). Baseline Key settings GRU dropout=0.0 S5 dstate=64d_state=64, internal Δtssm=0.01 t_ssm=0.01 LinOSS internal Δtlinoss=1.0 t_linoss=1.0 (absorbed into learned eigenvalues) D-LinOSS internal state dim =⌊d/2⌋= d/2 , learnable per-oscillator timestep Transformer nheads=4n_heads=4, F dim =4d=4d, dropout=0.1 VPT nff_sublayers=2n_f\_sublayers=2, skew-symmetric attention Dataset parameters (Windy Pendulum). Trajectories are generated by the simulator with d0=0.3d_0=0.3, Δd=0.5 d=0.5, Δt=0.05 t=0.05, length T=200T=200, initial conditions θ0∼Unif[−π,π] _0 [-π,π], p0∼(0,42)p_0 (0,4^2), and Ntrain/Nval/Ntest=1000/200/200N_train/N_val/N_test=1000/200/200. We report mean ± std over 5 random model seeds using a fixed dataset shared across models. Initial-condition distributions (all benchmarks). All q-only benchmarks use trajectories of length T=200T=200 with Ntrain/Nval/Ntest=1000/200/200N_train/N_val/N_test=1000/200/200 (Sec. 4.1.1). Table 15 summarizes the initial-condition distributions used by the simulators. Table 15: Simulator initial conditions for the thirteen q-only benchmarks. Benchmark Δt t Initial conditions Damping Pendulum (cons) 0.05 θ0∼Unif[−π,π] _0 [-π,π], p0∼(0,32)p_0 (0,3^2) none Pendulum (damped) 0.05 θ0∼Unif[−π,π] _0 [-π,π], p0∼(0,32)p_0 (0,3^2) constant (γ=0.5γ=0.5) Pendulum (windy) 0.05 θ0∼Unif[−π,π] _0 [-π,π], p0∼(0,42)p_0 (0,4^2) d(θ)=d0+Δd|sinθ|d(θ)=d_0+ d| θ| (d0=0.3d_0=0.3, Δd=0.5 d=0.5) Cart-Pole (windy) 0.02 (x0,θ0)∼Unif([−1,1]×[−π,π])(x_0, _0) ([-1,1]×[-π,π]), p0∼Unif[−2,2]2p_0 [-2,2]^2 windy on θ (d0=0.3d_0=0.3, Δd=0.5 d=0.5) Oscillator (cons) 0.02 q0∼(0,I)q_0 (0,I), p0∼(0,I)p_0 (0,I) none Oscillator (damped) 0.02 q0∼(0,I)q_0 (0,I), p0∼(0,I)p_0 (0,I) constant (γ=0.1γ=0.1) Double Pendulum (cons) 0.01 θ1,0,θ2,0∼Unif[−π,π] _1,0, _2,0 [-π,π], ω1,0,ω2,0∼Unif[−2,2] _1,0, _2,0 [-2,2] none Double Pendulum (damped) 0.01 θ1,0,θ2,0∼Unif[−π,π] _1,0, _2,0 [-π,π], ω1,0,ω2,0∼Unif[−2,2] _1,0, _2,0 [-2,2] viscous (b=0.2b=0.2 on both joints) RLC circuit (damped) 0.02 q0∼Unif[−2,2]q_0 [-2,2], φ0∼Unif[−2,2] _0 [-2,2] resistive (R/L=0.5R/L=0.5) LJ-3 cluster (damped) 0.002 equilateral triangle at req=21/6σr_eq=2^1/6σ, perturbation σ=0.05σ=0.05, p0∼(0,0.12)p_0 (0,0.1^2) Langevin (γ=0.1γ=0.1) Heat exchange (damped) 0.02 T1,0,T2,0∼0.5+1.5⋅Unif[0,1]T_1,0,T_2,0 0.5+1.5·Unif[0,1], p0∼(0,0.32)p_0 (0,0.3^2) heat loss (κloss=0.1 _loss=0.1) N-body 3 (damped) 0.01 equilateral triangle R=2.0R=2.0, perturbation σ=0.3σ=0.3, p0∼(0,0.52)p_0 (0,0.5^2) drag (γ=0.05γ=0.05) Predator–prey (damped) 0.1 x0∼Unif[10,50]x_0 [10,50] (prey), y0∼Unif[5,20]y_0 [5,20] (predator) intra-species (α/Kα/K, μ) Appendix F Evaluation Metrics We summarize the main metrics reported in Sec. 4.1.2 and in the forecasting/identifiability tables (e.g., Tables 17 and 5). All metrics are computed on the test split under the evaluation protocol described in Sec. 4.1.1. One-step wrapped-angle MSE. In the q-only setting, the model predicts q^t+1 q_t+1 from q0:tq_0:t. We report the mean squared wrapped angular error: WrapMSE=1B(T−1)∑b=1B∑t=0T−2(wrap(q^b,t+1−qb,t+1))2,WrapMSE= 1B(T-1) _b=1^B _t=0^T-2 (wrap( q_b,t+1-q_b,t+1) )^2, (89) where wrap(⋅)wrap(·) maps angle differences to [−π,π][-π,π], B is the number of trajectories, and T is the sequence length. One-step Euclidean MSE. For Euclidean q-only environments (e.g., Oscillator), we report the standard mean squared error: MSE=1B(T−1)∑b=1B∑t=0T−2‖q^b,t+1−qb,t+1‖22.MSE= 1B(T-1) _b=1^B _t=0^T-2\| q_b,t+1-q_b,t+1\|_2^2. (90) Mixed-manifold MSE. For product manifolds with both Euclidean and angular coordinates (e.g., Cart-Pole with qt=(xt,θt)∈ℝ×1q_t=(x_t, _t) ×S^1), we report a mixed metric that averages translation MSE and wrapped-angle MSE: MixedMSE=12(MSEx+WrapMSEθ),MSEx=1B(T−1)∑b=1B∑t=0T−2(x^b,t+1−xb,t+1)2.MixedMSE= 12 (MSE_x+WrapMSE_θ ), _x= 1B(T-1) _b=1^B _t=0^T-2( x_b,t+1-x_b,t+1)^2. (91) Open-loop rollouts use the analogous metric at horizon H. Open-loop rollouts. For burn-in length K and horizon H, we condition on a ground-truth prefix q0:K−1q_0:K-1 and then predict open-loop for H steps. Baselines produce autoregressive rollouts; PHAST additionally supports takeover rollouts (burn-in then integrate). We report wrapped rollout error and (where applicable) takeover rollout error at horizon H. For Euclidean environments, we report rollout MSE at horizon H (and takeover analogues where available). For reproducibility, Table 16 maps paper metrics to the identifiers used in our implementation. Damping recovery. When a model exposes a diagonal damping prediction d^t≈D(qt) d_t≈ D(q_t), we compute MAED _D =1BT∑b=1B∑t=0T−1|d^b,t−db,t|, = 1BT _b=1^B _t=0^T-1 |\, d_b,t-d_b,t |, (92) RD2 R^2_D =1−∑b=1B∑t=0T−1(d^b,t−db,t)2∑b=1B∑t=0T−1(db,t−d¯)2+10−12, =1- _b=1^B _t=0^T-1( d_b,t-d_b,t)^2 _b=1^B _t=0^T-1(d_b,t- d)^2+10^-12, (93) where db,t=D(qb,t)d_b,t=D(q_b,t) is the ground-truth damping and d¯=1BT∑b=1B∑t=0T−1db,t d= 1BT _b=1^B _t=0^T-1d_b,t is its mean over the test set. Passivity violations on rollouts. Using an energy estimate H H formed from the predicted rollout (e.g., finite-difference velocity estimates in q-only settings), we report the fraction of rollout steps with an energy increase: PassViol(H)=1BH∑b=1B∑h=0H−1[H^b,h+1−H^b,h>ε],PassViol(H)= 1BH _b=1^B _h=0^H-1I\! [ H_b,h+1- H_b,h> ], (94) where ε is a small tolerance (we use ε=10−6 =10^-6 in code). Energy budget residual on rollouts. For our input-free benchmarks (u=0u=0), the continuous-time energy identity is dHenvdt=−Denv(q)‖v‖22, dH_envdt=-D_env(q)\,\|v\|_2^2, (95) where HenvH_env is the benchmark’s analytic energy/Hamiltonian (Appendix B) and Denv(⋅)D_env(·) is the simulator damping law (a scalar damping coefficient applied isotropically in our benchmarks). On predicted rollouts q^b,hh=0H\ q_b,h\_h=0^H, we form a finite-difference velocity v^b,h≈(q^b,h+1−q^b,h)/Δt v_b,h≈( q_b,h+1- q_b,h)/ t (wrapping angular components) and an energy estimate H^b,h H_b,h using HenvH_env (e.g., Henv(q^b,h+1,v^b,h)H_env( q_b,h+1, v_b,h) in q-only). We then report the mean absolute discrete-time residual ρb,h=H^b,h+1−H^b,hΔt−(−Denv(q^b,h)‖v^b,h‖22),EbudRes(H)=1BH∑b=1B∑h=0H−1|ρb,h|. _b,h= H_b,h+1- H_b,h t- (-D_env( q_b,h)\,\| v_b,h\|_2^2 ), (H)= 1BH _b=1^B _h=0^H-1| _b,h|. (96) Evaluating DenvD_env on the predicted configuration makes this diagnostic comparable across baselines that do not learn an explicit D(q)D(q); damping recovery metrics are reported separately when a model exposes D^(q) D(q). Table 16: Paper metric names and their corresponding metric identifiers used in our implementation. Paper metric Implementation key(s) One-step wrapped-angle MSE (WrapMSEWrapMSE) theta_wrap_mse One-step Euclidean MSE (MSEMSE) mse Rollout wrapped-angle MSE at horizon H (WrapMSEθroll(H)WrapMSE_θ^roll(H)) rollout_theta_wrap_mse_hH Rollout Euclidean MSE at horizon H (MSEroll(H)MSE^roll(H)) rollout_mse_hH Rollout mixed-manifold MSE at horizon H (MixedMSEroll(H)MixedMSE^roll(H)) rollout_mixed_mse_hH Takeover rollouts (PHAST only; Eq. (8)) rollout_takeover_* Damping recovery (RD2R^2_D, MAEDMAE_D) damping_r2, damping_mae Energy-budget residual (EbudRes(H)EbudRes(H)) rollout_energy_budget_resid_hH Passivity violation rate on rollouts (PassViol(H)PassViol(H)) rollout_passivity_violations_hH F.1 Error Sources and Diagnostics (q-only) The q-only pipeline introduces multiple amplification points beyond discretization error: (i) velocity inference from position-only measurements, (i) the semantics of the inferred canonical state (q,p^)(q, p), and (i) open-loop compounding when the model consumes its own predictions at test time. Figure 21 summarizes the computation graph and diagnostic comparisons that help localize these error sources. Teacher-forced training (next-step)q0:T−1q_0:T-1 (ground truth)observer oϕo_φq˙^0:T−1 q_0:T-1canonicalizer CψC_ψx^t=(qt,p^t) x_t=(q_t, p_t)PHAST core transition ΦΔt _ t(qt,p^t)↦(q^t+1,p^t+1)(q_t, p_t) ( q_t+1, p_t+1)loss on q^t+1 q_t+1 Observer error proxies MAE of inferred p p (when p is available) MAE of FD proxy (from q history) Train–test mismatch teacher forcing uses ground-truth history; open-loop drift can still grow. Open-loop evaluation (rollouts)burn-in contextq0:K−1q_0:K-1autoregressive rollout(observer-in-loop)takeover rollout(infer once, integrate)oϕ,Cψo_φ,C_ψ each stepΦΔt _ t onlyrollout metricsrollout metrics Observer-in-loop gap autoregressive vs. takeover rollouts (large gap ⇒ observer) Discrete-time stability passivity-violation rate and stiffness proxy for damping half-step Core/physics diagnostics rollout error and energy-budget residual Figure 21: Q-only computation graph and diagnostic probes. Training is teacher-forced (ground-truth q history), while evaluation is open-loop and therefore sensitive to compounding error. Autoregressive rollouts keep the observer in the loop; takeover rollouts isolate the PHAST core transition given a single inferred boundary state. Error sources and diagnostic comparisons. • Observer amplification (q-only): if autoregressive rollouts are much worse than takeover rollouts, repeated velocity re-inference is injecting error (noise mismatch / insufficient context). Compare inferred p p (when available) to the true p, and compare to a finite-difference proxy. • Discrete-time stiffness: energy increases or numerical instabilities can occur when the explicit damping half-step is too stiff for the chosen Δt t. Remedies include smaller Δt t, more substeps L, or bounding the damping strength (Sec. 3.2). • Damping as an error-sink (identifiability breakdown): accurate rollouts with poor damping recovery (RD2≪0R^2_D 0) suggests dissipation is compensating for errors in V/M. Bounded Householder damping reduces this degeneracy; when the isotropic base term d0d_0 is known, fixing d0d_0 further improves identifiability. • Time-rescaling ambiguity: a learned internal timestep Δtmodel≠Δt t_model≠ t can improve forecasting by reparameterizing time, but it confounds parameter recovery; for identifiability studies, we recommend fixed Δt t. • Error compounding / phase sensitivity: small one-step errors can yield large long-horizon drift. Horizon sweeps (e.g., H∈10,50,100H∈\10,50,100\) and the full-state ablation help localize whether the dominant source is the observer pipeline or the PHAST core transition. Implementation note: Table 16 lists the metric identifiers used in our evaluation pipeline; additional PHAST-specific stability diagnostics (e.g., stiffness proxies and timestep statistics) are reported when applicable. Appendix G Additional Benchmark Tables We report the full set of q-only forecasting tables across all thirteen benchmarks. All results follow the protocol in Sec. 4.1.1 (mean ± std over 5 seeds; dataset seed fixed to 42). Tables 17–19 cover the eight mechanical benchmarks (rollout MSE at horizon H=100H=100). Tables 20–24 cover the five non-mechanical systems from Table 1 (next-step MSE). Table 17: Q-only open-loop forecasting (1-DOF pendulum environments). Mean ± std of wrapped-angle rollout MSE at horizon H=100H=100 over 5 seeds. Parameter counts correspond to the Windy setting (largest; position-dependent damping introduces additional parameters). Lower is better. Model Params Conservative Damped Windy PHAST (ours) PHAST (KNOWN) 3,364 0.738±0.1430.738± 0.143 0.017±0.0050.017± 0.005 0.106±0.0200.106± 0.020 PHAST (PARTIAL) 13,736 0.680±0.0430.680± 0.043 0.018±0.0010.018± 0.001 0.092±0.0140.092± 0.014 PHAST (UNKNOWN) 13,738 0.993±0.1280.993± 0.128 0.275±0.0830.275± 0.083 0.298±0.0480.298± 0.048 Baselines GRU 37,889 2.939±0.3022.939± 0.302 1.310±0.7431.310± 0.743 1.796±0.6251.796± 0.625 S5 17,089 2.941±0.0422.941± 0.042 0.657±0.0650.657± 0.065 0.600±0.0470.600± 0.047 LinOSS 17,089 2.849±0.4752.849± 0.475 2.286±0.3242.286± 0.324 1.458±0.3241.458± 0.324 D-LinOSS 33,793 2.738±0.6302.738± 0.630 0.450±0.2410.450± 0.241 0.435±0.2390.435± 0.239 Transformer 100,161 2.320±0.2242.320± 0.224 0.493±0.1050.493± 0.105 0.824±0.1340.824± 0.134 VPT 16,833 2.875±0.1822.875± 0.182 2.111±0.2222.111± 0.222 2.218±0.1352.218± 0.135 Table 18: Q-only open-loop forecasting (2-DOF double pendulum environments). Mean ± std of wrapped-angle rollout MSE at horizon H=100H=100 over 5 seeds. Model Params Double (cons) Double (damped) PHAST (ours) PHAST (KNOWN) 3,588 0.421±0.0610.421± 0.061 0.332±0.0300.332± 0.030 PHAST (PARTIAL) 12,166 0.402±0.0470.402± 0.047 0.320±0.0320.320± 0.032 PHAST (UNKNOWN) 13,067 0.629±0.0600.629± 0.060 0.529±0.0470.529± 0.047 Baselines GRU 38,146 1.300±0.1431.300± 0.143 1.346±0.1271.346± 0.127 S5 17,218 0.618±0.0280.618± 0.028 0.630±0.0310.630± 0.031 LinOSS 17,218 1.640±0.2141.640± 0.214 1.573±0.1291.573± 0.129 D-LinOSS 33,922 1.501±0.1241.501± 0.124 1.298±0.1051.298± 0.105 Transformer 100,290 1.033±0.2481.033± 0.248 0.846±0.1320.846± 0.132 VPT 16,962 2.848±0.2132.848± 0.213 2.721±0.1562.721± 0.156 Table 19: Q-only open-loop forecasting (Oscillator environments). Mean ± std of rollout MSE at horizon H=100H=100 over 5 seeds. Model Params Oscillator (cons) Oscillator (damped) PHAST (ours) PHAST (KNOWN) 3,587 0.0010±0.00030.0010± 0.0003 0.0012±0.00040.0012± 0.0004 PHAST (PARTIAL) 12,165 0.0010±0.00020.0010± 0.0002 0.0011±0.00030.0011± 0.0003 PHAST (UNKNOWN) 12,171 0.0072±0.00310.0072± 0.0031 0.0101±0.00460.0101± 0.0046 Baselines GRU 38,146 1.8470±0.07541.8470± 0.0754 1.4565±0.08111.4565± 0.0811 S5 17,218 2.0200±0.08992.0200± 0.0899 1.7079±0.05701.7079± 0.0570 LinOSS 17,218 1.7241±0.17871.7241± 0.1787 1.3927±0.11891.3927± 0.1189 D-LinOSS 33,922 1.8603±0.27501.8603± 0.2750 1.6275±0.18091.6275± 0.1809 Transformer 100,290 1.0873±0.29871.0873± 0.2987 0.9259±0.25370.9259± 0.2537 VPT 16,962 3.0941±1.11403.0941± 1.1140 2.5710±0.79902.5710± 0.7990 Table 20: Q-only open-loop forecasting (RLC circuit, KNOWN regime). Mean ± std of next-step MSE over 5 seeds. The damped series RLC circuit (L=1L=1, C=1C=1, R=0.5R=0.5) is the electrical analog of a damped spring–mass system. Lower is better. Model Params RLC (damped) PHAST (ours) PHAST (KNOWN) 3,364 2.64×−±1.5×−2.64× 10^-5± 1.5× 10^-8 PHAST (PARTIAL) 11,878 2.64×10−5±1.1×10−82.64× 10^-5± 1.1× 10^-8 PHAST (UNKNOWN) 11,879 2.63×10−5±3.1×10−82.63× 10^-5± 3.1× 10^-8 Baselines S5 17,089 1.84×10−3±6.69×10−41.84× 10^-3± 6.69× 10^-4 LinOSS 17,089 8.57×10−4±1.26×10−48.57× 10^-4± 1.26× 10^-4 GRU 37,889 3.22×10−3±3.00×10−43.22× 10^-3± 3.00× 10^-4 LSTM 50,497 7.53×10−3±2.78×10−47.53× 10^-3± 2.78× 10^-4 VPT 16,833 8.29×10−4±4.38×10−58.29× 10^-4± 4.38× 10^-5 Transformer 100,161 4.81×10−4±1.32×10−44.81× 10^-4± 1.32× 10^-4 Table 21: Q-only open-loop forecasting (Lennard–Jones 3-particle cluster, KNOWN regime). Mean ± std of next-step MSE over 5 seeds. Three particles in 2D with pairwise LJ potential (ϵ=1ε=1, σ=1σ=1) and Langevin friction (γ=0.1γ=0.1). Lower is better. Model Params LJ-3 (damped) PHAST (ours) PHAST (KNOWN) 4,490 4.61×−±4.20×−4.61× 10^-10± 4.20× 10^-12 PHAST (PARTIAL) 13,324 4.59×10−10±2.04×10−124.59× 10^-10± 2.04× 10^-12 PHAST (UNKNOWN) 13,355 6.93×10−10±1.38×10−116.93× 10^-10± 1.38× 10^-11 Baselines S5 17,734 2.05×10−4±4.33×10−52.05× 10^-4± 4.33× 10^-5 LinOSS 17,734 3.07×10−4±6.64×10−53.07× 10^-4± 6.64× 10^-5 GRU 39,174 1.96×10−3±1.80×10−41.96× 10^-3± 1.80× 10^-4 LSTM 52,102 3.39×10−3±3.15×10−43.39× 10^-3± 3.15× 10^-4 VPT 17,478 1.49×10−3±8.07×10−41.49× 10^-3± 8.07× 10^-4 Transformer 100,806 1.20×10−3±1.62×10−41.20× 10^-3± 1.62× 10^-4 Table 22: Q-only open-loop forecasting (Heat exchange, KNOWN regime). Mean ± std of next-step MSE over 5 seeds. Two coupled thermal masses (c1=c2=1c_1=c_2=1, coupling κ=0.5κ=0.5, loss κloss=0.1 _loss=0.1). Lower is better. Model Params Heat Exch. (damped) PHAST (ours) PHAST (KNOWN) 3,592 2.42×−±9.74×−2.42× 10^-6± 9.74× 10^-10 PHAST (PARTIAL) 12,170 2.42×10−6±2.63×10−92.42× 10^-6± 2.63× 10^-9 PHAST (UNKNOWN) 12,173 2.48×10−6±1.85×10−82.48× 10^-6± 1.85× 10^-8 Baselines S5 17,218 1.00×10−3±1.62×10−41.00× 10^-3± 1.62× 10^-4 LinOSS 17,218 4.46×10−4±3.11×10−54.46× 10^-4± 3.11× 10^-5 GRU 38,146 3.63×10−3±3.31×10−43.63× 10^-3± 3.31× 10^-4 LSTM 50,818 8.03×10−3±1.39×10−48.03× 10^-3± 1.39× 10^-4 VPT 16,962 1.71×10−3±1.53×10−41.71× 10^-3± 1.53× 10^-4 Transformer 100,290 7.74×10−4±2.14×10−47.74× 10^-4± 2.14× 10^-4 Table 23: Q-only open-loop forecasting (N-body gravity, PARTIAL regime). Mean ± std of next-step MSE over 5 seeds. Three bodies in 2D with softened gravity (G=1G=1, ϵsoft=0.1 _soft=0.1) and drag (γ=0.05γ=0.05). Gravitational potential template is given; masses and damping are learned. Lower is better. Model Params N-body 3 (damped) PHAST (ours) PHAST (KNOWN) 4,489 4.28×−±1.28×−4.28× 10^-8± 1.28× 10^-11 PHAST (PARTIAL) 13,323 4.27×10−8±2.25×10−114.27× 10^-8± 2.25× 10^-11 PHAST (UNKNOWN) 13,355 4.28×10−8±6.17×10−114.28× 10^-8± 6.17× 10^-11 Baselines S5 17,734 3.03×10−3±5.98×10−43.03× 10^-3± 5.98× 10^-4 LinOSS 17,734 2.00×10−3±7.74×10−52.00× 10^-3± 7.74× 10^-5 GRU 39,174 7.18×10−3±1.79×10−37.18× 10^-3± 1.79× 10^-3 LSTM 52,102 4.64×10−2±8.17×10−34.64× 10^-2± 8.17× 10^-3 VPT 17,478 1.47×10−2±1.30×10−21.47× 10^-2± 1.30× 10^-2 Transformer 100,806 1.83×10−3±6.88×10−41.83× 10^-3± 6.88× 10^-4 Table 24: Q-only open-loop forecasting (Predator–prey, UNKNOWN regime). Mean ± std of next-step MSE over 5 seeds. Dissipative Lotka–Volterra (α=1α=1, β=0.1β=0.1, γ=0.4γ=0.4, δ=0.1δ=0.1, K=100K=100, μ=0.01μ=0.01). All PHAST components are neural (non-canonical Hamiltonian structure). Lower is better. Model Params Predator–Prey (damped) PHAST (ours) PHAST (KNOWN) 3,364 0.0203±0.00020.0203± 0.0002 PHAST (PARTIAL) 11,878 0.0203±0.00020.0203± 0.0002 PHAST (UNKNOWN) 11,879 0.0199±0.00040.0199± 0.0004 Baselines S5 17,089 6.12±0.586.12± 0.58 LinOSS 17,089 3.53±0.703.53± 0.70 GRU 37,889 4.72±0.234.72± 0.23 LSTM 50,497 5.85±0.145.85± 0.14 VPT 16,833 0.226±0.0120.226± 0.012 Transformer 100,161 0.179±0.0140.179± 0.014 Appendix H Additional Ablations H.1 Full-State Baseline (Isolating PHAST Core Transition) Motivation. Our main benchmarks are q-only and therefore include an observer+canonicalizer pipeline (Appendix C.3). To isolate the contribution of the port-Hamiltonian transition itself (PHAST core transition), we include a full-state ablation where the state (q,p)(q,p) is observed and no observer/canonicalizer is used. Setup. We evaluate on the full-state Windy Pendulum benchmark and compare PHAST variants against a dissipative Hamiltonian baseline (DHNN). All methods use the same training protocol as Sec. 4.1.1 (AdamW with a cosine schedule; 50 epochs), and we report mean ± std over 5 model seeds with a fixed dataset seed (4242). Table 25: Full-state Windy Pendulum (q,p observed): PHAST vs DHNN. Mean ± std over 5 seeds (dataset seed fixed to 42). Lower is better. DHNN does not expose an explicit damping field D(q)D(q), so RD2R^2_D/MAEDMAE_D are not reported. Model Params WrapMSEθWrapMSE_θ ↓ WrapMSEθroll(100)WrapMSE_θ^roll(100) ↓ RD2R^2_D ↑ MAEDMAE_D ↓ PHAST (PARTIAL) 10,376 1.49×10−5±2.86×10−61.49×10^-5± 2.86×10^-6 0.089±0.0050.089± 0.005 0.798±0.0340.798± 0.034 0.045±0.0050.045± 0.005 PHAST (UNKNOWN) 10,378 6.35×10−5±1.47×10−56.35×10^-5± 1.47×10^-5 0.128±0.0030.128± 0.003 0.583±0.0820.583± 0.082 0.081±0.0100.081± 0.010 DHNN 107,124 0.186±0.0610.186± 0.061 2.264±0.1602.264± 0.160 — — H.2 Configuration-Dependent Mass (Nonseparable Hamiltonian) Motivation. Several of our benchmarks (e.g., Cart-Pole and Double Pendulum; Appendix B) have a configuration-dependent inertia M(q)M(q), yielding a nonseparable Hamiltonian. Our main results intentionally use a separable constant-mass approximation (Sec. 1) to enable a lightweight leapfrog core. Here we provide a targeted ablation showing that using the true M(q)M(q) together with a nonseparable Hamiltonian integrator can materially improve long-horizon q-only rollouts. Setup. We use Windy Cart-Pole (q-only) in the KNOWN regime and compare: (i) a constant-mass approximation vs. (i) the true configuration-dependent M(q)M(q) (Eq. (61)). Both variants use Strang splitting with a Hamiltonian implicit midpoint conservative core and a fixed timestep, and follow the training protocol of Sec. 4.1.1. We report mean ± std over 3 model seeds (dataset seed fixed to 42). Table 26: Nonseparable mass ablation (Windy Cart-Pole, q-only). Mean ± std over 3 seeds (dataset seed fixed to 42). Both rows use the same implicit-midpoint conservative core; only the mass model differs. Lower is better. Model Params MixedMSEroll(H=100)MixedMSE^roll(H=100) ↓ EbudResroll(H=100)EbudRes^roll(H=100) ↓ PassViolroll(H=100)PassViol^roll(H=100) ↓ RD2R^2_D ↑ MAEDMAE_D ↓ PHAST (KNOWN, constant M) 3,589 0.096±0.0160.096± 0.016 3.31±0.393.31± 0.39 0.124±0.0070.124± 0.007 0.956±0.0040.956± 0.004 0.028±0.0010.028± 0.001 PHAST (KNOWN, M(q)M(q)) 3,589 0.031±0.0060.031± 0.006 1.30±0.061.30± 0.06 0.016±0.0070.016± 0.007 0.986±0.0020.986± 0.002 0.013±0.0010.013± 0.001 Table 27: Nonseparable mass ablation (Damped Double Pendulum, q-only). Mean ± std over 3 seeds (dataset seed fixed to 42). Both rows use the same implicit-midpoint conservative core; only the mass model differs. Lower is better. Model Params WrapMSEθroll(H=100)WrapMSE_θ^roll(H=100) ↓ EbudResroll(H=100)EbudRes^roll(H=100) ↓ PassViolroll(H=100)PassViol^roll(H=100) ↓ PHAST (KNOWN, constant M) 3,589 0.367±0.0830.367± 0.083 32.71±4.4032.71± 4.40 0.408±0.0300.408± 0.030 PHAST (KNOWN, M(q)M(q)) 3,589 0.284±0.0290.284± 0.029 17.16±1.0017.16± 1.00 0.343±0.0260.343± 0.026 H.3 PHAST Core Transition Substeps and Learnable Timestep What we vary. The PHAST core transition can compose L substeps per environment step (Eq. (78)) and can optionally learn per-substep timesteps δtss=1L\δ t_s\_s=1^L. In the main experiments (Sec. 4.1.1) we use L=1L=1 and initialize the internal timestep to the dataset sampling interval Δt t; unless otherwise noted, this timestep is learnable. Why it may help. Increasing L primarily reduces discretization error in stiff or coarsely sampled regimes; it does not introduce additional (V,M,D)(V,M,D) parameters (Table 12). Learning δtsδ t_s can improve forecasting by time-rescaling, but it can also confound identifiability since time scaling is partially interchangeable with force/damping magnitudes. Results (coarse-Δt t; fixed timestep). Table 28 summarizes a small multi-seed sweep for PHAST (PARTIAL, q-only) in two coarse-sampling settings. Increasing substeps from L=1L=1 to L=4L=4 yields modest but consistent improvements on H=100H=100 rollouts, at roughly 4×4× per-step compute. Table 28: Effect of PHAST core transition substeps at coarse sampling. Mean ± std over 5 seeds for PHAST (PARTIAL, q-only) with fixed timesteps δts=Δt/Lδ t_s= t/L. Metrics are rollout error at horizon H=100H=100. Environment Metric Δt t L=1L=1 L=4L=4 Windy Cart-Pole MixedMSEroll(H=100)MixedMSE^roll(H=100) 0.05 0.134±0.0200.134± 0.020 0.126±0.0160.126± 0.016 Double Pendulum WrapMSEθroll(H=100)WrapMSE_θ^roll(H=100) 0.04 2.412±0.0712.412± 0.071 2.366±0.0882.366± 0.088 Recommended reporting. When sweeping L and learnable δtsδ t_s, we recommend reporting both rollout error and wall-clock per step, and including at least one coarse-Δt t setting to stress discretization. H.4 Runtime Microbenchmark (PHAST Primitives) To support the efficiency claims in Sec. 1, we microbenchmark the two dominant structured primitives on CPU: (i) applying Householder damping to a vector (v↦D(q)v D(q)v, O(nr)O(nr)) and (i) computing a velocity via Woodbury (p↦M−1p M^-1p, O(nr2+r3)O(nr^2+r^3)). Figure 22 shows near-linear scaling in n for fixed rank r=2r=2. Figure 22: CPU microbenchmark of PHAST structured primitives. Average wall-clock time per call for Householder damping application and Woodbury-based mass solve, sweeping degrees of freedom n with fixed rank r=2r=2. H.5 Base Damping Policy (d0d_0) Choice in our experiments. In grey-box settings where the base damping scale is known (e.g., Windy Pendulum and Windy Cart-Pole), we fix the isotropic base term d0d_0 in the PARTIAL regime to avoid an identifiability loophole where d0d_0 becomes an error sink when ∑i=1rβi(q) _i=1^r _i(q) is bounded (Sec. 3.2). We also considered learning or bounding d0d_0, but we keep it fixed in the reported PARTIAL windy experiments for interpretability. Ablation: learnable timestep and d0d_0 policy (Windy Pendulum, q-only). Table 29 reports a small sweep that varies whether the internal timestep is fixed (Δtmodel=Δt t_model= t) or learnable, and whether the base damping term d0d_0 is fixed or learned, holding the damping-strength cap at ∑iβi(q)≤β¯ _i _i(q)≤ β with β¯=0.5 β=0.5. Learning Δtmodel t_model improves long-horizon rollouts and reduces discrete-time passivity violations, while learning d0d_0 has only a minor effect under this bound; we therefore keep d0d_0 fixed for interpretability in the reported PARTIAL windy experiments. Table 29: Learnable timestep and base damping policy (Windy Pendulum, q-only; PARTIAL). Mean ± std over 3 seeds on CPU (20 epochs; dataset seed fixed to 42). All runs use Householder damping with a fixed cap ∑iβi(q)≤β¯ _i _i(q)≤ β (β¯=0.5 β=0.5) and report horizon-H=100H=100 diagnostics under the standard q-only rollout protocol (Sec. 4.1.1). Setting WrapMSEθroll(H=100)WrapMSE_θ^roll(H=100) ↓ RD2R^2_D ↑ MAEDMAE_D ↓ EbudRes(H=100)EbudRes(H=100) ↓ PassViolroll(H=100)PassViol^roll(H=100) ↓ fixed Δtmodel t_model, fixed d0d_0 0.175±0.0090.175± 0.009 0.465±0.0380.465± 0.038 0.084±0.0030.084± 0.003 1.34±0.041.34± 0.04 0.020±0.0040.020± 0.004 learnable Δtmodel t_model, fixed d0d_0 0.153±0.0130.153± 0.013 0.472±0.0400.472± 0.040 0.085±0.0030.085± 0.003 1.28±0.041.28± 0.04 0.015±0.0040.015± 0.004 fixed Δtmodel t_model, learned d0d_0 0.175±0.0090.175± 0.009 0.476±0.0470.476± 0.047 0.082±0.0040.082± 0.004 1.34±0.041.34± 0.04 0.020±0.0040.020± 0.004 learnable Δtmodel t_model, learned d0d_0 0.153±0.0130.153± 0.013 0.478±0.0470.478± 0.047 0.084±0.0040.084± 0.004 1.28±0.041.28± 0.04 0.015±0.0040.015± 0.004 H.6 Baseline Capacity Scaling (Q-Only) What we vary. To address concerns that large margins might be driven by baseline under-capacity, we run a small capacity scaling check on Windy Pendulum by doubling the baseline hidden dimension (from 64 to 128) while keeping depth fixed (2 layers) and using the same training protocol as Sec. 4.1.1. Table 30: Baseline capacity scaling (Windy Pendulum, q-only). Mean ± std of wrapped-angle rollout MSE at horizon H=100H=100 over 3 seeds (dataset seed fixed to 42) on CPU. Increasing baseline hidden dimension does not close the gap to PHAST (Table 17). Model (scaled) Params WrapMSEθroll(H=100)WrapMSE_θ^roll(H=100) ↓ GRU (d=128d=128) 149,505 2.130±0.5662.130± 0.566 S5 (d=128d=128) 34,049 0.654±0.0150.654± 0.015 Transformer (d=128d=128) 396,929 0.868±0.2370.868± 0.237 H.7 Bounding Damping Strength What we vary. To study the impact of spectral control on identifiability (Sec. 3.2), we vary a single scalar bound on the total damping strength, ∑i=1rβi(q)≤β¯ _i=1^r _i(q)≤ β (Eq. 11). We also consider per-term bounds βi(q)≤βmax _i(q)≤ _ . d0Id_0I++k1k_1β1k1k1⊤ _1k_1k_1 ++k2k_2β2k2k2⊤ _2k_2k_2 ==D(q)D(q) (PSD)Why bound ∑i=1rβi≤β¯ _i=1^r _i≤ β ?Unboundedtrue D absorbs model errorR2<0R^2<0 (non-identifiable)⟶ (ours)β¯ βV, M must be correctR2R^2 improves ✓ Figure 23: Low-rank PSD damping with spectral bound. Top: D(q)=d0I+∑i=1rβi(q)ki(q)ki(q)⊤D(q)=d_0I+ _i=1^r _i(q)\,k_i(q)k_i(q) is PSD by construction. Bottom: Bounding ∑i=1rβi(q)≤β¯ _i=1^r _i(q)≤ β certifies λmax(D(q))≤d0+β¯ _ (D(q))≤ d_0+ β (Eq. (11)). Ablation: cap value sensitivity (Windy Pendulum, q-only; PARTIAL). Table 31 sweeps the total-strength cap ∑iβi(q)≤β¯ _i _i(q)≤ β at a reduced training budget with fixed Δtmodel=Δt t_model= t and fixed d0d_0. We observe a clear identifiability sweet spot around β¯=0.5 β=0.5: loosening the cap (unbounded or β¯=1.0 β=1.0) makes damping recovery fail (RD2≪0R^2_D 0), while tightening to β¯=0.25 β=0.25 slightly increases rollout error but improves energy-residual and passivity diagnostics. Table 31: Cap sensitivity sweep (Windy Pendulum, q-only; PARTIAL). Mean ± std over 3 seeds on CPU (20 epochs; dataset seed fixed to 42). All runs use Householder damping with fixed Δt=0.05 t=0.05, fixed Δtmodel=Δt t_model= t, fixed d0d_0, and vary only the total-strength cap ∑iβi(q)≤β¯ _i _i(q)≤ β. We report horizon-H=100H=100 diagnostics under the standard q-only rollout protocol (Sec. 4.1.1). Cap on ∑i=1rβi(q) _i=1^r _i(q) WrapMSEθroll(H=100)WrapMSE_θ^roll(H=100) ↓ RD2R^2_D ↑ MAEDMAE_D ↓ EbudRes(H=100)EbudRes(H=100) ↓ PassViolroll(H=100)PassViol^roll(H=100) ↓ Unbounded 0.162±0.0150.162± 0.015 −40.940±8.630-40.940± 8.630 0.897±0.0890.897± 0.089 1.43±0.041.43± 0.04 0.022±0.0060.022± 0.006 ∑i=1rβi(q)≤0.25 _i=1^r _i(q)≤ 0.25 0.180±0.0120.180± 0.012 −0.080±0.081-0.080± 0.081 0.105±0.0040.105± 0.004 1.32±0.031.32± 0.03 0.019±0.0040.019± 0.004 ∑i=1rβi(q)≤0.5 _i=1^r _i(q)≤ 0.5 0.175±0.0090.175± 0.009 0.465±0.0380.465± 0.038 0.084±0.0030.084± 0.003 1.34±0.041.34± 0.04 0.020±0.0040.020± 0.004 ∑i=1rβi(q)≤1.0 _i=1^r _i(q)≤ 1.0 0.175±0.0140.175± 0.014 −1.506±0.582-1.506± 0.582 0.211±0.0260.211± 0.026 1.38±0.051.38± 0.05 0.022±0.0040.022± 0.004 Table 32: UNKNOWN regime: effect of bounding total damping strength (Windy Pendulum, q-only). Mean ± std over 5 seeds on CPU (same data protocol as Sec. 4.1.1). Bounding ∑i=1rβi(q) _i=1^r _i(q) substantially improves damping recovery (identifiability) but can degrade long-horizon open-loop forecasting, illustrating a forecasting–identifiability trade-off in the UNKNOWN regime. Cap on ∑i=1rβi(q) _i=1^r _i(q) WrapMSEθroll(H=100)WrapMSE_θ^roll(H=100) ↓ RD2R^2_D ↑ MAEDMAE_D ↓ EbudRes(H=100)EbudRes(H=100) ↓ Unbounded 0.298±0.0540.298± 0.054 −96.546±17.774-96.546± 17.774 1.343±0.1291.343± 0.129 2.581±0.2142.581± 0.214 ∑i=1rβi(q)≤0.5 _i=1^r _i(q)≤ 0.5 0.416±0.0710.416± 0.071 −0.234±0.268-0.234± 0.268 0.134±0.0190.134± 0.019 2.713±0.2072.713± 0.207 ∑i=1rβi(q)≤1.0 _i=1^r _i(q)≤ 1.0 0.396±0.0760.396± 0.076 −6.457±1.376-6.457± 1.376 0.363±0.0340.363± 0.034 2.666±0.2142.666± 0.214 H.8 Gauge Freedom and Parameter Identification Why identifiability can fail. Even in noise-free settings, recovering physical parameters from trajectories can be ill-posed: the map from a mechanical model (M(⋅),V(⋅))(M(·),V(·)) to observed trajectories need not be injective. Two common sources of non-uniqueness are: (i) Lagrangian gauge—adding a total time derivative to the Lagrangian does not change the Euler–Lagrange equations, L′(q,q˙,t)=L(q,q˙,t)+dtF(q,t),L (q, q,t)=L(q, q,t)+ ddtF(q,t), (97) and (i) coordinate/canonical freedom—changes of variables can alter the representation of (M,V)(M,V) (and the meaning of p) without changing the underlying configuration-space trajectories. In data-driven learning, this appears as parameter trade-offs: distinct (M,V)(M,V) (and, in dissipative settings, also D) can yield similarly accurate rollouts. Connection to PHAST regimes. PHAST’s knowledge regimes can be viewed as gauge-fixing choices that trade flexibility for identifiability: KNOWN anchors (V,M)(V,M) and learns only D(q)D(q); PARTIAL anchors the form of V and uses calibrated damping bounds; UNKNOWN is intentionally flexible but can be non-identifiable (Sec. 4.1.2). The damping-strength bound in Eq. (11) plays an analogous role to gauge fixing: it prevents dissipation from absorbing model mismatch and improves physical recovery (Tables 6–32). Illustrative controlled study (double pendulum). To make the gauge-freedom issue concrete, we include a small controlled experiment on a conservative double pendulum. We compare learning M(q)M(q) and/or V(q)V(q) under different anchors using an Euler–Lagrange residual loss: r(q,q˙,q¨):=M(q)q¨+CM(q,q˙)q˙+∇V(q),ℒres=‖r(q,q˙,q¨)‖22,r(q, q, q):=M(q) q+C_M(q, q) q+∇ V(q), _res=\|r(q, q, q)\|_2^2, (98) where CM(q,q˙)q˙C_M(q, q) q is the Coriolis/centrifugal term induced by M(q)M(q). In components, [CM(q,q˙)q˙]i=∑j,k=1nΓijk(q)q˙jq˙k,Γijk(q)=12(∂qkMij+∂qjMik−∂qiMjk). [C_M(q, q) q ]_i= _j,k=1^n _ijk(q) q_j q_k, _ijk(q)= 12 ( _q_kM_ij+ _q_jM_ik- _q_iM_jk ). (99) The residual loss directly constrains M when V is anchored. Empirically, when both M and V are fully free (neural V), parameter recovery fails despite low trajectory MSE, consistent with gauge freedom. Anchoring V (exactly, or via a structured form with learnable parameters) makes M identifiable and yields accurate recovery. We report results for two mass parameterizations: (i) a full SPD (Cholesky) model used in earlier canonicalization studies, and (i) the current PHAST UNKNOWN-style mass NeuralMass (diagonal ++ low-rank outer products with Woodbury solves). Double pendulum physics. The true configuration-dependent mass and potential for the double pendulum (unit masses and lengths) are: M(q)=(2cos(θ1−θ2)cos(θ1−θ2)1),V(q)=−2gcosθ1−gcosθ2,M(q)= pmatrix2& ( _1- _2)\\ ( _1- _2)&1 pmatrix, V(q)=-2g _1-g _2, (100) where g=9.81g=9.81. The Euler–Lagrange residual loss enforces M(q)q¨+CM(q,q˙)q˙+∇V(q)=0M(q) q+C_M(q, q) q+∇ V(q)=0 via the Coriolis term CMC_M. Structured potential (PARTIAL analogue). “Structured V” means the functional form is fixed to the true physics, but coefficients are learned: V(q;g1,g2)=−g1cosθ1−g2cosθ2,V(q;g_1,g_2)=-g_1 _1-g_2 _2, (101) where g1,g2g_1,g_2 are trainable parameters (true values: g1=2g=19.62g_1=2g=19.62, g2=g=9.81g_2=g=9.81). This is analogous to PHAST’s PARTIAL regime where the potential template is known but parameters are learned. NeuralMass parameterization. The NeuralMass module uses the same Householder-style low-rank form as PHAST damping (Sec. 3.2): M(q)=diag(d(q))+∑i=1rαi(q)ki(q)ki(q)⊤,dj>0,αi≥0,‖ki‖=1,M(q)=diag(d(q))+ _i=1^r _i(q)\,k_i(q)k_i(q) , d_j>0,\ _i≥ 0,\ \|k_i\|=1, (102) where positivity is enforced via softplus and directions are normalized. Inversion uses the Woodbury identity in O(nr2+r3)O(nr^2+r^3). The SPD (Cholesky) baseline parameterizes M=LL⊤M=L with lower-triangular L having positive diagonal. Table 33: Gauge-freedom ablation (conservative double pendulum). Mean over 2 runs (50 epochs). Lower is better. “Structured V” means the functional form is fixed but coefficients are learned (analogous to PHAST PARTIAL). Errors are averaged over held-out configurations: FrobErrM:=q‖M(q)−M^(q)‖FFrobErr_M:=E_q\|M(q)- M(q)\|_F and MAE∇V:=1nq‖∇V(q)−∇V^(q)‖1MAE_∇ V:= 1nE_q\|∇ V(q)-∇ V(q)\|_1. Setup Mass FrobErrMFrobErr_M ↓ MAE∇VMAE_∇ V ↓ Conclusion Learn M (known V) + ℒresL_res SPD 0.0110.011 — Identifiable Learn M (known V) + ℒresL_res NeuralMass 0.3060.306 — Identifiable (high var.) Learn V (known M) + ℒresL_res Known M — 0.0940.094 Identifiable Learn both (M + neural V) + ℒresL_res SPD 2.372.37 6.936.93 Non-identifiable Learn both (M + neural V) + ℒresL_res NeuralMass 2.232.23 6.736.73 Non-identifiable Learn M + structured V + ℒresL_res SPD 0.0350.035 0.0820.082 Identifiable Learn M + structured V + ℒresL_res NeuralMass 0.0890.089 0.0660.066 Identifiable Learn both (M + sparse-basis V) + ℒresL_res NeuralMass 1.791.79 12.0512.05 Non-identifiable Appendix I Notation Reference Configuration space. We write q∈q for the configuration, where Q is typically a product manifold =ℝne×(1)naQ=R^n_e×(S^1)^n_a with ne+na=n_e+n_a=n. In implementation, q is stored in ℝnR^n and angular errors use wrap(⋅)wrap(·) to respect periodicity (Appendix F). We use x to denote the phase-space state x=(q,p)x=(q,p), but in the Cart-Pole environment q=(x,θ)q=(x,θ) also uses x for the cart position coordinate; context disambiguates. Table 34: Index ranges. Index Range Meaning b 1,…,B1,…,B Trajectory index (batching) t 0,…,T−10,…,T-1 Time index (transitions use t=0,…,T−2t=0,…,T-2) h 0,…,H−10,…,H-1 Rollout step index (horizon H) s 1,…,L1,…,L PHAST substep index i 1,…,r1,…,r Rank-1 term index in low-rank updates j 1,…,n1,…,n Coordinate index Table 35: Symbol reference. Symbols are grouped by the section in which they are introduced. Gray rows mark group headers indicating the originating section. Symbol Space / shape Meaning Problem setting & knowledge regimes Sec. 1 n ℕN Degrees of freedom Q manifold Configuration space (e.g., ℝne×(1)naR^n_e×(S^1)^n_a) q Q Generalized configuration (forecast target) p ℝnR^n Conjugate momentum (latent in q-only setting) x=(q,p)x=(q,p) ×ℝnQ×R^n Phase-space state yty_t Q Observation (q-only: yt=qty_t=q_t) Port-Hamiltonian dynamics Sec. 3.1 H(q,p)H(q,p) ×ℝn→ℝQ× ^n→ Hamiltonian (total energy) V(q)V(q) →ℝQ Potential energy M(q)M(q) →ℝn×nQ ^n× n Mass/inertia matrix (SPD; general, configuration-dependent) M ℝn×nR^n× n Constant-mass approximation M(q)≈M(q)≈ M (used in main experiments) v=M(q)−1pv=M(q)^-1p ℝnR^n Generalized velocity D(q)D(q) →ℝn×nQ ^n× n Damping matrix (PSD) J ℝ2n×2nR^2n× 2n Interconnection matrix (skew-symmetric) R ℝ2n×2nR^2n× 2n Dissipation matrix (PSD; block-diagonal with D) I ℝn×nR^n× n Identity matrix (size inferred from context) m ℕN Input/port dimension (when control is present) u ℝmR^m Port input (control) G ℝ2n×mR^2n× m Input matrix (port mapping) yporty^port ℝmR^m Port output yport=G⊤∇H(x)y^port=G ∇ H(x) Low-rank parameterizations Sec. 3.2 r ℕN Rank of low-rank expansion d0d_0 ℝ≥0R_≥ 0 Isotropic baseline damping βi(q) _i(q) →ℝ≥0Q _≥ 0 Strength of rank-11 damping term ki(q)k_i(q) →n−1Q ^n-1 Direction of rank-11 damping term (‖ki(q)‖2=1\|k_i(q)\|_2=1) β¯ β ℝ≥0R_≥ 0 Damping strength bound (∑iβi≤β¯ _i _i≤ β) Λ ℝn×nR^n× n Diagonal base matrix Λ=diag(d) =diag(d) (Woodbury) U ℝn×rR^n× r Low-rank factor matrix (Woodbury: M=Λ+UU⊤M= +U ) IrI_r ℝr×rR^r× r Identity matrix in rank-r Woodbury updates Structure-preserving integration Sec. 3.3 Δt t ℝ>0R_>0 Time step ΦΔt _ t ×ℝn→×ℝnQ× ^n→ × ^n Full integrator map (Strang splitting) ΦHΔt _H t ×ℝn→×ℝnQ× ^n→ × ^n Conservative step (symplectic core) ΦDΔt/2 _D t/2 ×ℝn→×ℝnQ× ^n→ × ^n Dissipation half-step L ℕN PHAST core transition substeps per step δtsδ t_s ℝ>0R_>0 Per-substep timestep (optional; s=1,…,Ls=1,…,L) Q-only observer pipeline Eqs. (6)–(8) q˙tfd q_t^fd ℝnR^n Finite-difference velocity estimate from q-only data oϕo_φ causal map Velocity observer (FD+TCN): (q0:t,q˙0:tfd)↦δt(q_0:t, q^fd_0:t) _t δt _t ℝnR^n Observer correction to finite-difference velocity q˙^t q_t ℝnR^n Observer-estimated velocity CψC_ψ ×ℝn→×ℝnQ× ^n→ × ^n Canonicalizer: (q,q˙^)↦(q,p^)(q, q) (q, p) p^t p_t ℝnR^n Canonicalizer output (identity: p^t=q˙^t p_t= q_t) Πq(q,p) _q(q,p) ×ℝn→Q× ^n→ Position projection Πq(q,p)=q _q(q,p)=q Experiments & evaluation Sec. 4.1.1–4.1.2 T ℕN Trajectory length (sequence length) K ℕN Burn-in context length (q-only rollouts) H ℕN Rollout horizon B ℕN Number of trajectories used in metric averaging Denv(q)D_env(q) →ℝ≥0Q _≥ 0 Simulator (ground-truth) damping coefficient wrap(⋅)wrap(·) ℝ→[−π,π]R→[-π,π] Wrap angular differences mod 2π2π