Poisson-Dirac Neural Networks for Modeling Coupled Dynamical Systems across Domains

## Abstract

Abstract Deep learning has achieved great success in modeling dynamical systems, providing data-driven simulators to predict complex phenomena, even without known governing equations.
However, existing models have two major limitations: their narrow focus on mechanical systems and their tendency to treat systems as monolithic.
These limitations reduce their applicability to dynamical systems in other domains, such as electrical and hydraulic systems, and to coupled systems.
To address these limitations, we propose Poisson-Dirac Neural Networks (PoDiNNs), a novel framework based on the Dirac structure that unifies the port-Hamiltonian and Poisson formulations from geometric mechanics.
This framework enables a unified representation of various dynamical systems across multiple domains as well as their interactions and degeneracies arising from couplings.
Our experiments demonstrate that PoDiNNs offer improved accuracy and interpretability in modeling unknown coupled dynamical systems from data.

## 1 Introduction

Deep learning has achieved great success in modeling dynamical systems (Chen et al., ; Anandkumar et al., ), following its successes in image processing and natural language processing (He et al., ; Vaswani et al., ).
It provides a data-driven approach for predicting the behavior of complex dynamical systems, even when governing equations are unknown (Kasim & Lim, ; Matsubara & Yaguchi, ).
These models function as computational simulators and show promise in applications across diverse fields, including weather forecasting, mechanical design, and system control (Lam et al., ; Pfaff et al., ; Horie et al., ).

However, Greydanus et al. () identified a key limitation of data-driven models: they accumulate modeling errors in long-term predictions, leading to rapid failure.
Hamiltonian Neural Networks (HNNs) were proposed to overcome this limitation by incorporating Hamiltonian mechanics.
Inspired by this approach, many studies have developed neural network models that not only learn the dynamics from data, but also adhere to fundamental laws of physics.
Examples include Lagrangian Neural Networks (LNNs) (Cranmer et al., ), Neural Symplectic Forms (NSFs) (Chen et al., ), Constrained HNNs (CHNNs) (Finzi et al., ), and Dissipative SymODENs (Zhong et al., ), as shown in Table 1.
We emphasize that our focus is on modeling unknown dynamical systems from data, which differs from solving known symbolic equations, such as in Physics Informed Neural Networks (Sirignano & Spiliopoulos, ; Raissi et al., ; Du & Zaki, ).

Despite recent progress, two key limitations remain in modeling dynamical systems, especially those described by ordinary differential equations (ODEs).
The first limitation is the narrow focus on mechanical systems.
Models applied to other domains, such as electric circuits or magnetic fields, often fail to leverage the governing principles of those systems, such as Kirchhoff’s current and voltage laws (Jin et al., ; Matsubara & Yaguchi, ).
The second limitation is that most methods treat the system as a single, monolithic entity.
In reality, many systems consist of interacting components, such as robot arms with multiple joints or electric circuits with various elements (Yoshimura & Marsden, ).
Although the port-Hamiltonian formulation theoretically addresses these interactions, no prior work has fully leveraged its potential.
These limitations prevent current methods from effectively handling multiphysics scenarios, where systems from different domains interact, such as those involving through DC motors (van der Schaft & Jeltsema, ; Gay-Balmaz & Yoshimura, ).

Figure: Figure 1: Conceptual diagram of PoDiNNs.
Refer to caption: x1.png

To address these limitations, we propose *Poisson-Dirac Neural Networks (PoDiNNs)*, which leverage the Dirac structure to unify the port-Hamiltonian and Poisson formulations (Courant, ; Duindam et al., ; van der Schaft, ; van der Schaft & Jeltsema, ; Yoshimura & Marsden, ).
The Dirac structure explicitly represents the coupling between internal and external components.
A conceptual diagram is shown in Fig. 1.
The advantages of PoDiNNs are summarized below and in Table 1, and are validated through experiments on seven simulation datasets spanning mechanical, rotational, electro-magnetic, and hydraulic domains.

Identifying Coupling Patterns Many real-world systems are coupled systems, consisting of interacting components (Yoshimura & Marsden, ).
While existing methods represent these as single vector fields or energy functions, PoDiNNs explicitly learn internal couplings as bivectors on vector bundles, separate from component-wise characteristics modeled by neural networks.
This approach improves both interpretability and modeling accuracy.

Applicable to Various Domains of Systems The Dirac structure allows PoDiNNs to describe various dynamical systems across multiple domains.
While most deep learning models of ODEs focus on mechanical systems, PoDiNNs can model interactions of mechanical, rotational, electrical, and hydraulic systems, addressing real-world multiphysics scenarios.

Unifying Various Aspects of Dynamical Systems Previous models have addressed specific aspects of dynamical systems, such as degenerate dynamics (Finzi et al., ; Jin et al., ), energy dissipation, external inputs (Zhong et al., ), and coordinate transformations (Chen et al., ; Jin et al., ).
PoDiNNs offer a unified framework to address all these aspects together by leveraging the Dirac structure.

**Table 1: Comparison of Methods for Modeling Dynamical Systems**
| Model | Formulation | Identifying<br>Coupling Patterns | Applicable to<br>Multiphysics | Degeneracy | Dissipation | External Inputs | Coordinate-Free |  |
| --- | --- | --- | --- | --- | --- | --- | --- | --- |
| HNNs (Greydanus et al., ) | canonical Hamiltonian | ✗ | ✗ | ✗ | ✗ | ✗ | ✗ |  |
| LNNs (Cranmer et al., ) | Lagrangian | ✗ | ✗ | ✗ | ✗ | ✗ | ✗ |  |
| NSFs (Chen et al., ) | general Hamiltonian | ✗ | ✗ | ✗ | ✗ | ✗ | ✓ |  |
| CHNNs (Finzi et al., ) | constrained canonical Hamiltonian | ✗ | ✗ | ✻ | ✗ | ✗ | ✗ |  |
| CLNNs (Finzi et al., ) | constrained Lagrangian | ✗ | ✗ | ✻ | ✗ | ✗ | ✗ |  |
| PNNs (Jin et al., ) | Poisson | ✗ | ✗ | ✓ | ✗ | ✗ | ✓ |  |
| Dis. SymODENs (Zhong et al., ) | port-Hamiltonian on Darboux coordinates | ✗ | ✗ | ✗ | ✓ | $\dagger$ | ✗ |  |
| PoDiNNs (proposed) | poisson-Dirac with ports | ✓ | ✓ | ✓ | ✓ | ✓ | ✓ |  |

## 2 Background Theory and Related Work

Attempts to learn ODEs have a long history (Nelles, ), but Neural Ordinary Differential Equations (Neural ODEs) transformed this field (Chen et al., ).
Neural ODEs train networks to approximate the vector field (the right-hand side of the ODE), which is then integrated by a numerical integrator.
However, these models are often overly general and fail to capture dynamical systems governed by principles such as energy conservation (Greydanus et al., ).
To address this, many refinements have been proposed by integrating insights from analytical mechanics.

### 2.1 Hamiltonian Systems

Canonical Hamiltonian Systems The Hamiltonian formulation of mechanics defines Hamilton’s equations of motion, using a pair of $n$-dimensional generalized coordinates ${\bm{q}}=(q_{1},\dots,q_{n})$ and momenta ${\bm{p}}=(p_{1},\dots,p_{n})$ as the state.
Given an energy function $H$, the system dynamics evolves as:

$$ $\textstyle\dot{q}_{i}=\frac{\partial H}{\partial p_{i}},\quad\dot{p}_{i}=- \frac{\partial H}{\partial q_{i}},\text{ or equivalently, }\begin{pmatrix}\dot {\bm{q}}\\ \dot{\bm{p}}\end{pmatrix}=\begin{pmatrix}O&I\\ -I&O\end{pmatrix}\nabla H({\bm{q}},{\bm{p}})$ (1) $$

where $\dot{\ }$ denotes the time derivative, $I$ and $O$ are the $n$-dimensional identity and zero matrices, respectively, and $\nabla H$ denotes the gradient of $H$.
Systems describable by Hamilton’s equations are Hamiltonian systems.
See Appendix A.1 for an example.

HNNs were introduced to learn such systems, with a neural network $H_{N\!N}$ serving as the energy function $H$ (Greydanus et al., ).
HNNs have shown greater long-term robustness than Neural ODEs (Chen et al., ; Gruver et al., ).

Hamiltonian Systems of Coordinate-Free Form From the differential geometry perspective, the equations of motion in Eq. (1) are defined on the cotangent bundle ${\mathcal{M}}={T^{*}{\mathcal{Q}}}$ of the configuration space ${\mathcal{Q}}=\mathbb{R}^{n}$ with the standard symplectic form $\Omega=\sum_{i}{\mathrm{d}}q_{i}\wedge{\mathrm{d}}p_{i}$, where ${\bm{q}}\in{\mathcal{Q}}$ and ${\bm{p}}\in{T^{*}_{\bm{q}}{\mathcal{Q}}}$, respectively.
Here, ${T^{*}_{\bm{q}}{\mathcal{Q}}}$ denotes the cotangent space to ${\mathcal{Q}}$ at ${\bm{q}}$, and $\wedge$ denotes the exterior product.
The symplectic form $\Omega$ leads to a skew-symmetric linear bundle map $\Omega^{\flat}_{\bm{u}}:{T_{\bm{u}}{\mathcal{M}}}\to{T^{*}_{\bm{u}}{\mathcal{M
}}}$ at each point ${\bm{u}}\in{\mathcal{M}}$, which satisfies $\langle\Omega^{\flat}_{\bm{u}}({\bm{v}}),{\bm{w}}\rangle=\Omega_{\bm{u}}({\bm{
v}},{\bm{w}})$ for any ${\bm{u}}\in{\mathcal{M}}$ and ${\bm{v}},{\bm{w}}\in{T_{\bm{u}}{\mathcal{M}}}$, where $\langle\cdot,\cdot\rangle$ denotes the natural pairing.
We will denote by the subscript u an assignment at point ${\bm{u}}$ and omit the subscripts if an equation holds independently of the point ${\bm{u}}$.
Given a smooth function $H:{\mathcal{M}}\to\mathbb{R}$, its differential is denoted by ${\mathrm{d}}H$.
A vector field $X$ on ${\mathcal{M}}$ assigns a tangent vector $X_{\bm{u}}\in{T_{\bm{u}}{\mathcal{M}}}$ to each point ${\bm{u}}\in{\mathcal{M}}$.
Then, Hamilton’s equations in the coordinate-free form defines a vector field $X_{H}$ called the Hamiltonian vector field by

$$ $\Omega^{\flat}(X_{H})={\mathrm{d}}H.$ (2) $$

$X_{H}$ defines the time evolution of the state ${\bm{u}}$ as $\dot{\bm{u}}(t)=(X_{H})_{{\bm{u}}(t)}$ at time $t$, and ${\bm{u}}(t)$ over a certain period is a solution of the system.
Then, the tuple $(H,\Omega,X_{H})$ is called a Hamiltonian system.
A Hamiltonian system with the standard symplectic form $\Omega$ is said to be canonical, and its coordinate system is known as Darboux coordinates.

Hamilton’s equations also describe the same dynamics using generalized velocities ${\bm{v}}\in{T_{\bm{q}}{\mathcal{Q}}}$ instead of generalized momenta ${\bm{p}}\in{T^{*}_{\bm{q}}{\mathcal{Q}}}$, i.e., on the tangent bundle ${T{\mathcal{Q}}}$ rather than the cotangent bundle ${T^{*}{\mathcal{Q}}}$.
In this case, the symplectic form $\Omega$ is replaced by a Lagrangian 2-form (Marsden & Ratiu, ).
From another viewpoint, the symplectic form $\Omega$ defines the coordinate system.
Lagrangian Neural Networks (LNNs) learn the dynamics on the tangent bundle ${T{\mathcal{Q}}}$ (Cranmer et al., ).
Neural Symplectic Forms (NSFs) learn the symplectic form directly from data, generalizing HNNs and LNNs for arbitrary coordinate systems (Chen et al., ).
In general, Hamiltonian systems preserve the energy $H$, as ${\mathcal{L}}_{X_{H}}H=\langle{\mathrm{d}}H,X_{H}\rangle=\Omega(X_{H},X_{H})=0$, where ${\mathcal{L}}_{X_{H}}$ denotes the Lie derivative along $X_{H}$, and the last equality follows from the skew-symmetry of $\Omega$.

Poisson Systems Because the symplectic form $\Omega$ is non-degenerate in the sense that the bundle map $\Omega^{\flat}_{\bm{u}}$ is non-degenerate, it leads to a 2-tensor $B$ called a Poisson bivector satisfying $B({\mathrm{d}}H,{\mathrm{d}}G)=\Omega(X_{H},X_{G})$ for any smooth functions $H,G$ on ${\mathcal{M}}$.
$B$ leads to a skew-symmetric linear bundle map $B^{\sharp}_{\bm{u}}:{T^{*}_{\bm{u}}{\mathcal{M}}}\to{T_{\bm{u}}{\mathcal{M}}}$ at each point ${\bm{u}}\in{\mathcal{M}}$, which satisfies $B_{\bm{u}}({\bm{\alpha}},{\bm{\beta}})=\langle{\bm{\alpha}},B^{\sharp}_{\bm{u}
}({\bm{\beta}})\rangle$ for any ${\bm{u}}\in{\mathcal{M}}$ and ${\bm{\alpha}},{\bm{\beta}}\in{T^{*}_{\bm{u}}{\mathcal{M}}}$.
Then, it holds that $B^{\sharp}=(\Omega^{\flat})^{-1}$.
Hamilton’s equations, Eq. (2), are rewritten as

$$ $X_{H}=B^{\sharp}({\mathrm{d}}H).$ (3) $$

The tuple $(H,B,X_{H})$ is called a Poisson system.
On the Darboux coordinates, $B=\sum_{i}{\frac{\partial}{\partial p_{i}}}\wedge{\frac{\partial}{\partial q_{
i}}}$.

### 2.2 Degeneracy, Dissipation, and External Inputs

Degenerate Systems If the state ${\bm{u}}$ is constrained to a submanifold $\tilde{\mathcal{M}}\subset{\mathcal{M}}$, Hamilton’s equations do not directly describe the dynamics.
For example, consider a pair of mass-spring systems, indexed by $i\in\{1,2\}$, with a constraint $q_{1}=q_{2}$ such that the two masses are coupled and always have the same displacement and velocity.
While the Hamiltonian is the sum of those of coupled systems, Hamilton’s equations $\Omega^{\flat}(X_{H})={\mathrm{d}}H$ with the standard symplectic form $\Omega$ do not describe the dynamics.
There are two primary methods for handling such degenerate Hamiltonian systems.

The first approach uses coordinate transformations to reduce the system’s degrees of freedom and define a submanifold $\tilde{\mathcal{M}}\subset{\mathcal{M}}$, where Hamilton’s equations describe the dynamics using the standard symplectic form $\Omega$ or Poisson bivector $B$.
The Darboux-Lie theorem ensures the local existence of such transformations.
The standard Poisson bivector $B$ on $\tilde{\mathcal{M}}$ can also be expressed on ${\mathcal{M}}$ while its bundle map is degenerate, and it represents the coordinate transformation.
Hence, degenerate Hamiltonian systems are included into Poisson systems.
See Appendix A.2 for an example.
Jin et al. () introduced Poisson neural networks (PNNs) that combine neural network-based coordinate transformations with SympNets for degenerate systems (Dinh et al., ; Jin et al., ).
However, this approach lacks interpretability due to non-unique, nonlinear transformations, and extending it to systems with external inputs or dissipation is challenging.

The other approach introduces constraint forces from the coupling, defining constrained Hamiltonian systems.
Finzi et al. () proposed constrained HNNs (CHNNs) for such systems.
However, this method is limited to holonomic constraints (e.g., constraints on configurations).
Furthermore, as the constraints are predetermined, it remains unclear how to learn these constraints directly from data.

Port-Hamiltonian Systems To incorporate dissipation and external inputs, some studies have explored the port-Hamiltonian formulation with neural networks (Zhong et al., ).
Although some have mentioned the underlying Dirac structure (Eidnes et al., ; Neary & Topcu, ; Di Persio et al., ), they rely on the following canonical form without fully leveraging its flexibility:

$$ $\textstyle\begin{pmatrix}\dot{\bm{q}}\\ \dot{\bm{p}}\end{pmatrix}=\left(\begin{pmatrix}O&I\\ -I&O\end{pmatrix}-\begin{pmatrix}O&O\\ O&D({\bm{q}})\end{pmatrix}\right)\nabla H({\bm{q}},{\bm{p}})+\begin{pmatrix}{ \bm{0}}\\ G({\bm{q}})\end{pmatrix}{\bm{f}},$ (4) $$

where $D\in\mathbb{R}^{n\times n}$ represents dissipation, ${\bm{f}}\in\mathbb{R}^{m}$ is the control input vector, and $G\in\mathbb{R}^{n\times m}$ is the control input matrix.
This formulation has several disadvantages.
The matrix $\begin{pmatrix}O&I\\
-I&O\end{pmatrix}$ represents the standard symplectic form $\Omega=\sum_{i}{\mathrm{d}}q_{i}\wedge{\mathrm{d}}p_{i}$ on ${\mathbb{R}}^{2n}$, which cannot handle degeneracies.
Dissipation from multiple sources, such as dampers or friction, is condensed into a single term $D$, limiting interpretability and modeling performance.
Additionally, this formulation cannot handle externally defined velocities, as seen in models of buildings shaken by the ground.

## 3 Poisson-Dirac Neural Networks

To overcome the limitations of existing methods, we propose Poisson-Dirac Neural Networks (PoDiNNs), which leverage the Dirac structure to unify Poisson and port-Hamiltonian systems.
Detailed background theory and proofs of the following theorems are provided in Appendix C.

### 3.1 Dirac Structure

Let $V$ be an $n$-dimensional vector space and $V^{*}$ its dual space, with the natural pairing $\langle\cdot,\cdot\rangle$ between them.
Define the symmetric pairing $\mathopen{\hbox{\set@color${\langle}$}\kern-1.94444pt\leavevmode\hbox{
\set@color${\langle}$}}\cdot,\cdot\mathclose{\hbox{\set@color${\rangle}$}\kern
-1.94444pt\leavevmode\hbox{\set@color${\rangle}$}}$ on $V\oplus V^{*}$ as

$$ $\mathopen{\hbox{\set@color${\langle}$}\kern-1.94444pt\leavevmode\hbox{ \set@color${\langle}$}}({\bm{v}},{\bm{\alpha}}),(\bar{{\bm{v}}},\bar{{\bm{ \alpha}}})\mathclose{\hbox{\set@color${\rangle}$}\kern-1.94444pt\leavevmode \hbox{\set@color${\rangle}$}}=\langle{\bm{\alpha}},\bar{{\bm{v}}}\rangle+ \langle\bar{{\bm{\alpha}}},{\bm{v}}\rangle\quad\text{for}\quad({\bm{v}},{\bm{ \alpha}}),(\bar{{\bm{v}}},\bar{{\bm{\alpha}}})\in V\oplus V^{*}$ $$

where $\oplus$ denotes the direct sum of two vector spaces (Courant, ).

###### Definition 1 ( Courant ( 1990 ); Yoshimura & Marsden ( 2006a ) ) .

A Dirac structure on a vector space $V$ is a vector subspace $D\subset V\oplus V^{*}$ such that $D=D^{\perp}$, where $D^{\perp}$ is the orthogonal complement of $D$ with respect to the pairing $\mathopen{\hbox{\set@color${\langle}$}\kern-1.94444pt\leavevmode\hbox{
\set@color${\langle}$}}\cdot,\cdot\mathclose{\hbox{\set@color${\rangle}$}\kern
-1.94444pt\leavevmode\hbox{\set@color${\rangle}$}}$.

Since $D=D^{\perp}$, $\langle{\bm{\alpha}},\bar{{\bm{v}}}\rangle+\langle\bar{{\bm{\alpha}}},{\bm{v}}
\rangle=0$ for any $({\bm{v}},{\bm{\alpha}}),(\bar{{\bm{v}}},\bar{{\bm{\alpha}}})\in D$, and hence $\langle{\bm{\alpha}},{\bm{v}}\rangle=0$ for any $({\bm{v}},{\bm{\alpha}})\in D$.
A typical Dirac structure can be constructed as follows:

###### Theorem 1 ( Courant ( 1990 ); Yoshimura & Marsden ( 2006a ) ) .

Let $V$ be a vector space and $\Delta$ a vector subspace.
Define the annihilator $\Delta^{\circ}$ of $\Delta$ as
$\Delta^{\circ}=\{{\bm{\alpha}}\in V^{*}\mid\langle{\bm{\alpha}},{\bm{v}}
\rangle=0\text{ for all }{\bm{v}}\in\Delta\}\subset V^{*}.$
Then, $D_{V}=\Delta\oplus\Delta^{\circ}\subset V\oplus V^{*}$ is a Dirac structure on $V$.

A vector bundle ${\mathcal{F}}$ over a manifold ${\mathcal{M}}$ is a collection of vector spaces ${\mathcal{F}}_{\bm{u}}$, called fibers, smoothly assigned to points ${\bm{u}}\in{\mathcal{M}}$, where the manifold ${\mathcal{M}}$ is called the base space.
A typical example is the tangent bundle ${T{\mathcal{M}}}$, where the fiber at ${\bm{u}}$ is the tangent space ${T_{\bm{u}}{\mathcal{M}}}$, and its dual bundle is the cotangent bundle ${T^{*}{\mathcal{M}}}$.
The Whitney sum $\oplus$ defines a vector bundle whose fiber at each point is the direct sum of the fibers of the two bundles at that point.
Given a vector bundle ${\mathcal{F}}$ over ${\mathcal{M}}$ and its dual ${\mathcal{E}}={\mathcal{F}}^{*}$, the Dirac structure is defined as a subbundle of ${\mathcal{F}}\oplus{\mathcal{E}}$.

###### Definition 2 ( Courant ( 1990 ); Yoshimura & Marsden ( 2006a ); van der Schaft & Jeltsema ( 2014 ) ) .

Consider a vector bundle ${\mathcal{F}}$ over a manifold ${\mathcal{M}}$.
A distribution $\Delta$ is a collection of vector subspaces $\Delta_{\bm{u}}\subset{\mathcal{F}}_{\bm{u}}$, each assigned smoothly to ${\mathcal{M}}$ at point ${\bm{u}}$, forming a vector subbundle of ${\mathcal{F}}$.
The annihilator $\Delta^{\circ}$ of $\Delta$ is a collection of the annihilators $\Delta^{\circ}_{\bm{u}}$ of $\Delta_{\bm{u}}$, also forming a subbundle of ${\mathcal{E}}={\mathcal{F}}^{*}$.
Then, a Dirac structure $D$ is constructed as $D=\Delta\oplus\Delta^{\circ}$, which is a subbundle of ${\mathcal{F}}\oplus{\mathcal{E}}$.

If ${\mathcal{F}}\oplus{\mathcal{E}}=TM\oplus{T^{*}{\mathcal{M}}}$, the Dirac structure $D$ can reformulate Hamiltonian, Poisson, and constrained Hamiltonian systems (see Appendix C).
Here, we assume that the fibers ${\mathcal{F}}_{\bm{u}}$ and ${\mathcal{E}}_{\bm{u}}$ of the vector bundles ${\mathcal{F}}$ and ${\mathcal{E}}$ at ${\bm{u}}$ are decomposed as

$$ ${\mathcal{F}}_{\bm{u}}={\mathcal{F}}^{S}_{\bm{u}}\oplus{\mathcal{F}}^{R}_{\bm{ u}}\oplus{\mathcal{F}}^{I}_{\bm{u}}\text{\quad and\quad}{\mathcal{E}}_{\bm{u}} ={\mathcal{E}}^{S}_{\bm{u}}\oplus{\mathcal{E}}^{R}_{\bm{u}}\oplus{\mathcal{E}} ^{I}_{\bm{u}},$ (5) $$

where ${\mathcal{F}}^{S}_{\bm{u}}={T_{\bm{u}}{\mathcal{M}}}$ and ${\mathcal{E}}^{S}_{\bm{u}}={T^{*}_{\bm{u}}{\mathcal{M}}}$.
A point on the fibers ${\mathcal{F}}_{\bm{u}}$ and ${\mathcal{E}}_{\bm{u}}$ is denoted by ${\bm{f}}=({\bm{f}}^{S},{\bm{f}}^{R},{\bm{f}}^{I})$ and ${\bm{e}}=({\bm{e}}^{S},{\bm{e}}^{R},{\bm{e}}^{I})$, respectively.
We refer to ${\bm{f}}\in{\mathcal{F}}_{\bm{u}}$ as *flows*, ${\bm{e}}\in{\mathcal{E}}_{\bm{u}}$ as *efforts*, and both collectively as *port variables*.
Note that our definitions mainly followed those in Duindam et al. (); van der Schaft & Jeltsema (), while we can find other definitions in Yoshimura & Marsden ().

###### Theorem 2 .

Consider vector bundles ${\mathcal{F}}$ and ${\mathcal{E}}={\mathcal{F}}^{*}$ over ${\mathcal{M}}$.
The collection of

$$ $D_{\bm{u}}=\{({\bm{f}},{\bm{e}})\in{\mathcal{F}}_{\bm{u}}\times{\mathcal{E}}_{ \bm{u}}\mid{\bm{f}}=B^{\sharp}_{\bm{u}}({\bm{e}})\}$ $$

for the bundle map $B^{\sharp}_{\bm{u}}:{\mathcal{E}}_{\bm{u}}\to{\mathcal{F}}_{\bm{u}}$ of a bivector $B$ is a Dirac structure $D\subset{\mathcal{F}}\oplus{\mathcal{E}}$.

Here, we define PoDiNNs as a special case of the Poisson-Dirac formulation (Courant, ).

###### Definition 3 (Poisson-Dirac Neural Network) .

Let ${\mathcal{F}}\oplus{\mathcal{G}}$ be a vector bundle over a manifold ${\mathcal{M}}$ defined in Eq. (5), which assigns to each ${\bm{u}}\in{\mathcal{M}}$ a vector space $({\mathcal{F}}^{S}_{\bm{u}}\oplus{\mathcal{F}}^{R}_{\bm{u}}\oplus{\mathcal{F}}
^{I}_{\bm{u}})\times({\mathcal{E}}^{S}_{\bm{u}}\oplus{\mathcal{E}}^{R}_{\bm{u}
}\oplus{\mathcal{E}}^{I}_{\bm{u}})$.
Let $D\subset{\mathcal{F}}\oplus{\mathcal{E}}$ be a Dirac structure defined in Theorem 2.
Let $H:{\mathcal{M}}\to\mathbb{R}$ be an energy function, which determines the effort ${\bm{e}}^{S}={\mathrm{d}}H\in{\mathcal{E}}^{S}_{\bm{u}}$.
The effort ${\bm{e}}^{R}\in{\mathcal{E}}^{R}_{\bm{u}}$ is determined by a mapping $R_{\bm{u}}:{\mathcal{F}}^{R}_{\bm{u}}\to{\mathcal{E}}^{R}_{\bm{u}}$ of the flow ${\bm{f}}^{R}\in{\mathcal{F}}^{R}_{\bm{u}}$.
The effort ${\bm{e}}^{I}(t)\in{\mathcal{E}}^{I}_{\bm{u}}$ is a time-dependent function.
The functions $H$ and $R_{\bm{u}}$ are implemented using neural networks.
If for each ${\bm{u}}(t)$ and $t\in[a,b]$, it holds that

$$ $(({\bm{f}}^{S}(t),{\bm{f}}^{R}(t),{\bm{f}}^{I}(t)),({\bm{e}}^{S}(t),{\bm{e}}^{ R}(t),{\bm{e}}^{I}(t)))\in D_{{\bm{u}}(t)},$ $$

the tuple $(H,B,R,{\bm{e}}^{I},{\bm{f}}^{S})$ is called Poisson-Dirac Neural Networks (PoDiNNs).

**Table 2: Categorization of Components in Different Domains**
| Domain | Mechanical | Rotational | Electro-Magnetic | Hydraulic |  |  |  |
| --- | --- | --- | --- | --- | --- | --- | --- |
| Subdomain | Potential | Kinetic | Potential | Kinetic | Electric | Magnetic | Potential |
| flow (input) | velocity | force | angular velocity | torque | current | voltage | volume flow rate |
| effort (output) | force | velocity | torque | angular velocity | voltage | current | pressure |
| state | displacement | momentum | angle | angular momentum | electric charge | magnetic flux | volume |
| energy-storing | spring | mass | (potential) | inertia | capacitor | inductor | hydraulic tank |
| energy-dissipating | damper | – | friction | – | resistor | resistor | — |
| external input | external force | moving boundary | external torque | – | voltage source | current source | incoming fluid flow |

###### Definition 1 ( Courant ( 1990 ); Yoshimura & Marsden ( 2006a ) ) .

###### Theorem 1 ( Courant ( 1990 ); Yoshimura & Marsden ( 2006a ) ) .

###### Definition 2 ( Courant ( 1990 ); Yoshimura & Marsden ( 2006a ); van der Schaft & Jeltsema ( 2014 ) ) .

###### Theorem 2 .

###### Definition 3 (Poisson-Dirac Neural Network) .

### 3.2 Flows and Efforts for Components

The point ${\bm{u}}$ at the base space ${\mathcal{M}}$ represents the states of the dynamics, which include the displacement of a spring, the momentum of a mass, the angle and angular momentum of a rotating rod, the electric charge of a capacitor.
Intuitively, flows ${\bm{f}}$ are the inputs to components, while efforts ${\bm{e}}$ are the outputs from the components.
Components considered in our formulation are summarized in Table 2, with concrete examples in Appendices B and D.

Similar to Hamiltonian and Poisson systems, the flow ${\bm{f}}^{S}\in{T_{\bm{u}}{\mathcal{M}}}$ represents the vector field $X_{H}$ on ${\mathcal{M}}$, defining the time evolution $\dot{{\bm{u}}}$ of the state ${\bm{u}}$.
The effort ${\bm{e}}^{S}\in{T^{*}_{\bm{u}}{\mathcal{M}}}$ corresponds to the differential ${\mathrm{d}}H$ of the Hamiltonian $H:{\mathcal{M}}\to{\mathbb{R}}$.
Superscript ^S indicates that these components store energy.
In electric circuits, capacitors and inductors are examples, with states as electric charge and magnetic flux, and efforts as the voltages across and currents through them, respectively.
In PoDiNNs, these components are modeled using neural networks that replace the energy functions $H$, similar to HNNs and Dissipative SymODENs.

The flow ${\bm{f}}^{R}\in{\mathcal{F}}^{R}_{\bm{u}}$ and effort ${\bm{e}}^{R}\in{\mathcal{E}}^{R}_{\bm{u}}$ represent energy-dissipating components, such as dampers and resistors.
A damper’s flow $f^{R}$ is the velocity (i.e., the rate of extension or compression), and its effort $e^{R}$ is the force, which are related as $e^{R}=-df^{R}$ for a linear damper with a damping coefficient $d$.
Unlike a spring, a damper does not store its own energy, so its flow $f^{R}$ is not on ${\mathcal{F}}^{S}_{\bm{u}}$ but on ${\mathcal{F}}^{R}_{\bm{u}}$.
Superscript ^R indicates that these components are “resistive.”
However, components that supply energy, such as diodes, can also be classified into this category.
In any case, each component is implemented in PoDiNNs by a neural network that approximates its characteristic $R_{\bm{u}}:{\bm{f}}^{R}\mapsto{\bm{e}}^{R}$.

The flow ${\bm{f}}^{I}\in{\mathcal{F}}^{I}_{\bm{u}}$ and effort ${\bm{e}}^{I}\in{\mathcal{E}}^{I}_{\bm{u}}$ represent external inputs, such as an external force (where force is the effort) or a moving boundary (where velocity is the effort).
Their efforts ${\bm{e}}^{I}$ depend only on time $t$, not on other components.
Their flows ${\bm{f}}^{I}$ are not required for determining the system’s dynamics but represent the outcomes of external inputs, such as the reaction force exerted on the moving boundary.
In PoDiNNs, these external inputs are fed into the neural networks.

### 3.3 Bivector for Representing Coupled Systems

For the coordinates $q_{i}$ and $p_{i}$ on ${\mathcal{M}}$, the basis vectors of the tangent space ${T_{\bm{u}}{\mathcal{M}}}$ are ${\frac{\partial}{\partial q_{i}}}$ and ${\frac{\partial}{\partial p_{i}}}$, respectively.
The $i$-th basis vectors of ${\mathcal{F}}^{R}_{\bm{u}}$, ${\mathcal{F}}^{I}_{\bm{u}}$, ${\mathcal{E}}^{R}_{\bm{u}}$, and ${\mathcal{E}}^{I}_{\bm{u}}$ are denoted by $\xi_{i}^{R}$, $\xi_{i}^{I}$, $\xi^{R*}_{i}$, and $\xi^{I*}_{i}$, respectively.
The bivector $B$ assigns to each point ${\bm{u}}\in{\mathcal{M}}$ wedge products of the basis vectors of the flow space ${\mathcal{F}}_{\bm{u}}$, such as ${\frac{\partial}{\partial p_{i}}}\wedge{\frac{\partial}{\partial q_{j}}}$, ${\frac{\partial}{\partial p_{i}}}\wedge\xi_{j}^{R}$, and ${\frac{\partial}{\partial p_{i}}}\wedge\xi_{j}^{I}$, thereby defining the coupling patterns among components.
For example, ${\frac{\partial}{\partial p_{i}}}\wedge\xi_{j}^{I}$ couples the $j$-th external input with the $i$-th mass $m_{i}$ with the state $p_{i}$.
The effort of the external flow is expressed as $e^{I}_{j}\xi_{j}^{I*}$ when the basis is explicit.
This is fed to the bivector ${\frac{\partial}{\partial p_{i}}}\wedge\xi_{j}^{I}=-\xi_{j}^{I}\wedge{\frac{
\partial}{\partial p_{i}}}$, resulting in $-e^{I}_{j}{\frac{\partial}{\partial p_{i}}}$, which forms part of the flow $f^{S}_{i}{\frac{\partial}{\partial p_{i}}}$ for mass $m_{i}$.
Therefore, we can make the following remarks.

###### Remark 1 (Coupling as Non-zero Elements of Bivector) .

Coupling between two components is represented by a non-zero bivector element, which links the effort of one to the flow of the other.
Thus, by learning the bivector $B$ from the observations of the target system, we can identify the coupling patterns between the system’s components.

###### Remark 2 (Degeneracy of Dynamics as Degeneracy of Bundle Map) .

Constraints between components that cause degenerate dynamics are reflected in degeneracy of the bundle map $B^{\sharp}$.
Thus, by learning the bivector $B$ from the observations of the target system and examining how it degenerates, we can identify the constraints imposed on the system.

###### Remark 3 (Coordinate Transformation by Bivector) .

The bivector $B$ defines the coordinate system, which allows PoDiNNs to learn system dynamics regardless of the coordinate system used for the observations.

###### Remark 4 (Multiphysics) .

Our formulation can represent systems across various domains, as summarized in Table 2.

These remarks are fundamental in system identification and can aid in reverse engineering, as the coupling patterns of circuit elements serve as representations of the circuit diagrams.
Also, PoDiNNs are the first neural-network method to handle multiple domains of dynamical systems and their interactions by leveraging the Dirac structure.
See Appendix B for concrete examples.

###### Remark 1 (Coupling as Non-zero Elements of Bivector) .

###### Remark 2 (Degeneracy of Dynamics as Degeneracy of Bundle Map) .

###### Remark 3 (Coordinate Transformation by Bivector) .

###### Remark 4 (Multiphysics) .

### 3.4 Discussions for Comparisons and Limitations

As discussed above, PoDiNNs are the first model to cover degenerate dynamics, dissipation, external inputs, and coordinate transformations in a unified manner.

Previous models, including HNNs, LNNs, NSFs, PNNs, and Dissipative SymODENs, approximate the Hamiltonian $H$, Lagrangian $L$, or dissipative term $D$ using neural networks, which implicitly learn the relationships between variables.
However, these relationships are difficult to extract due to the implicit nature of the learning and the high nonlinearity of the networks.
In contrast, PoDiNNs separate the coupling patterns $B$ from the energy functions $H$, improving interpretability and generalization performance.

In an absolute coordinate system, the displacements of springs were computed from the absolute positions of springs’ edges internally within the energy function.
Although PoDiNNs can handle such system, the bivector for energy-storing components takes the standard form $B=\sum_{i}{\frac{\partial}{\partial p_{i}}}\wedge{\frac{\partial}{\partial q_{
i}}}$, which hinders the identification of their coupling patterns.
PoDiNNs are still able to identify coupling patterns involving energy-dissipating components and external inputs.
To identify the coupling patterns between energy-storing components, it is necessary to employ a relative coordinate system based on displacements, rather than absolute positions.

LNNs and NSFs are designed for non-canonical Hamiltonian systems, where the Hamiltonian vector field $X_{H}$ is implicitly defined as $\Omega^{\flat}(X_{H})={\mathrm{d}}H$ (Cranmer et al., ; Chen et al., ).
This approach requires matrix inversion to compute $X_{H}$, which is computationally expensive, numerically unstable, and not directly applicable to degenerate dynamics.
In contrast, PoDiNNs explicitly define $X_{H}$ as part of ${\bm{f}}=B^{\sharp}({\bm{e}})$, offering faster and more stable computations that also handle degenerate dynamics.

In mechanical systems, energy-dissipating components such as dampers and friction are characterized by first setting the velocities set, from which the corresponding forces are then derived.
As a result, flow is always defined as velocity, and effort as force.
In the electric circuits, however, the flow for resistors and diodes can be either current or voltage, depending on their coupling with other components.
In practice, it is advisable to include an abundance of both types of components.
Any excess components will either be ignored or exhibit redundant characteristics, as demonstrated in the experiments.

The dynamics of electric circuits are generally described by differential-algebraic equations (DAEs).
For example, when a capacitor is connected to a direct voltage source in parallel, infinite current instantaneously flows into the capacitor, and its voltage matches that of the direct voltage source.
This behavior cannot be represented by ODEs alone.
While PoDiNNs cannot fully represent such systems, they can still describe a wide range of systems and expand the scope of modeling unknown dynamical systems from observations.

## 4 Experiments and Results

### 4.1 Experimental Settings

Figure: Figure 2: Diagrams of systems that provide datasets. Detailed definitions are found in Appendix D.
Refer to caption: x2.png

Datasets We evaluated PoDiNNs and related methods to demonstrate their modeling performance using seven simulation datasets, as shown in Fig. 2.
Due to page limitations, we briefly overview their characteristics.
The full explanations can be found in Appendix D.

For the mechanical domain, we prepared three mass-spring(-damper) systems (a)–(c), with the velocities of the masses as observations.
In the absolute coordinate system, we used the absolute positions of the springs’ ends as the springs’ states, and in the relative coordinate system, their displacements.
As external inputs, an external force $F$ is applied to system (a), and a moving boundary $b$ is coupled with system (b).
Dissipative SymODENs cannot directly account for the latter.
System (c) has redundant observations; while springs are coupled only with two masses, the displacements of all five springs were provided as observations.
CHNNs cannot model this redundancy, as it does not arise from holonomic constraints, but PNNs can.
In all three systems, the springs and dampers exhibit nonlinear characteristics.

We selected two nonlinear electric circuits, (d) the FitzHugh-Nagumo model and (e) Chua’s circuit, from the electro-magnetic domain (Izhikevich & FitzHugh, ; Chua, ).
The capacitor voltage and inductor current were used as the observations.

We also used two multiphysics systems (f) and (g).
In system (f), a DC motor bridges an electric circuit in the electro-magnetic domain and a pendulum in the rotational domain.
In system (g), a hydraulic tank in the hydraulic domain is connected to two cylinders with pistons in the mechanical domain, each of which is also connected to a fixed wall via a spring and damper.
An external force is applied to the smaller piston, which moves the larger piston through the fluid in the tank.

Implementation Details We implemented all experimental code from scratch using Python v3.11.9, along with numpy v1.26.4, scipy v1.12.1, pytorch v2.3.1, and desolver v4.1.1 (Paszke et al., ).
See also Appendix D for more details.

We compared Neural ODEs, Dissipative SymODENs, PNNs, and PoDiNNs.
Neural ODEs were evaluated on all datasets.
Dissipative SymODENs were tested on systems (a) and (b), while PNNs were evaluated on system (c), in the absolute coordinate system.
Other combinations were out of scope of the original studies.
For PNNs, we used HNNs in place of SympNets for a fair comparison.
For Dissipative SymODENs and PoDiNNs, we assumed that the kinetic energy in the mechanical domain could be expressed as $\frac{1}{2m}p^{2}$, where $m$ is the mass and $p$ is the momentum.
Therefore, we employed this form with a learnable parameter $m$, rather than a neural network.
The same approach was applied to capacitors and inductors for PoDiNNs.
All other components were assumed to be nonlinear and were modeled using neural networks.
In the absolute coordinate system, the potential energy was modeled for all configurations $q_{1},\dots,q_{n}$ plus the position $q_{b}$ of the moving boundary collectively by a single neural network.
In the relative coordinate system or other domains, the potential energy was modeled separately for each component.
We assumed that the number of energy-dissipating components and the nature of their flows (e.g., current or voltage in electric circuits) are known, and also examined the impact of inaccurate assumptions.

Each model was trained using one-step predictions on the training subset.
Specifically, given a random state snapshot ${\bm{u}}^{(n)}$ at $n$-th step, each model predicted the next state $\tilde{{\bm{u}}}^{(n+1)}$ after a time step $\Delta t$.
Then, all parameters were updated to minimize the squared error between the predicted state $\tilde{{\bm{u}}}^{(n+1)}$ and the ground truth ${\bm{u}}^{(n+1)}$, normalized by state standard deviations.
It is known that longer prediction steps can improve robustness against noise (Chen et al., ).
However, as we aimed to purely compare the representational performance of models, no noise was added to the datasets.
We confirmed that longer prediction steps only led to performance degradation.

Evaluation Metrics We evaluated models using the accuracy of the solution to the initial value problem on the test subset.
Starting from the initial value of each trajectory, each model predicted the entire trajectory of $N$ steps and calculated the mean squared error (MSE) between the predicted state $\tilde{\bm{u}}^{(n)}$ and the ground truth ${\bm{u}}^{(n)}$ at each step indexed by $n$.
The mean of these MSEs across all trajectories and all time steps was computed as the evaluation metric, referred to as the overall MSE;

$$ $\textstyle MSE(\tilde{\bm{u}};{\bm{u}})=\frac{1}{N}\sum_{n=1}^{N}[\mathrm{MSE} (\tilde{\bm{u}}^{(n)},{\bm{u}}^{(n)})].$ (6) $$

Lower values indicate better performance.
Additionally, we defined the valid prediction time (VPT) as the ratio of the number of steps taken before the MSE first exceeds a certain threshold $\theta$ to the total length $N$ of the test trajectory (Botev et al., ; Jin et al., ; Vlachas et al., );

$$ $\textstyle VPT(\tilde{\bm{u}};{\bm{u}})=\frac{1}{N}\arg\max_{n_{f}}\{n_{f}| \mathrm{MSE}(\tilde{\bm{u}}^{(n)},{\bm{u}}^{(n)})<\theta\mbox{ for all }n\leq n _{f}\}.$ (7) $$

Higher values indicate better performance.
The threshold $\theta$ was set to ensure that most models achieved VPTs between 0.1 and 0.9.
For each dataset, models were trained and evaluated from scratch over 10 trials per dataset.

**Table 3: Experimental Results.**
|  | Mass-Spring-Damper Systems |  |  |  |  |  |  |  |  |  |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
| Dataset | with external force | with moving boundary | (c) with redundancy |  |  |  |  |  |  |  |
| Coordinate | (a) relative | (a’) absolute | (b) relative | (b’) absolute | relative |  |  |  |  |  |
| Model | MSE$\downarrow$ | VPT$\uparrow$ | MSE$\downarrow$ | VPT$\uparrow$ | MSE$\downarrow$ | VPT$\uparrow$ | MSE$\downarrow$ | VPT$\uparrow$ | MSE$\downarrow$ | VPT$\uparrow$ |
| Neural ODEs | 4.90<br>​$\pm$<br>0.27 | 0.128<br>​$\pm$<br>0.021 | 7.68<br>​$\pm$<br>1.07 | 0.097<br>​$\pm$<br>0.008 | 7.43<br>​$\pm$<br>1.19 | 0.153<br>​$\pm$<br>0.039 | 5.02<br>​$\pm$<br>0.56 | 0.135<br>​$\pm$<br>0.052 | 2490.61<br>​$\pm$<br>1847.24 | 0.099<br>​$\pm$<br>0.004 |
| HNN Variants^∗ | — | 8.31<br>​$\pm$<br>0.56 | 0.104<br>​$\pm$<br>0.017 | — | 5.92<br>​$\pm$<br>0.12 | 0.001<br>​$\pm$<br>0.000 | 634.22<br>​$\pm$<br>300.01 | 0.000<br>​$\pm$<br>0.000 |  |  |
| PoDiNNs | 4.33<br>​$\pm$<br>0.26 | 0.622<br>​$\pm$<br>0.002 | 7.02<br>​$\pm$<br>0.49 | 0.437<br>​$\pm$<br>0.053 | 0.26<br>​$\pm$<br>1.12 | 0.856<br>​$\pm$<br>0.015 | 3.74<br>​$\pm$<br>0.84 | 0.581<br>​$\pm$<br>0.040 | 0000.11<br>​$\pm$<br>0.02000 | 0.863<br>​$\pm$<br>0.017 |
|  | $\times 10^{-1}$ | $\theta=10^{-3}$ | $\times 10^{-1}$ | $\theta=10^{-3}$ | $\times 10^{-2}$ | $\theta=10^{-4}$ | $\times 10^{-2}$ | $\theta=10^{-4}$ | $\times 10^{-1}$ | $\theta=10^{-3}$ |

Figure: Figure 3: Visualizations of example results. Each top panel shows ground truth trajectories, while the other panels show the absolute errors of all 10 trials in semi-transparent color. See also Fig. A1.
Refer to caption: x9.png

### 4.2 Results

Numerical Performance The results, summarized in Table 3, show that PoDiNNs consistently outperformed all other models across all datasets and metrics.
VPTs show that PoDiNNs provided stable predictions for durations 1.5 to 20 times longer.
These differences were statistically significant, with p-values less than 0.0005, as evaluated by the Mann-Whitney U test.
In some cases, the difference in MSE is small, while that in VPT is substantial.
This happens when a model ignores oscillations in the trajectory and outputs an average path, which keeps the MSE low but results in poor VPT.
This suggests that VPT is a more reliable metric than MSE (Botev et al., ; Jin et al., ; Vlachas et al., ).
We also show example trajectories and absolute errors in Figs. 3 and A1.

PoDiNNs perform well both in the absolute and relative coordinate systems, showing their adaptability to different coordinate systems.
Even in the absolute coordinate system, PoDiNNs decompose dissipations and external inputs into coupling patterns and individual characteristics, providing a more effective inductive bias, whereas Dissipative SymODENs treat dissipations and external inputs using black-box functions $D$ and $G$.

PoDiNNs demonstrate excellent accuracy in system (c) because they successfully simplify the dynamics by correctly identifying the degeneracy from high-dimensional observations.
In contrast, PNNs struggled to learn the appropriate coordinate transformation.
Although the Darboux-Lie theorem guarantees the existence of such a transformation, it does not imply that learning it is straightforward.
Also, Neural ODEs lack guarantees for energy conservation, resulting in diverging trajectories.

Identifying Coupling Patterns and Component Characteristics (e) We examined how the bivector $B$ identifies coupling patterns in system (e), Chua’s circuit.
When two coefficients differed by a factor of 1000 or more, we considered the larger one as a detected coupling and the smaller one as effectively zero, indicating no coupling.
In all 10 trials, we obtained the bivector $B=-a{\frac{\partial}{\partial\psi}}\wedge{\frac{\partial}{\partial Q_{2}}}+b
\xi^{R}_{1}\wedge{\frac{\partial}{\partial Q_{2}}}-c_{1}\xi^{R}_{1}\wedge{
\frac{\partial}{\partial Q_{1}}}+c_{2}\xi^{R}_{2}\wedge{\frac{\partial}{
\partial Q_{1}}}$, with trial-wise positive parameters $a$, $b$, $c_{1}$, and $c_{2}$.
We emphasize that the coefficients of other bivector elements, such as ${\frac{\partial}{\partial\psi}}\wedge{\frac{\partial}{\partial Q_{1}}}$ and $\xi^{R}_{2}\wedge{\frac{\partial}{\partial Q_{2}}}$, were effectively zero.
Because $R_{1}$ and $R_{2}$ cannot be directly observed and their indices are interchangeable, they were appropriately reordered for analysis.

We normalized the coefficients of the bivector elements and the characteristics of system components, as their scales cancel each other out.
Then, from the learned bivector $B$, we can derive Kirchhoff’s current laws, $I_{C_{1}}=-I_{R_{1}}+I_{R_{2}}$ and $I_{C_{2}}=-I_{L}+I_{R_{1}}$, and Kirchhoff’s voltage laws, $V_{L}=V_{C_{2}}$, $V_{R_{1}}=V_{C_{1}}-V_{C_{2}}$, and $V_{R_{2}}=-V_{C_{1}}$.
This allows us to construct a circuit diagram, which perfectly matches that of Chua’s circuit, shown in Fig. 2 (e).

Figure: Figure 4: The identified characteristics of $R_{2}$.
Refer to caption: x13.png

While coupling patterns in more complex circuits may not be unique due to Norton’s and Thevenin’s theorems, our formulation can learn one valid realization.
Figure 4 illustrates the input-output relationship of the neural network representing $R_{2}$.
The results from all 10 trials, shown as colored dashed lines, almost perfectly overlap the true relationship, shown as a black dashed line, within the range of $\pm 2.5$.
This demonstrates that PoDiNNs accurately identified the component characteristics, even though its response was never observed directly.
The region beyond this range was not included in the training subset, so it is expected that PoDiNNs did not learn the relationship there.

Identifying Coupling Patterns of System (g) We also examined the bivector $B$ identified in system (g), hydraulic tank.
In all 10 trials, we consistently obtained the bivector $B={\frac{\partial}{\partial p_{1}}}\wedge{\frac{\partial}{\partial V}}-0.3({
\frac{\partial}{\partial p_{2}}}\wedge{\frac{\partial}{\partial V}})+{\frac{
\partial}{\partial p_{1}}}\wedge{\frac{\partial}{\partial q_{1}}}+{\frac{
\partial}{\partial p_{2}}}\wedge{\frac{\partial}{\partial q_{2}}}$ between the energy-storing components, with the coefficients $1.0$ or $-0.3$, accurate to five significant figures.
The coefficients of the first two terms, $1.0$ and $-0.3$, correspond to the cross-sectional areas $a_{1}$ and $a_{2}$ of two cylinders attached at the bottom of the tank, with the negative sign indicating that the flow directions are opposite.
The remaining two terms represent the couplings between the masses and springs.
All other couplings—between the dampers and masses, ${\frac{\partial}{\partial p_{1}}}\wedge{\frac{\partial}{\partial d_{1}}}$ and ${\frac{\partial}{\partial p_{2}}}\wedge{\frac{\partial}{\partial d_{2}}}$, and between the external force and the mass, ${\frac{\partial}{\partial p_{2}}}\wedge{\frac{\partial}{\partial\xi^{S}}}$—were also identified with non-zero coefficients, even though the overall scale is indeterminate due to the cancellation of gravitational acceleration, fluid density, and masses.

Figure: Table 4: Impact of # Components and VPT.

Impact of Number of Hidden Components (b) The states of energy-storing components are provided as observations, and external inputs are typically known (or inferred using methods like neural CDE (Kidger et al., )).
PoDiNNs also require specifying the number and type of energy-dissipating components, which are usually unknown.
To assess the impact of the assumed number of components, $n_{d}$, we tested system (b) in the relative coordinate system (see Table 4).
The correct number of dampers is $n_{d}=3$.
When $n_{d}<3$, performance was extremely poor, indicating incorrect dynamics due to missing dampers.
When $n_{d}>3$, performance was similar to the case for $n_{d}=3$.
Redundant dampers provided the extra parameters, which sometimes made optimization easier, but were often ignored by learning identical properties to existing dampers, by adopting a zero damping coefficient, or by having zero coupling strength.
Interestingly, this trend appeared in both the training and test subsets, which suggests that assessing performance on the training subset can help identify the correct number of dampers.
In this way, PoDiNNs offer interpretable insights into the system’s internal structure.

See Appendix E for additional visualizations and analyses.

## 5 Conclusion

In this study, we proposed *Poisson-Dirac Neural Networks* (PoDiNNs), which use a Dirac structure to unify the port-Hamiltonian and Poisson formulations.
PoDiNNs offer a unified framework capable of handling various domains of dynamical systems, identifying internal coupling patterns, learning component-wise characteristics, and effectively modeling multiphysics systems.
Our experiments with mechanical, rotational, electro-magnetic, and hydraulic systems validate these capabilities.
Developing methods to address dynamical systems that are described by DAEs and partial differential equations remains for future research.

## 6 Ethics Statement

This study is purely focused on dynamical systems and modeling, and it is not expected to have any direct negative impact on society or individuals.

## 7 Reproducibility Statement

The environment, datasets, methods, evaluation metrics, and other experimental settings are given in Section 4.1 and Appendix D.
For full reproducibility, it is recommended to run the source code attached as supplementary material.

## Appendix A Examples of Formulations

Consider the time evolution of a point ${\bm{u}}$ on a manifold ${\mathcal{M}}$.
If a local coordinate on ${\mathcal{M}}$ is denoted by $x_{i}$, the corresponding basis vector of the tangent space ${T_{\bm{u}}{\mathcal{M}}}$ is denoted by ${\frac{\partial}{\partial x_{i}}}$.
A vector field $X$ on ${\mathcal{M}}$ is expressed as $X=\sum_{i}X_{i}{\frac{\partial}{\partial x_{i}}}$.
If the curve ${\bm{u}}(t)$ satisfies $\dot{\bm{u}}(t)=X_{{\bm{u}}(t)}$, $X_{i}$ represents the local rate of change of the state ${\bm{u}}(t)$ in the $x_{i}$-direction, i.e., $\dot{u}_{i}=X_{i}$.
For a function $F:{\mathcal{M}}\to{\mathbb{R}}$, its differential ${\mathrm{d}}F$ is given by ${\mathrm{d}}F=\sum_{i}{\frac{\partial F}{\partial x_{i}}}{\mathrm{d}}x_{i}$.
The differential is a covector field, that is, a collection of points on cotangent spaces ${T^{*}_{\bm{u}}{\mathcal{M}}}$.
The basis vector of the cotangent space ${T^{*}_{\bm{u}}{\mathcal{M}}}$ is ${\mathrm{d}}x_{i}$, and the following relation holds: ${\mathrm{d}}x_{i}{\frac{\partial}{\partial x_{j}}}=1$ if $i=j$ and 0 otherwise.

When a function $F$ defines dynamics, a two-tensor field links a vector field $X$ and the covector field ${\mathrm{d}}F$.
Such a two-tensor field can be a symplectic form, Poisson bivector, or Riemannian metric.
For example, a link with the negative of a Riemannian metric defines a gradient flow.
A symplectic form and Poisson bivector are defined using the wedge product $\wedge$, which satisfies the skew-symmetry ${\mathrm{d}}x_{i}\wedge{\mathrm{d}}x_{j}=-{\mathrm{d}}x_{j}\wedge{\mathrm{d}}x
_{i}$, and the relation $({\mathrm{d}}x_{i}\wedge{\mathrm{d}}x_{j})({\frac{\partial}{\partial x_{i}}})=
{\mathrm{d}}x_{j}$.
Marsden & Ratiu () and Hairer et al. () have thoroughly discussed how to describe dynamical systems using these geometric objects.
While theoretical details are left to these textbooks, this section introduces a mass-spring system as a concrete example.

### A.1 Mass-Spring System

Consider a mass-spring system with spring constant $k$ and mass $m$.
We will write its dynamics by several formulations.

Canonical Hamiltonian Systems In the Darboux coordinates (i.e., on the cotangent bundle ${T^{*}{\mathcal{Q}}}$), the generalized coordinate $q$ is the displacement of spring $k$, and the generalized momentum $p$ is obtained as $p=mv$ for the velocity $v$ of mass $m$.
The manifold ${\mathcal{M}}$ is a 2-dimensional Euclidean space $\mathbb{R}^{2}$.
The Hamiltonian $H$ is $H(q,p)=\frac{1}{2m}p^{2}+\frac{1}{2}kq^{2}$, and the symplectic form $\Omega$ is standard, i.e., $\Omega={\mathrm{d}}q\wedge{\mathrm{d}}p$.

Hamilton’s equations state that $(\Omega^{\flat})(X_{H})=({\mathrm{d}}q\wedge{\mathrm{d}}p)(\dot{q}{\frac{
\partial}{\partial q}}+\dot{p}{\frac{\partial}{\partial p}})=\dot{q}{\mathrm{d
}}p-\dot{p}{\mathrm{d}}q$ equals ${\mathrm{d}}H={\frac{\partial H}{\partial q}}{\mathrm{d}}{q}+{\frac{\partial H
}{\partial p}}{\mathrm{d}}{p}$, leading to the equations of motion, $\dot{q}={\frac{\partial H}{\partial p}}=p/m,\ \dot{p}=-{\frac{\partial H}{
\partial q}}=-kq$.

Non-Canonical Hamiltonian Systems On the tangent bundle ${\mathcal{M}}={T{\mathcal{Q}}}$, the velocity $v$ is used in place of the momentum $p$, and the symplectic form $\Omega$ is the Lagrangian 2-form $\Omega=m({\mathrm{d}}q\wedge{\mathrm{d}}v)$.
The Hamiltonian $H$ is $H(q,v)=\frac{1}{2}mv^{2}+\frac{1}{2}kq^{2}$.

Hamilton’s equations state that $(\Omega^{\flat})(X_{H})=m({\mathrm{d}}q\wedge{\mathrm{d}}v)(\dot{q}{\frac{
\partial}{\partial q}}+\dot{v}{\frac{\partial}{\partial v}})=m\dot{q}{\mathrm{
d}}v-m\dot{v}{\mathrm{d}}q$ equals ${\mathrm{d}}H={\frac{\partial H}{\partial q}}{\mathrm{d}}{q}+{\frac{\partial H
}{\partial v}}{\mathrm{d}}{v}$, leading to the equations of motion, $\dot{q}=\frac{1}{m}{\frac{\partial H}{\partial v}}=v,\ \dot{v}=-\frac{1}{m}{
\frac{\partial H}{\partial q}}=-kq/m$.

Poisson Systems In the Darboux coordinates, the Poisson bivector $B$ is $B={\frac{\partial}{\partial p}}\wedge{\frac{\partial}{\partial q}}$.
Hamilton’s equations state that $X_{H}=\dot{q}{\frac{\partial}{\partial q}}+\dot{p}{\frac{\partial}{\partial p}}$ equals $B^{\sharp}({\mathrm{d}}H)=({\frac{\partial}{\partial p}}\wedge{\frac{\partial}
{\partial q}})({\frac{\partial H}{\partial q}}{\mathrm{d}}{q}+{\frac{\partial H
}{\partial p}}{\mathrm{d}}{p})={\frac{\partial H}{\partial p}}{\frac{\partial}
{\partial q}}-{\frac{\partial H}{\partial q}}{\frac{\partial}{\partial p}}$.
The equations of motion are $\dot{q}={\frac{\partial H}{\partial p}}=p/m,\ \dot{p}=-{\frac{\partial H}{
\partial q}}=-kq$.

On the tangent bundle ${\mathcal{M}}={T{\mathcal{Q}}}$, the Poisson bivector $B$ is $B=\frac{1}{m}({\frac{\partial}{\partial v}}\wedge{\frac{\partial}{\partial q}})$, and the Hamiltonian is $H(q,v)=\frac{1}{2}mv^{2}+\frac{1}{2}kq^{2}$.
Hamilton’s equations state that $X_{H}=\dot{q}{\frac{\partial}{\partial q}}+\dot{v}{\frac{\partial}{\partial v}}$ equals $B^{\sharp}({\mathrm{d}}H)=\frac{1}{m}({\frac{\partial}{\partial v}}\wedge{
\frac{\partial}{\partial q}})({\frac{\partial H}{\partial q}}{\mathrm{d}}{q}+{
\frac{\partial H}{\partial v}}{\mathrm{d}}{v})=\frac{1}{m}{\frac{\partial H}{
\partial v}}{\frac{\partial}{\partial q}}-\frac{1}{m}{\frac{\partial H}{
\partial q}}{\frac{\partial}{\partial v}}$.
The equations of motion are $\dot{q}=\frac{1}{m}{\frac{\partial H}{\partial v}}=v,\ \dot{v}=-\frac{1}{m}{
\frac{\partial H}{\partial q}}=-kq/m$.

### A.2 Constrained Mass-Spring System

Consider a pair of mass-spring systems, indexed by $i\in\{1,2\}$, and introduce a constraint such that the two masses are coupled and always have the same displacement and velocity; the dynamics is degenerate.
This constraint is expressed as $q_{1}=q_{2}$.
The Hamiltonian is given simply by the sum of two systems, $H({\bm{q}},{\bm{p}})=\frac{1}{2m_{1}}p_{1}^{2}+\frac{1}{2m_{2}}p_{2}^{2}+\frac
{1}{2}k_{1}q_{1}^{2}+\frac{1}{2}k_{2}q_{2}^{2}$.
However, Hamilton’s equations $\Omega^{\flat}(X_{H})={\mathrm{d}}H$ with the standard symplectic form $\Omega$ cannot describe the dynamics.

Degenerate Systems with Coordinate Transformation Define new coordinates $q=q_{1}=q_{2}$ and $p=\frac{m_{1}+m_{2}}{m_{1}}p_{1}=\frac{m_{1}+m_{2}}{m_{2}}p_{2}$, and consider the submanifold $\tilde{\mathcal{M}}$ spanned by $(q,p)$.
The Hamiltonian $H$ is unchanged but rewritten as $H(q,p)=\frac{1}{2(m_{1}+m_{2})}p^{2}+\frac{1}{2}(k_{1}+k_{2})q^{2}$.
With the standard symplectic form $\Omega={\mathrm{d}}q\wedge{\mathrm{d}}p$ on $\tilde{\mathcal{M}}$, the equations of motion are written as $\dot{q}=\frac{p}{m_{1}+m_{2}},\ \dot{p}=-(k_{1}+k_{2})q$.

Degenerate Systems with Constraint Force The constraint on the coordinates, $q_{1}-q_{2}=0$, is naturally satisfied by the constraint on the velocities, $\dot{q}_{1}-\dot{q}_{2}=0$, which is written as the constraint on momenta, $p_{1}/m_{1}-p_{2}/m_{2}=0$.
The constants of motion are $q_{1}-q_{2}$ and $p_{1}/m_{1}-p_{2}/m_{2}$.
A constraint force can be introduced to satisfy these constraints, yielding the same equation as above.
Finzi et al. () proposed CHNNs by combining this approach with neural networks.
However, the automatic derivation of the constraint on momenta from that on the coordinates is non-trivial, and it remains unclear how to learn constraints from data.

Degenerate Systems as Poisson Systems The equations of motion on the submanifold $\tilde{\mathcal{M}}$, $\dot{q}=\frac{p}{m_{1}+m_{2}},\ \dot{p}=-(k_{1}+k_{2})q$, can be rewritten in the original coordinate system as $\dot{q}_{1}=\frac{p_{1}+p_{2}}{m_{1}+m_{2}},\ \dot{q}_{2}=\frac{p_{1}+p_{2}}{m
_{1}+m_{2}},\ \dot{p}_{1}=-\frac{m_{1}}{m_{1}+m_{2}}(k_{1}q_{1}+k_{2}q_{2}),\
\dot{p}_{2}=-\frac{m_{2}}{m_{1}+m_{2}}(k_{1}q_{1}+k_{2}q_{2})$.
Even on the original manifold ${\mathcal{M}}$, these equations are obtained from Hamilton’s equations $X=B^{\sharp}({\mathrm{d}}H)$ with the Poisson bivector $B=\frac{1}{m_{1}+m_{2}}(m_{1}{\frac{\partial}{\partial p_{1}}}+m_{2}{\frac{
\partial}{\partial p_{2}}})\wedge({\frac{\partial}{\partial q_{1}}}+{\frac{
\partial}{\partial q_{2}}})$ on ${\mathcal{M}}$, while its bundle map $B^{\sharp}$ is degenerate.
This fact implies that, by adjusting the Poisson bivector from data, we can learn the Hamiltonian systems with constraints and identify how the dynamics is degenerate.

On the tangent bundle ${\mathcal{M}}={T{\mathcal{Q}}}$, the Hamiltonian is $H(q_{1},q_{2},v_{1},v_{2})=\frac{1}{2}m_{1}v_{1}^{2}+\frac{1}{2}m_{2}v_{2}^{2}
+\frac{1}{2}k_{1}q_{1}^{2}+\frac{1}{2}k_{2}q_{2}^{2}$, and the Poisson bivector $B$ is $B=\frac{1}{m_{1}+m_{2}}({\frac{\partial}{\partial v_{1}}}+{\frac{\partial}{
\partial v_{2}}})\wedge({\frac{\partial}{\partial q_{1}}}+{\frac{\partial}{
\partial q_{2}}})$.
Then, the equations of motion are $\dot{q}_{1}=\frac{m_{1}v_{1}+m_{2}v_{2}}{m_{1}+m_{2}},\ \dot{q}_{2}=\frac{m_{1
}v_{1}+m_{2}v_{2}}{m_{1}+m_{2}},\ \dot{v}_{1}=-\frac{k_{1}q_{1}+k_{2}q_{2}}{m_
{1}+m_{2}},\ \dot{v}_{2}=-\frac{k_{1}q_{1}+k_{2}q_{2}}{m_{1}+m_{2}}$.

## Appendix B Poisson Systems

### B.1 Mechanical Systems

Using the mass-spring systems described above, we provide concrete examples of how the bivector represents coupling and constraints in mechanical systems.

Coupled Systems as Poisson Systems Consider two masses and two springs, indexed by $i\in\{1,2\}$, coupled in sequence from a fixed wall.
Let $q_{i}$ denote the displacement of $i$-th spring.
The equations of motion are given by $\dot{q}_{1}=p_{1}/m_{1}$, $\dot{q}_{2}=p_{2}/m_{2}-p_{1}/m_{1}$, $\dot{p}_{1}=-k_{1}q_{1}+k_{2}q_{2}$, and $\dot{p}_{2}=-k_{2}q_{2}$.
The bivector $B$ leading to these equations is $B={\frac{\partial}{\partial p_{1}}}\wedge{\frac{\partial}{\partial q_{1}}}-{
\frac{\partial}{\partial p_{1}}}\wedge{\frac{\partial}{\partial q_{2}}}+{\frac
{\partial}{\partial p_{2}}}\wedge{\frac{\partial}{\partial q_{2}}}$ for the Hamiltonian $H({\bm{q}},{\bm{p}})=\frac{1}{2m_{1}}p_{1}^{2}+\frac{1}{2m_{2}}p_{2}^{2}+\frac
{1}{2}k_{1}q_{1}^{2}+\frac{1}{2}k_{2}q_{2}^{2}$.
These terms indicate the couplings between $p_{1}$ and $q_{1}$, $p_{1}$ and $q_{2}$, and $p_{2}$ and $q_{2}$.
The negative coefficient indicates that the coupling between $p_{1}$ and $q_{2}$ is in the opposite direction.

Degenerate Systems as Poisson Systems As shown in Appendix A.2, degenerate systems can be expressed as Poisson systems.
Recall the case of a pair of mass-spring systems, indexed by $i\in\{1,2\}$, constrained so that the two masses always have the same displacement and velocity.
This system is expressed with the Poisson bivector $B=\frac{1}{m_{1}+m_{2}}(m_{1}{\frac{\partial}{\partial p_{1}}}+m_{2}{\frac{
\partial}{\partial p_{2}}})\wedge({\frac{\partial}{\partial q_{1}}}+{\frac{
\partial}{\partial q_{2}}})$, which is degenerate in the sense that its bundle map $B^{\sharp}$ is degenerate.
This indicates the absence of a corresponding symplectic form.
Conversely, by learning the bivector $B$ and examining how it degenerates, one can identify the constraints imposed on the system.

The coefficient $\frac{m_{1}}{m_{1}+m_{2}}$ for the term ${\frac{\partial}{\partial p_{1}}}\wedge{\frac{\partial}{\partial q_{1}}}$ indicates the coupling strength between mass $p_{1}$ and spring $q_{1}$.
The effort $e^{S}_{q_{1}}={\frac{\partial H}{\partial q_{1}}}$ from spring $q_{1}$ is distributed to the masses such that $\frac{m_{1}}{m_{1}+m_{2}}e_{q_{1}}$ goes to mass $m_{1}$ and $\frac{m_{2}}{m_{1}+m_{2}}e_{q_{1}}$ goes to mass $m_{2}$.

Coordinate Transformation by Bivector As shown in Appendix A.2, the elements of the Poisson bivector $B$ in a Poisson system depend on whether the dynamics are defined on the tangent bundle ${T{\mathcal{Q}}}$ (using generalized velocities as part of the state) or the cotangent bundle ${T^{*}{\mathcal{Q}}}$ (using generalized momenta as part of the state).
Despite this difference, both coordinate systems are represented by Poisson bivectors.
Thus, by learning the bivector $B$ to approximate the dynamics of a given system, one can implicitly learn the coordinate system employed by that system.

### B.2 Electric Circuits

Our formulation is applicable to electric circuits, as summarized in Table 2.
Let $I_{X}$ and $V_{X}$ denote the current through and voltage across a circuit element $X$, respectively.

A capacitor with capacitance $C$ has the electric charge $Q$ as its state, and stores the energy $H_{C}=\frac{1}{2C}Q^{2}$.
Its flow is the current $I_{C}$, which leads to the change in the electric charge $Q$ as $\dot{Q}=I_{C}$.
Its effort is the voltage $V_{C}$ because the stored electric charge $Q$ generates the voltage $V_{C}$ as $V_{C}=\frac{Q}{C}={\frac{\partial H_{C}}{\partial Q}}$.

An inductor with inductance $L$ has the magnetic flux $\varphi$ as its state, and stores the energy $H_{L}=\frac{1}{2L}\varphi^{2}$.
Its flow is the voltage $V_{L}$, which leads to the change in the magnetic flux $\varphi$ as $\dot{\varphi}=V_{L}$.
Its effort is the current $I_{L}$ because the magnetic flux $\varphi$ generates the current $I_{L}$ as $I_{L}=\frac{\varphi}{L}={\frac{\partial H_{L}}{\partial\varphi}}$.

Consider a system composed of an inductor $L$ and a capacitor $C$ coupled in series.
The state space ${\mathcal{M}}$ is the space of the magnetic flux $\varphi$ and electric charge $Q$.
The total energy is $H=H_{L}+H_{C}$, and its differential is ${\mathrm{d}}H=\frac{\varphi}{L}{\mathrm{d}}\varphi+\frac{Q}{C}{\mathrm{d}}Q$.
Define a bivector $B={\frac{\partial}{\partial\varphi}}\wedge{\frac{\partial}{\partial Q}}$, which leads to

$$ $B^{\sharp}({\mathrm{d}}H)=\frac{\varphi}{L}{\frac{\partial}{\partial Q}}-\frac {Q}{C}{\frac{\partial}{\partial\varphi}}.$ $$

The vector field on ${\mathcal{M}}$ is $X=\dot{\varphi}{\frac{\partial}{\partial\varphi}}+\dot{Q}{\frac{\partial}{
\partial Q}}$.
Hamilton’s equations $X=B^{\sharp}({\mathrm{d}}H)$ lead to the equations of motion:

$$ $V_{L}=\dot{\varphi}=-\frac{Q}{C}=-V_{C}\text{ and }I_{C}=\dot{Q}=\frac{\varphi }{L}=I_{L}.$ $$

However, electrical circuits that can be described as Poisson systems are limited to energy-conservative LC circuits.
Our formulation extends this to include resistors, diodes, voltage sources, and current sources.

## Appendix C Dirac Structure

### C.1 Dirac Structure on a Vector Space

###### Proof of Theorem 1 (Yoshimura & Marsden, 2006a ) .

By Definition 1,

$$ $D_{V}^{\perp}=\{({\bm{w}},{\bm{\beta}})\in V\times V^{*}\mid\langle{\bm{\alpha }},{\bm{w}}\rangle+\langle{\bm{\beta}},{\bm{v}}\rangle=0\text{ for all }{\bm{v }}\in\Delta\text{ and }{\bm{\alpha}}\in\Delta^{\circ}\}.$ $$

Let $(\bar{\bm{v}},\bar{\bm{\alpha}})\in D_{V}$.
Then, $\bar{\bm{v}}\in\Delta$ and $\bar{\bm{\alpha}}\in\Delta^{\circ}$, so $\langle{\bm{\alpha}},\bar{\bm{v}}\rangle+\langle\bar{\bm{\alpha}},{\bm{v}}
\rangle=0$ for all $({\bm{v}},{\bm{\alpha}})\in D_{V}$.
This implies $(\bar{\bm{v}},\bar{\bm{\alpha}})\in D_{V}^{\perp}$
Thus, $D_{V}\subset D_{V}^{\perp}$.

Let $({\bm{w}},{\bm{\beta}})\in D_{V}^{\perp}$.
In the above definition of $D_{V}^{\perp}$, setting ${\bm{\alpha}}=0$ gives $\langle{\bm{\beta}},{\bm{v}}\rangle=0$ for all ${\bm{v}}\in\Delta$.
Hence, ${\bm{\beta}}\in\Delta^{\circ}$.
Similarly, setting ${\bm{v}}={\bm{0}}$ gives $\langle{\bm{\alpha}},{\bm{w}}\rangle=0$ for all ${\bm{\alpha}}\in\Delta^{\circ}$, which implies ${\bm{w}}\in\Delta$.
Hence, $({\bm{w}},{\bm{\beta}})\in D_{V}$.
Thus, $D_{V}^{\perp}\subset D_{V}$.

Therefore, $D_{V}=D_{V}^{\perp}$.
∎

###### Proof of Theorem 1 (Yoshimura & Marsden, 2006a ) .

### C.2 Dirac Structure on a Manifold

###### Definition 4 ( Courant ( 1990 ); Yoshimura & Marsden ( 2006a ) ) .

A distribution $\Delta$ on a manifold ${\mathcal{M}}$ is a collection of vector subspaces $\Delta_{\bm{u}}$ of tangent spaces ${T_{\bm{u}}{\mathcal{M}}}$, each assigned smoothly to ${\mathcal{M}}$ at point ${\bm{u}}$, forming a vector subbundle of ${T{\mathcal{M}}}$.
The annihilator $\Delta^{\circ}$ of $\Delta$ is a collection of the annihilator $\Delta^{\circ}_{\bm{u}}$ of $\Delta_{\bm{u}}$, also forming a subbundle of ${T^{*}{\mathcal{M}}}$.
Then, a Dirac structure $D$ is constructed as $D=\Delta\oplus\Delta^{\circ}$, which is a subbundle of the Pontryagin bundle ${T{\mathcal{M}}}\oplus{T^{*}{\mathcal{M}}}$.

A type of Dirac structure on a manifold ${\mathcal{M}}$ can be defined using a symplectic form $\Omega$ on ${\mathcal{M}}$.

###### Theorem 3 ( Yoshimura & Marsden ( 2006a ) ) .

Let $\Delta_{\mathcal{M}}$ be a distribution on a manifold ${\mathcal{M}}$.
Let $\Omega$ be a symplectic form on a manifold ${\mathcal{M}}$.
$\Omega$ is restricted to $\Delta_{\mathcal{M}}$ and denoted by $\Omega_{\Delta_{\mathcal{M}}}$.
Then, the collection of

$$ $\displaystyle(D_{\Delta_{\mathcal{M}}})_{\bm{u}}=\{({\bm{v}},{\bm{\alpha}})\in {T_{\bm{u}}{\mathcal{M}}}\times{T^{*}_{\bm{u}}{\mathcal{M}}}\mid$ $\displaystyle{\bm{v}}\in(\Delta_{\mathcal{M}})_{\bm{u}},$ $\displaystyle\text{ and }\langle{\bm{\alpha}},{\bm{w}}\rangle=(\Omega_{\Delta_ {\mathcal{M}}})_{\bm{u}}({\bm{v}},{\bm{w}})\text{ for all }{\bm{w}}\in(\Delta_ {\mathcal{M}})_{\bm{u}}\}$ $$

is a Dirac structure $D_{\Delta_{\mathcal{M}}}\subset{T{\mathcal{M}}}\oplus{T^{*}{\mathcal{M}}}$ on ${\mathcal{M}}$, which is said to be induced by $\Omega$.

###### Proof of Theorem 3 (Yoshimura & Marsden, 2006a ) .

By definition, the orthogonal of $D_{\Delta_{\mathcal{M}}}$ at ${\bm{u}}\in{\mathcal{M}}$ is given by

$$ $\displaystyle(D_{\Delta_{\mathcal{M}}}^{\perp})_{\bm{u}}=\{({\bm{w}},{\bm{ \beta}})$ $\displaystyle\in{T_{\bm{u}}{\mathcal{M}}}\times{T^{*}_{\bm{u}}{\mathcal{M}}} \mid\langle{\bm{\alpha}},{\bm{w}}\rangle+\langle{\bm{\beta}},{\bm{v}}\rangle=0 \text{ for all }{\bm{v}}\in(\Delta_{\mathcal{M}})_{\bm{u}},$ $\displaystyle\text{ and }\langle{\bm{\alpha}},{\bm{w}}\rangle=(\Omega_{\Delta_ {\mathcal{M}}})_{\bm{u}}({\bm{v}},{\bm{w}})\text{ for all }{\bm{w}}\in(\Delta_ {\mathcal{M}})_{\bm{u}}\}.$ $$

Let $({\bm{v}},{\bm{\alpha}}),(\bar{{\bm{v}}},\bar{{\bm{\alpha}}})\in(D_{\Delta_{
\mathcal{M}}})_{\bm{u}}$,
$\langle{\bm{\alpha}},\bar{{\bm{v}}}\rangle+\langle\bar{{\bm{\alpha}}},{\bm{v}}
\rangle=\Omega_{\Delta_{\mathcal{M}}}({\bm{v}},\bar{{\bm{v}}})+\Omega_{\Delta_
{\mathcal{M}}}(\bar{{\bm{v}}},{\bm{v}})=0$.
The latter equality holds because of the skew-symmetry of $\Omega$.
This implies $({\bm{v}},{\bm{\alpha}})\in(D_{\Delta_{\mathcal{M}}}^{\perp})_{\bm{u}}$.
Thus, $(D_{\Delta_{\mathcal{M}}})_{\bm{u}}\subset(D_{\Delta_{\mathcal{M}}}^{\perp})_{
\bm{u}}$.

Let $({\bm{w}},{\bm{\beta}})\in(D_{\Delta_{\mathcal{M}}}^{\perp})_{\bm{u}}$.
Then, $\langle{\bm{\alpha}},{\bm{w}}\rangle+\langle{\bm{\beta}},{\bm{v}}\rangle=0$ for all $({\bm{v}},{\bm{\alpha}})\in(D_{\Delta_{\mathcal{M}}})_{\bm{u}}$.
First, setting ${\bm{v}}={\bm{0}}$ gives $\langle{\bm{\alpha}},{\bm{w}}\rangle=(\Omega_{\Delta_{\mathcal{M}}})_{\bm{u}}(
{\bm{0}},{\bm{w}})=0$ for any ${\bm{\alpha}}\in(\Delta_{\mathcal{M}}^{\circ})_{\bm{u}}$.
Thus, ${\bm{w}}\in({\Delta_{\mathcal{M}}^{\circ}})^{\circ}_{\bm{u}}=(\Delta_{\mathcal
{M}})_{\bm{u}}$.
Second, if ${\bm{\alpha}}$ satisfies $\langle{\bm{\alpha}},{\bm{w}}\rangle=(\Omega_{\Delta_{\mathcal{M}}})_{\bm{u}}(
{\bm{v}},{\bm{w}})$ for any ${\bm{v}}\in(\Delta_{\mathcal{M}})_{\bm{u}}$, then $\langle{\bm{\alpha}},{\bm{w}}\rangle+\langle{\bm{\beta}},{\bm{v}}\rangle=(
\Omega_{\Delta_{\mathcal{M}}})_{\bm{u}}({\bm{v}},{\bm{w}})+\langle{\bm{\beta}}
,{\bm{v}}\rangle=0$ for any ${\bm{v}}\in(\Delta_{\mathcal{M}})_{\bm{u}}$.
This implies $\langle{\bm{\beta}},{\bm{v}}\rangle=(\Omega_{\Delta_{\mathcal{M}}})_{\bm{u}}({
\bm{w}},{\bm{v}})$ for any ${\bm{v}}\in(\Delta_{\mathcal{M}})_{\bm{u}}$.
Because it has been proved that ${\bm{w}}\in(\Delta_{\mathcal{M}})_{\bm{u}}$, $({\bm{w}},{\bm{\beta}})\in(D_{\Delta_{\mathcal{M}}})_{\bm{u}}$.
Thus, $(D_{\Delta_{\mathcal{M}}}^{\perp})_{\bm{u}}\subset(D_{\Delta_{\mathcal{M}}})_{
\bm{u}}$.

Therefore, $(D_{\Delta_{\mathcal{M}}})_{\bm{u}}=(D_{\Delta_{\mathcal{M}}}^{\perp})_{\bm{u}}$.
∎

###### Definition 5 (Hamilton-Dirac System (Yoshimura & Marsden, 2006b ) ) .

Let ${\mathcal{M}}$ be a manifold, $H:{\mathcal{M}}\to\mathbb{R}$ be an energy function, and $D$ be a Dirac structure on ${\mathcal{M}}$.
Given a vector field $X_{H}$ on ${\mathcal{M}}$, if it holds for each ${\bm{u}}(t)\in{\mathcal{M}}$ and $t\in[a,b]$ that

$$ $((X_{H})_{{\bm{u}}(t)},{\mathrm{d}}H_{{\bm{u}}(t)})\in(D_{\Delta_{\mathcal{M}} })_{{\bm{u}}(t)},$ $$

the tuple $(H,D_{\Delta_{\mathcal{M}}},X_{H})$ is called a Hamilton-Dirac system (or implicit Hamiltonian system).

Note that the symplectic form restricted to a distribution $\Delta_{\mathcal{M}}$, denoted by $\Omega_{\Delta_{\mathcal{M}}}$, satisfies $(\Omega_{\Delta_{\mathcal{M}}})_{\bm{u}}({\bm{v}},{\bm{w}})=\Omega_{\bm{u}}({
\bm{v}},{\bm{w}})$ for any ${\bm{v}},{\bm{w}}\in(\Delta_{\mathcal{M}})_{\bm{u}}$.
The symplectic form $\Omega$ is degenerate in the sense that its bundle map $\Omega^{\flat}$ on the tangent space ${T_{\bm{u}}{\mathcal{M}}}$ is degenerate, but $\Omega_{\Delta_{\mathcal{M}}}$ is non-degenerate on the distribution $(\Delta_{\mathcal{M}})_{\bm{u}}$.
If there is no constraint, $\Delta_{\mathcal{M}}={T{\mathcal{M}}}$ and $\Omega_{\Delta_{\mathcal{M}}}=\Omega$.
When the Dirac structure $D$ is given as in Theorem 3 with no constraint, the Hamilton-Dirac system is identical to the Hamiltonian system $\Omega^{\flat}(X_{H})={\mathrm{d}}H$ in Eq. (2).
The distribution $\Delta_{\mathcal{M}}$ describes constraints on the velocities, and hence constrained Hamiltonian systems can be rewritten as Hamilton-Dirac systems.
Therefore, Hamiltonian-Dirac systems are generalizations of Hamiltonian and constrained Hamiltonian systems.
For example, the system in Appendix A.2 can be written as a Hamilton-Dirac system as follows:

Degenerate Systems as Hamilton-Dirac Systems The distribution $\Delta_{\mathcal{M}}$ is given by $(\Delta_{\mathcal{M}})_{\bm{u}}=\{(\dot{q}_{1},\dot{q}_{2},\dot{p}_{1},\dot{p}
_{2})\in{T_{\bm{u}}{\mathcal{M}}}\mid\dot{q}_{1}-\dot{q}_{2}=0,\ \dot{p}_{1}/m
_{1}-\dot{p}_{2}/m_{2}=0\}=\operatorname{span}\{{\frac{\partial}{\partial q_{1
}}}+{\frac{\partial}{\partial q_{2}}},\frac{1}{m_{1}}{\frac{\partial}{\partial
p
_{1}}}+\frac{1}{m_{2}}{\frac{\partial}{\partial p_{2}}}\}$.
Define the new coordinates $q=q_{1}=q_{2}$ and $p=\frac{m_{1}+m_{2}}{m_{1}}p_{1}=\frac{m_{1}+m_{2}}{m_{2}}p_{2}$, with $\Delta_{\mathcal{M}}=\operatorname{span}\{{\frac{\partial}{\partial q}},{\frac
{\partial}{\partial p}}\}$.
The symplectic form $\Omega_{\Delta_{\mathcal{M}}}$ restricted to the distribution $\Delta_{\mathcal{M}}$ should satisfy $(\Omega_{\Delta_{\mathcal{M}}})_{\bm{u}}({\bm{v}},{\bm{w}})=\Omega_{\bm{u}}({
\bm{v}},{\bm{w}})$ for any ${\bm{v}},{\bm{w}}\in(\Delta_{\mathcal{M}})_{\bm{u}}$ at ${\bm{u}}$.
Such form is given as $\Omega_{\Delta_{\mathcal{M}}}={\mathrm{d}}q\wedge{\mathrm{d}}p$.
Then, the equations of motion are written as $\dot{q}=\frac{p}{m_{1}+m_{2}},\ \dot{p}=-(k_{1}+k_{2})q$.

###### Definition 4 ( Courant ( 1990 ); Yoshimura & Marsden ( 2006a ) ) .

###### Theorem 3 ( Yoshimura & Marsden ( 2006a ) ) .

###### Proof of Theorem 3 (Yoshimura & Marsden, 2006a ) .

###### Definition 5 (Hamilton-Dirac System (Yoshimura & Marsden, 2006b ) ) .

### C.3 Port-based Systems

###### Proof of Theorem 2 .

By definition, the orthogonal of $D$ at ${\bm{u}}\in{\mathcal{M}}$ is given by

$$ $D^{\perp}_{\bm{u}}=\{(\bar{\bm{f}},\bar{\bm{e}})\in{\mathcal{F}}_{\bm{u}} \times{\mathcal{E}}_{\bm{u}}\mid\langle{\bm{e}},\bar{\bm{f}}\rangle+\langle \bar{\bm{e}},{\bm{f}}\rangle=0\text{ for all }{\bm{e}}\in{\mathcal{E}}_{\bm{u} }\text{ and }{\bm{f}}=B^{\sharp}_{\bm{u}}({\bm{e}})\}.$ $$

Let $({\bm{f}},{\bm{e}}),(\bar{{\bm{f}}},\bar{{\bm{e}}})\in D_{\bm{u}}$.
Then, $\langle{\bm{e}},\bar{{\bm{f}}}\rangle+\langle\bar{{\bm{e}}},{\bm{f}}\rangle=
\langle{\bm{e}},B^{\sharp}_{\bm{u}}(\bar{\bm{e}})\rangle+\langle\bar{{\bm{e}}}
,B^{\sharp}_{\bm{u}}({\bm{e}})\rangle=B_{\bm{u}}({\bm{e}},\bar{\bm{e}})+B_{\bm
{u}}(\bar{\bm{e}},{\bm{e}})=0$ due to the skew-symmetry of $B$.
This implies $({\bm{f}},{\bm{e}})\in D_{\bm{u}}^{\perp}$.
Thus, $D_{\bm{u}}\subset D^{\perp}_{\bm{u}}$.

Let $(\bar{\bm{f}},\bar{\bm{e}})\in D^{\perp}_{\bm{u}}$.
Then, for any $({\bm{f}},{\bm{e}})\in D_{\bm{u}}$ (i.e., for any ${\bm{e}}\in{\mathcal{E}}_{\bm{u}}$ and ${\bm{f}}=B^{\sharp}_{\bm{u}}({\bm{e}})$), $0=\langle{\bm{e}},\bar{\bm{f}}\rangle+\langle\bar{\bm{e}},{\bm{f}}\rangle=
\langle{\bm{e}},\bar{\bm{f}}\rangle+\langle\bar{\bm{e}},B^{\sharp}_{\bm{u}}({
\bm{e}})\rangle=\langle{\bm{e}},\bar{\bm{f}}\rangle+B_{\bm{u}}(\bar{\bm{e}},{
\bm{e}})=\langle{\bm{e}},\bar{\bm{f}}\rangle-B_{\bm{u}}({\bm{e}},\bar{\bm{e}})
=\langle{\bm{e}},\bar{\bm{f}}\rangle-\langle{\bm{e}},B^{\sharp}_{\bm{u}}(\bar{
\bm{e}})\rangle$.
$\langle{\bm{e}},\bar{\bm{f}}\rangle=B_{\bm{u}}({\bm{e}},\bar{\bm{e}})$ for all ${\bm{e}}\in{\mathcal{E}}_{\bm{u}}$ implies that $\bar{\bm{f}}=B^{\sharp}_{\bm{u}}(\bar{\bm{e}})$.
Thus, $D^{\perp}_{\bm{u}}\subset D_{\bm{u}}$.

Therefore, $D_{\bm{u}}=D^{\perp}_{\bm{u}}$.
∎

By extending the above Hamilton-Dirac system from the Pontryagin bundle ${T{\mathcal{M}}}\oplus{T^{*}{\mathcal{M}}}$ to the vector bundle ${\mathcal{F}}\oplus{\mathcal{E}}$, we obtain the so-called port-Hamiltonian systems (Courant, ; van der Schaft & Jeltsema, ).
However, previous neural network-based methods employed the canonical form in Eq. (4) and have not attempted to identify coupling patterns or component-wise characteristics (Zhong et al., ; Eidnes et al., ; Neary & Topcu, ; Di Persio et al., ).

In a bond graph representation of a dynamical system, a bond represents a component, its ports (flow and effort) represent points of interaction with other bonds, and the arrows connected to the ports represent the interactions between them.
Both the Port-Hamiltonian and our formulations, as well as their terminology, are based on this bond graph structure (Duindam et al., ).

Instead of the symplectic form $\Omega$, our formulation employs (Poisson) bivector $B$ to define the Dirac structure $D$ on the tangent and cotangent bundles (${T_{\bm{u}}{\mathcal{M}}}={\mathcal{F}}^{S}$ and ${T^{*}_{\bm{u}}{\mathcal{M}}}={\mathcal{E}}^{S}$) plus the vector bundles for port variables (${\mathcal{F}}^{R}$, ${\mathcal{F}}^{I}$, ${\mathcal{E}}^{R}$, and ${\mathcal{E}}^{I}$).
As is the case with the symplectic form $\Omega$, we can constrain a Poisson bivector $B$ on a codistribution $\Delta^{*}_{\mathcal{M}}$, which is a subbundle of the cotangent bundle ${T^{*}{\mathcal{M}}}$, and then define the Dirac structure.
However, because the bivector inherently handles constraints, we did not adopt this formulation.
Thus, our formulation is a special case of the Poisson-Dirac formulation with ports, as summarized in Table 1.
Note that, for a bivector to be called a Poisson bivector, it must satisfy certain conditions, such as the Jacobi identity (Courant, ).
However, the learned bivector in our formulation does not necessarily satisfy these conditions, nor are these conditions required for the proofs presented earlier.
Thus, throughout this paper, we refer to $B$ simply as a bivector, without restricting it to being a Poisson bivector.

###### Proof of Theorem 2 .

## Appendix D Datasets and Experimental Settings

### D.1 Implementation Details

Following previous studies (Greydanus et al., ; Matsubara et al., ), we used fully-connected neural networks with two hidden layers to implement any vector fields and energy functions for all models.
Each hidden layer had 200 units, followed by a hyperbolic tangent activation function.
Weight matrices were initialized using PyTorch’s default algorithm.

In PoDiNNs, each element of the bivector $B$ related to energy-dissipating components was initialized from a uniform distribution $U(-0.1,0.1)$, while the remaining elements are set to zero.
These bivector elements were updated along with the neural network parameters.
Elements representing incompatible couplings are constrained to be zero.
For instance, masses cannot couple directly with each other but can couple with springs, dampers, or external forces.

Unless stated otherwise, the time step size was set to $\Delta t=0.1$.
Each training subset consisted of 1,000 trajectories of 1,000 steps, and each test subset consisted of 10 trajectories of 10,000 steps.
The Dormand–Prince method (dopri5) with absolute tolerance $\operatorname{atol}=10^{-7}$ and relative tolerance $\operatorname{rtol}=10^{-9}$ was used to integrate the ground truth ODEs and neural network models (Dormand & Prince, ).
The Adam optimization algorithm (Kingma & Ba, ) was applied with parameters $(\beta_{1},\beta_{2})=(0.9,0.999)$ and a batch size of 100 for updates.
The learning rate was initialized at $10^{-3}$ and decayed to zero using cosine annealing (Loshchilov & Hutter, ).
The number of training iterations was set to 100,000.

All experiments were conducted on a single NVIDIA A100 GPU.

### D.2 Mechanical Systems

Overview and Experimental Setting A spring $k$ generates a force that restores it to its original length when stretched or compressed.
Specifically, this force is given by $e^{S}=-k(\Delta q)$, where $\Delta q$ is the spring’s displacement.
In both the Hamiltonian and our formulations, the potential energy $U$ is obtained by integrating this force $e^{S}$ with respect to the displacement $\Delta q$.
Conversely, the force $e^{S}$ is obtained as to the partial derivative of the potential energy $U$ with respect to the displacement $\Delta q$.
The spring’s flow $f^{S}$ is the rate of extension, $\Delta\dot{q}$.

A mass $m$ has a velocity $v$, and its momentum is given by $p=mv$.
The kinetic energy is $\frac{1}{2}mv^{2}=\frac{1}{2m}p^{2}$.
In general, velocity $v$ is more easily observed than momentum $p$ as the state of a mass.
Therefore, we provided velocity $v$ as the observation.

Coupling between two components means the output (effort) of one flows into the input (flow) of the other.
Hence, possible couplings are limited to interactions between potential and kinetic components, as shown in Table 2.
Specifically, the following couplings are possible: mass and spring ${\frac{\partial}{\partial p_{i}}}\wedge{\frac{\partial}{\partial q_{j}}}$, mass and damper ${\frac{\partial}{\partial p_{i}}}\wedge\xi^{R}_{j}$, mass and external force ${\frac{\partial}{\partial p_{i}}}\wedge\xi^{I}_{j}$, moving wall and spring $\xi_{i}^{I}\wedge{\frac{\partial}{\partial q_{j}}}$, and moving wall and damper $\xi_{i}^{I}\wedge\xi^{R}_{j}$.
These combinations were incorporated as elements of the learnable bivector $B$, while all other possible pairings were fixed to zero.

In the absolute coordinate system, the positions $q_{i}$ are used to represent the states of the springs, and their displacements $\Delta q_{i}$ are then calculated within the potential energy function $U$.
Therefore, the potential energy function $U$ depends on all positions $q_{i}$ of springs and moving boundaries collectively.
In the relative coordinate system, the displacements $\Delta q_{i}$ are used instead, and in this case, the potential energy function $U_{i}$ is defined individually for each displacement $\Delta q_{i}$.

In both Dissipative SymODENs and PoDiNNs, since we assume to know the symbolic expression $\frac{1}{2m}p^{2}$ for the kinetic energy, this was used, rather than a neural network.
Also, the velocity $v$ was transformed into momentum $p$ using the learnable parameter $m$ with $p=mv$.
This approach was originally adopted by Dissipative SymODENs and is a realistic and practical choice.

For Dissipative SymODENs and PoDiNNs in the absolute coordinate system, a single neural network was trained to approximate the potential energy function $U$ using all positions $q_{i}$ of springs and moving boundaries as inputs.
In the relative coordinate system, a separate neural network $U_{i}$ was used for the potential energy of each spring $k_{i}$.

For Dissipative SymODENs, the dissipative terms $D$ and input gain $G$ were modeled using neural networks.
In the original paper (Zhong et al., ), $D$ was defined as a symmetric matrix that depended solely on the positions $q$ of the springs, which is suitable for modeling linear dampers.
To extend this to nonlinear dampers, we modified $D$ to also depend on the velocity $v$ of the masses.

(a) Mass-Spring-Damper System with External Force This system consists of three springs $k_{i}$ and three masses $m_{i}$ arranged sequentially and indexed by $i\in\{1,2,3\}$ from a fixed wall.
Two dampers $d_{1}$ and $d_{3}$ are placed in parallel with springs $k_{1}$ and $k_{3}$, respectively.
An external force is applied to mass $m_{3}$.

The masses were set to $m_{i}=0.8+0.2i$.
The characteristics of the nonlinear spring were given by $k_{i}(\Delta q_{i})=(0.1+0.1i)\Delta q_{i}+0.1\Delta q_{i}^{3}$.
The characteristics of the nonlinear dampers were defined as $d_{i}(v_{i})=(0.1-0.02i)\,\operatorname{sgn}(v_{i})|v_{i}|^{1/3}$, where $v_{i}$ denotes the extension velocity of the $i$-th damper, and $\operatorname{sgn}$ is the sign function that returns $1$ for a positive value and $-1$ for a negative value.
The initial positions of the springs and the initial velocities of the masses were sampled from the uniform distributions $U(-0.5,0.5)$ and $U(-0.3,0.3)$, respectively.
The external force $e^{I}(t)$ was defined as the sum of three sine waves, with each wave’s amplitude, angular velocity, and initial phase sampled from the uniform distributions $U(0.2,0.5)$, $U(0.1\pi,0.2\pi)$, and $U(0,2\pi)$, respectively.

(b) Mass-Spring-Damper System with Moving Boundary This system consists of three springs $k_{i}$ and three masses $m_{i}$ arranged sequentially and indexed by $i\in\{1,2,3\}$ from a moving wall $b$.
Three dampers $d_{i}$ are placed in parallel with springs $k_{i}$.
This potentially represents a building’s response during an earthquake.

The masses were set to $m_{i}=1.6-0.2i$.
The characteristics of nonlinear spring were given by $k_{i}=(0.6-0.1i)\Delta q_{i}+0.1\Delta q_{i}^{3}$, where $\Delta q_{i}$ denotes the displacement of the $i$-th spring.
The characteristics of nonlinear damper were defined as $d_{i}(v_{i})=\tilde{d}_{i}\operatorname{sgn}(v_{i})|v_{i}|^{1/3}$ for $\tilde{d}_{1}=0.10$, $\tilde{d}_{2}=0.05$, and $\tilde{d}_{3}=0.02$.
The initial positions of springs $k_{i}$ and the initial velocities of masses $m_{i}$ were sampled from the uniform distributions $U(-0.5,0.5)$ and $U(-0.3,0.3)$, respectively.
The position $q_{b}$ of moving wall $b$ was set as the sum of three sine waves, with each wave’s amplitude, angular velocity, and initial phase sampled from the uniform distributions $U(0.2,0.4)$, $U(0.05\pi,0.2\pi)$, and $U(0,2\pi)$, respectively.

In the absolute coordinate system, it is necessary to represent the potential energy of spring $k_{1}$ connected to the moving wall $b$.
For both Dissipative SymODENs and PoDiNNs, the position $q_{b}$ of the moving wall $b$ was included as part of the inputs to the neural network $U$ for the potential energy, along with the positions $q_{i}$ of all springs $k_{i}$.
It was also fed to Neural ODEs.

The effort $e^{I}$ of the moving wall $b$ was its velocity $v_{b}$, which was fed to PoDiNNs and Neural ODEs as part of the external inputs.
However, no such mechanism exists for Dissipative SymODENs.
Without this, models cannot represent the force generated by damper $d_{1}$, connected to moving wall $b$.

(c) Mass-Spring System with Redundancy This system consists of five springs $k_{i}$ and two masses $m_{i}$ arranged in 2-dimensional space, as shown in Fig. 2 (c).
Masses $m_{1}$ and $m_{2}$ are connected to a fixed wall via springs $k_{1}$ and $k_{2}$, respectively.
These masses are also connected to each other by spring $k_{3}$.
Additionally, springs $k_{4}$ and $k_{5}$ diagonally connect masses $m_{2}$ and $m_{1}$ to the fixed wall, respectively, similar to cross braces.
Since there are no energy-dissipating components or external inputs, the total energy is conserved.
The masses were set to $m_{1}=5.0$ and $m_{2}=3.0$.
The natural lengths of spring $k_{1}$ and $k_{2}$ to $l_{1}=3.0$, that of $k_{3}$ to $l_{3}=4.0$, and those of $k_{4}$ and $k_{5}$ to $l_{4}=5.0$.
The characteristics of nonlinear spring were $k_{1}(\Delta q)=2.5\Delta q+3.4\Delta q^{3}$, $k_{2}(\Delta q)=3.0\Delta q+0.5\Delta q^{3}$, $k_{3}(\Delta q)=2.1\Delta q+4.1\Delta q^{3}$, $k_{4}(\Delta q)=3.5\Delta q+2.4\Delta q^{3}$, $k_{5}(\Delta q)=2.5\Delta q+1.6\Delta q^{3}$ for the displacement $\Delta q$.
The initial positions of masses $m_{1}$ and $m_{2}$ were sampled from uniform distributions $U(-0.5,0.5)$ in both $x$ and $y$ directions, and their velocities from $U(-0.1,0.1)$.

Because of two masses in 2-dimensional space, the system has 4 degrees of freedom for configuration, or 8 when including velocities.
However, the observations were composed of the velocities of both masses and the displacements of all five springs, resulting in a 14-dimensional observation space ${\mathcal{M}}$.
Hence, the dynamics is degenerate.

Due to the degeneracy, HNNs are not applicable (Greydanus et al., ).
Also, since this degeneracy does not come from a holonomic constraint, CHNNs are also not applicable (Finzi et al., ).
In PNNs, we used a Real-NVP consisting of four coupling layers for the coordinate transformation.
Each coupling layer was made of a fully-connected neural network with two hidden layers of 200 units, followed by a hyperbolic tangent activation function.
The 14-dimensional observations were transformed, and the eight dimensions were extracted (with the remaining six dimensions considered constant) as input to the energy function of HNNs.

### D.3 Electrical Systems

Overview and Experimental Settings Refer to Appendix B.2 for basic characteristics of capacitors and inductors.
Since the electric charge is analogous to displacement, a capacitor corresponds to a spring, while an inductor corresponds to a mass in mechanical systems.

The states of a capacitor $C$ and inductor $L$ are the electric charge $Q$ and magnetic flux $\varphi$, respectively, but these are difficult to observe directly.
We used the capacitor voltage $V_{C}=Q/C$ and inductor current $I_{L}=\varphi/L$ as the observations, which linearly correlate with their states.
PoDiNNs learned the element characteristics $C$ and $L$ as learnable parameters and internally performed the transformations $Q=CV_{C}$ and $\varphi=LI_{L}$.
This is the same approach used by Dissipative SymODENs for masses.
Strictly speaking, since PoDiNNs can learn coordinate transformations, they can approximate the dynamics even without such transformations.

In a mass-spring-damper system, the damper always has velocity as flow and force as effort, but the flow of a resistor can be either current or voltage.
This depends on the overall coupling pattern and the formulation step, though the total number remains constant.
If two resistors with different types of flow are coupled, PoDiNNs would result in a DAE, which cannot be solved explicitly.
Additionally, components that share the same effort type, such capacitors and current voltage sources, may couple.
In this case, an infinite current instantaneously flows into the capacitor, and its voltage matches that of the direct voltage source.
This dynamics also requires a DAE, rather than an ODE.
In this study, we assumed that the system does not include couplings requiring DAEs.
Methods capable of handling such circuits will be a topic for future research.

A separate neural network was used for approximating the characteristics of each energy-storing and dissipating component.

(d) FitzHugh-Nagumo Model FitzHugh-Nagumo model is a model for the electrical dynamics of a biological neuron, exhibiting oscillatory behavior when the external current $J$ is applied (Izhikevich & FitzHugh, ).
The governing equations are written as:

$$ $\displaystyle\dot{V}$ $\displaystyle=V-\frac{V^{3}}{3}-W+J,$ (A1) $\displaystyle\dot{W}$ $\displaystyle=0.08(V+0.7-0.8W),$ $$

where $V$ denotes the membrane potential, $W$ is the recovery variable, and $J$ is the input.

A circuit representation consists of a resistor $R_{2}$, an inductor $L$, a capacitor $C$, a current voltage source $E$, and a tunnel diode $R_{1}$, with an external current $J$.
The characteristics of these elements are defined as $L=1/0.08$, $R_{2}=0.8$, $C=1.0$, and $E=-0.7$.
The current $I_{R_{1}}$ through tunnel diode $R_{1}$ is characterized by $I_{R_{1}}=D(V_{R_{1}})=V_{R_{1}}^{3}/3-V_{R_{1}}$, where $V_{R_{1}}$ is the voltage across the diode.
The membrane potential $V$ and recovery variable $W$ are represented by the capacitor voltage $V_{C}$ and the inductor current $I_{L}$.
In our formulation, resistor $R_{2}$ and diode $R_{1}$ have the current and voltage as their flows, respectively.
Also, because the voltage generated by the current voltage source $E$ is unchanged for any trials, it can be treated as a part of resistor $R_{2}$.

The initial values of $V$ and $W$ were sampled from the uniform distribution $U(-3.0,3.0)$.
The external current $J$ was sampled at evenly spaced intervals within the range of 0.1 to 1.5 for each trajectory, and it was kept constant within each individual trajectory to evaluate whether each model can learn the oscillatory behavior.

Since the external current $J$ was kept constant, the variability of the dynamics was reduced, making the learning process easier.
To preserve the challenge of the task, only 30 trajectories were generated for the training subset.

(e) Chua’s Circuit Chua’s circuit is a nonlinear electronic circuit known for its chaotic behavior (Chua, ).
It consists of linear elements (a resistor $R_{1}$, two capacitors $C_{1}$ and $C_{2}$, and an inductor $L$) and a nonlinear element $R_{2}$, known as Chua’s diode.
With $R_{1}=1$, the governing equations are:

$$ $\displaystyle\dot{V}_{C_{1}}$ $\displaystyle=\alpha(V_{C_{2}}-V_{C_{1}}-f(V_{R_{2}})),$ (A2) $\displaystyle\dot{V}_{C_{2}}$ $\displaystyle=V_{C_{1}}-V_{C_{2}}+I_{L},$ $\displaystyle\dot{I}_{L}$ $\displaystyle=-\beta V_{C_{2}},$ $\displaystyle V_{R_{2}}$ $\displaystyle=V_{C_{1}}.$ $$

The parameters are $\alpha=1/C$ and $\beta=1/L$.
The nonlinear function $f(V_{R_{2}})$ describes the voltage-current characteristic of Chua’s diode:

$$ $f(V_{R_{2}})=m_{1}V_{R_{2}}+0.5(m_{0}-m_{1})(|V_{R_{2}}+1|-|V_{R_{2}}-1|),$ (A3) $$

where $V_{R_{2}}$ is the voltage across the diode, equal to $V_{C_{1}}$, and $m_{0}$ and $m_{1}$ are parameters.
In our formulation, both resistor $R_{1}$ and diode $R_{2}$ have voltages as flows.

We set $\alpha=15.6$, $\beta=28$, $m_{0}=-8/7$, and $m_{1}=-5/7$.
The initial values of $V_{C_{1}}$, $V_{C_{2}}$, and $I_{L}$ were sampled from the uniform distribution $U(-0.5,0.5)$.

Due to its chaotic behavior, we set the time step size to $\Delta t=0.01$ and the number of training iterations to 1,000,000.

### D.4 Multiphysics

Overview A multiphysics system involves components from different domains that are coupled together.
Our formulation inherently handles multiphysics systems as long as careful attention is paid to which components can and cannot be coupled.

(f) DC Motor

In this system, a DC motor bridges the electro-magnetic and the rotational domains.
In the electric circuit, an inductor $L$ and a resistor $R$ represent the inductance and resistance of the motor’s armature winding, respectively.
Additionally, a voltage source $E$ serves as an external input.
A massless pendulum rod of length $l$ is attached to the DC motor.
The pendulum’s angle $\theta$ is measured from the vertical position, and its angular velocity is denoted by $\omega$.
With a mass $m$ at the rod’s end and the gravitational acceleration $g$, the pendulum has a potential energy of $-mgl\cos\theta$ and experiences a torque of $-mgl\sin\theta$.
Additionally, friction $d$ occurs at the pivot point of the pendulum.

The DC motor generates torque $\tau_{DC}$ based on the current $I_{DC}$ through the armature.
We assume a linear relationship with constant $K$, such that $\tau_{DC}=KI_{DC}$.
Conversely, when the motor rotates at an angular velocity $\omega$, it produces a back electromotive force given by $V_{DC}=-K\omega$.
The governing equations are:

$$ $\displaystyle\dot{\theta}$ $\displaystyle=\omega,$ $\displaystyle ml^{2}\dot{\omega}$ $\displaystyle=-mgl\sin\theta+KI-d(\omega),$ $\displaystyle L\dot{I}$ $\displaystyle=-\omega K+E-R(I).$ $$

We set $L=2.5$, $m=2.0$, $l=1.5$, $g=1.0$, and $K=0.5$, respectively.
Also, we set $d(\omega)=0.02\operatorname{sgn}(\omega)|\omega|^{1/3}$ for friction $d$, and $R(I)=0.05I^{3}$ for resistor $R$.

PoDiNNs can model this multiphysics system seamlessly.
In the magnetic domain, inductor $L$ has magnetic flux $\varphi=LI$ as its state, voltage as its flow, and current as its effort.
For the pendulum’s motion in the rotational domain, the state is angular momentum $p=ml^{2}\omega$, the flow is torque, and the effort is angular velocity (see Table 2).
Because the properties of the armature are represented by inductor $L$ and resistor $R$, the remaining function of the DC motor can be represented by a coupling between inductor $L$ and the pendulum’s motion $p$, that is, the bivector element $K{\frac{\partial}{\partial\varphi}}\wedge{\frac{\partial}{\partial p}}$.
Both inductor $L$ and mass $m$ are considered to store kinetic energies, and this type of coupling is referred to as a gyrator.

Therefore, the only notable point for the implementation is that magnetic flux $\varphi$ and angular momentum $p$, which belong to different domains, can be coupled.
Since there are only a few elements in each domain, the coupling patterns are unique, and their strengths were set to 1.0 while keeping the coefficient $K$ of the bivector element ${\frac{\partial}{\partial\varphi}}\wedge{\frac{\partial}{\partial p}}$ learnable.

(g) Hydraulic Tank Consider a hydraulic tank storing incompressible fluid, which belongs to the hydraulic domain and can couple with components in the mechanical domain.
Let $V$ denote the volume of the fluid inside the tank.
The tank has a cross-sectional area $A$ and a fluid height $h$, giving the relationship $V=Ah$.
Let $\rho$ represent the density of the fluid per unit volume, and $g$ the gravitational acceleration.
The pressure $p$ exerted on the bottom of the tank per unit area is given by $p=\rho gh=\frac{\rho gV}{A}$, and the total force acting on the bottom surface is $\rho ghA$.

Assuming that fluid is supplied from the bottom of the tank, the potential energy stored in the tank can be expressed as:

$$ $U_{V}=\int_{0}^{h}\rho ghA\,{\mathrm{d}}h=\frac{1}{2}\rho gh^{2}A=\frac{\rho g }{2A}V^{2}.$ $$

Then, ${\frac{\partial U_{V}}{\partial V}}=\rho gV/A=p$.

Consider two cylinders attached in opposite directions at the bottom of the tank, with cross-sectional areas $a_{1}$ and $a_{2}$.
Each cylinder contains a piston, with masses $m_{1}$ and $m_{2}$, respectively.
The forces acting on the pistons are $-pa_{1}$ and $pa_{2}$, where the positive direction is towards $m_{2}$.
When piston $m_{1}$ moves by a displacement $q_{1}$, a volume of fluid $a_{1}q_{1}$ flows into the tank.
Similarly, when piston $m_{2}$ moves by $q_{2}$, a volume $a_{2}q_{2}$ flows out of the tank.
For simplicity, we assume that these inflows and outflows occur adiabatically, without any resistance.
However, compressible fluids, non-adiabatic process, fluid momentum, and fluid resistance could also be incorporated (Duindam et al., ).

Each piston $m_{i}$ moves with a velocity $v_{i}$, and is connected to a fixed wall via a spring $k_{i}$ and a damper $d_{i}$ for each $i$.
An external force $F$ is applied to piston $m_{2}$.

The equations of motion for the system can be written as:

$$ $\displaystyle\dot{V}$ $\displaystyle=a_{1}v_{1}-a_{2}v_{2},$ $\displaystyle m_{1}\dot{v}_{1}$ $\displaystyle=-pa_{1}-k(q_{1})-d_{1}(v_{1}),$ $\displaystyle m_{2}\dot{v}_{2}$ $\displaystyle=pa_{2}-k(q_{2})-d_{2}(v_{2}),$ $$

Each piston has kinetic energy, and each spring has potential energy.
The state of each piston can be expressed as momentum, $p_{i}=m_{i}v_{i}$, and the state of the spring can be described by the displacement $q_{i}$ for each $i$.
Thus, the total energy of the system is given by:

$$ $H=\sum_{i=1}^{2}\frac{p_{i}^{2}}{2m_{i}}+\sum_{i=1}^{2}U_{i}(q_{i})+\frac{\rho g }{2A}V^{2},$ $$

where $U_{i}$ denotes the potential energy stored in spring $k_{i}$, as well as the characteristics of dampers $d_{i}$.

The above equations of motion can be expressed using the energy function $H$ with the following bivector $B$:

$$ $\displaystyle B=$ $\displaystyle a_{1}\left({\frac{\partial}{\partial p_{1}}}\wedge{\frac{ \partial}{\partial V}}\right)-a_{2}\left({\frac{\partial}{\partial p_{2}}} \wedge{\frac{\partial}{\partial V}}\right)+{\frac{\partial}{\partial p_{1}}} \wedge{\frac{\partial}{\partial q_{1}}}+{\frac{\partial}{\partial p_{2}}} \wedge{\frac{\partial}{\partial q_{2}}}$ $\displaystyle+{\frac{\partial}{\partial p_{1}}}\wedge\xi^{R}_{1}+{\frac{ \partial}{\partial p_{2}}}\wedge\xi^{R}_{2}-{\frac{\partial}{\partial p_{2}}} \wedge\xi^{I},$ $$

where $\xi^{R}_{i}$ and $\xi^{I}$ represent the basis vectors of the spaces ${\mathcal{F}}^{R}$ and ${\mathcal{F}}^{I}$ for dampers $d_{i}$ and external force $F$, respectively.

From this, it seems that the hydraulic tank behaves similarly to a spring, but the key difference is that the coupling strength depends on the cross-sectional areas $a_{1}$ and $a_{2}$ of the cylinders.
This type of coupling is referred to as a transformer.
If the cross-sectional area of the tank is not constant, the tank would behave like a nonlinear spring.

We set the parameters as follows: $A=5.0$, $g=1.0$, $\rho=10.0$, $a_{1}=1.0$, $a_{2}=0.3$, $m_{1}=3.0$, and $m_{2}=1.0$.
The spring forces were defined as $k_{i}(q_{i})=0.1q_{i}+0.01q_{i}^{3}$, and the damping forces as $d_{i}(v_{i})=(0.1-0.04i)\,\operatorname{sgn}(v_{i})|v_{i}|^{1/3}$ for $i=1,2$.

As the system tends to reach equilibrium around $V=5$, $q_{1}=-10$, and $q_{2}=6$, the initial fluid volume $V$ and the springs’ displacements $q_{1}$ and $q_{2}$ were sampled from the uniform distributions $U(5-0.25,5+0.25)$, $U(-10-0.3,-10+0.3)$, and $U(6-0.3,6+0.3)$, respectively,
The initial velocities of the pistons were sampled from the uniform distribution $U(-0.3,0.3)$.
The external force $e^{I}(t)$ was defined as the sum of three sine waves, where the amplitude, angular velocity, and initial phase of each wave were sampled from the uniform distributions $U(0.05,0.2)$, $U(0.1\pi,0.3\pi)$, and $U(0,2\pi)$, respectively.

In Fig. 3, for clarity, we displayed the changes relative to $5$, $-10$, and $6$, rather than the actual values of $V$, $q_{1}$, and $q_{2}$.

## Appendix E Additional Results and Discussions

### E.1 Additional Results

Figure: Figure A1: Visualizations of example results. Each top panel shows ground truth trajectories, while the other panels show the absolute errors of all 10 trials in semi-transparent color. See also Fig. 3.
Refer to caption: x14.png

Additional Visualization We show example visualizations of trajectories and absolute errors in Fig. A1 in addition to Fig. 3.
Across all datasets, PoDiNNs demonstrate superior predictive performance compared to other methods.
However, in both systems (b) and (b’), there was one failure.
This was due to minor inaccuracies in the learned bivectors, as discussed later.
While this issue only occurred twice out of 90 trials, it indicates that there is room for further improvement of the initialization and regularization of bivector elements.

Identifying Coupling Patterns of System (b) We examined how the bivector $B$ identifies the coupling patterns using system (b), that is, the mass-spring-damper system with moving boundary.
When two coefficients differ by a factor of 1,000 or more, we considered the larger one as a detected coupling and the smaller one as effectively zero, indicating no coupling.
The basis vector ${\frac{\partial}{\partial q_{i}}}$ corresponds to spring $k_{i}$, ${\frac{\partial}{\partial p_{i}}}$ to mass $m_{i}$, and $\xi^{S}$ to moving boundary $b$.

In the relative coordinate system, we obtained the bivector
$B=({\frac{\partial}{\partial p_{1}}}-\xi^{S})\wedge({\frac{\partial}{\partial q
_{1}}}+a_{1}\xi^{R}_{1})+({\frac{\partial}{\partial p_{2}}}-{\frac{\partial}{
\partial p_{1}}})\wedge({\frac{\partial}{\partial q_{2}}}+a_{2}\xi^{R}_{2})+({
\frac{\partial}{\partial p_{3}}}-{\frac{\partial}{\partial p_{2}}})\wedge({
\frac{\partial}{\partial q_{3}}}+a_{3}\xi^{R}_{3})$
in all 10 trials.
Here, the basis vectors $\xi^{R}_{1}$, $\xi^{R}_{2}$, and $\xi^{R}_{3}$ for dampers $d_{1}$, $d_{2}$, and $d_{3}$ were appropriately reordered because their indices are interchangeable.
$a_{1}$, $a_{2}$, and $a_{3}$ are trial-wise positive parameters.
This result matches the ground truth coupling pattern in Fig. 2 (b).
Notably, the coefficients of other bivector elements, such as $\xi^{S}\wedge{\frac{\partial}{\partial q_{3}}}$ and ${\frac{\partial}{\partial p_{3}}}\wedge{\frac{\partial}{\partial q_{2}}}$, were effectively zero, indicating that PoDiNNs correctly identified the absence of non-existent couplings.

The coefficients expected to be $1$ or $-1$ were accurate to five significant figures in 9 out of 10 trials.
However, in one trial, there was a small error of approximately 0.002.
This trial corresponded with the failure case shown in Fig. A1.
This suggests that even minor errors in the coefficients of the bivector elements can lead to significant prediction failures.
Nonetheless, as mentioned earlier, such occurrences are extremely rare.

Identifying Coupling Patterns of System (b’) In the absolute coordinate system, we obtained the bivector
$B={\frac{\partial}{\partial p_{1}}}\wedge{\frac{\partial}{\partial q_{1}}}+{
\frac{\partial}{\partial p_{2}}}\wedge{\frac{\partial}{\partial q_{2}}}+{\frac
{\partial}{\partial p_{3}}}\wedge{\frac{\partial}{\partial q_{3}}}+a_{1}({
\frac{\partial}{\partial p_{1}}}-\xi^{S})\wedge\xi^{R}_{1}+a_{2}({\frac{
\partial}{\partial p_{2}}}-{\frac{\partial}{\partial p_{1}}})\wedge\xi^{R}_{2}
+a_{3}({\frac{\partial}{\partial p_{3}}}-{\frac{\partial}{\partial p_{2}}})
\wedge\xi^{R}_{3}$
in 9 out of 10 trials, where coefficients with a difference of less than 1% were considered identical.
The relationship between the masses and the springs is represented by the standard Poisson bivector $\sum_{i}{\frac{\partial}{\partial p_{i}}}\wedge{\frac{\partial}{\partial q_{i}}}$, as indicated by the first three terms.
Due to this, PoDiNNs cannot identify their specific coupling pattern.
The springs’ displacements are computed internally within the potential energy function.
However, we can still identify how the dampers are coupled with the masses and the moving boundary, as indicated by the latter three terms.
In one trial, the coefficients of bivector elements between the dampers and masses were disorganized, and this trial corresponded with the failure case shown in Fig. A1.

Identifying Coupling Patterns of System (c) We also examined case of system (c).
We denote the displacement of $i$-th spring in $x$-direction by $q_{xi}$.
The learned bivector $B$ was $B={\frac{\partial}{\partial m_{x1}}}\wedge{\frac{\partial}{\partial q_{x1}}}+{
\frac{\partial}{\partial m_{y1}}}\wedge{\frac{\partial}{\partial q_{y1}}}+{
\frac{\partial}{\partial m_{x2}}}\wedge{\frac{\partial}{\partial q_{x2}}}+{
\frac{\partial}{\partial m_{y2}}}\wedge{\frac{\partial}{\partial q_{y2}}}-{
\frac{\partial}{\partial m_{x1}}}\wedge{\frac{\partial}{\partial q_{x3}}}-{
\frac{\partial}{\partial m_{y1}}}\wedge{\frac{\partial}{\partial q_{y3}}}+{
\frac{\partial}{\partial m_{x2}}}\wedge{\frac{\partial}{\partial q_{x3}}}+{
\frac{\partial}{\partial m_{y2}}}\wedge{\frac{\partial}{\partial q_{y3}}}+{
\frac{\partial}{\partial m_{x2}}}\wedge{\frac{\partial}{\partial q_{x4}}}+{
\frac{\partial}{\partial m_{y2}}}\wedge{\frac{\partial}{\partial q_{y4}}}+{
\frac{\partial}{\partial m_{x1}}}\wedge{\frac{\partial}{\partial q_{x5}}}+{
\frac{\partial}{\partial m_{y1}}}\wedge{\frac{\partial}{\partial q_{y5}}}$ in all 10 trials.
The coefficients were exactly 1.0 or -1.0, accurate to five significant figures.
This is because the velocities of the masses and the displacements of the springs are all observable, and their scales are known.
However, the scales of spring constants $k_{i}$ and masses $m_{i}$ still cancel each other out, which means the parameters cannot be uniquely determined for the system as a whole.

The masses and springs are coupled in either the $x$- or the $y$-directions, which implies that the data is neatly separated along the $x$- and $y$-axes.
From the bivector $B$, we can say the following; mass $m_{1}$ is coupled with springs $k_{1}$ and $k_{5}$, while mass $m_{2}$ is coupled with springs $k_{2}$ and $k_{4}$; and spring $k_{3}$ is coupled with both masses $m_{1}$ and $m_{2}$, but the coupling to $m_{1}$ is in the opposite direction.
This means that the coupling pattern shown in Fig. 2 (c) was fully identified.
By reorganizing the expression for $B$, we get
$B={\frac{\partial}{\partial m_{x1}}}\wedge({\frac{\partial}{\partial q_{x1}}}-
{\frac{\partial}{\partial q_{x3}}}+{\frac{\partial}{\partial q_{x5}}})+{\frac{
\partial}{\partial m_{y1}}}\wedge({\frac{\partial}{\partial q_{y1}}}-{\frac{
\partial}{\partial q_{y3}}}+{\frac{\partial}{\partial q_{y5}}})+{\frac{
\partial}{\partial m_{x2}}}\wedge({\frac{\partial}{\partial q_{x2}}}+{\frac{
\partial}{\partial q_{x3}}}+{\frac{\partial}{\partial q_{x4}}})+{\frac{
\partial}{\partial m_{y2}}}\wedge({\frac{\partial}{\partial q_{y2}}}+{\frac{
\partial}{\partial q_{y3}}}+{\frac{\partial}{\partial q_{y4}}}),$
which indicates that system (c) has eight degrees of freedom.
In this way, PoDiNNs can identify the constraints and handle degenerate dynamics.

Figure: Table A1: Impact of # Components and VPT.

Impact of Number of Hidden Components for System (d) In the electric circuits, the flow for resistors and diodes can be either current or voltage, depending on their coupling with other components.
Let $n_{d}$ denote the assumed number of resistors whose flow is current, and $n_{g}$ the assumed number of resistors whose flow is voltage.
We tested system (d), the FitzHugh-Nagumo model, to explore the impact of the assumed numbers $n_{d}$ and $n_{g}$, and summarized the results in Table A1.
The correct numbers for this model are $n_{d}=1$ and $n_{g}=1$.
When fewer components are assumed, VPT values became significantly poor.
On the other hand, assuming more components yields similar accuracy to the correct configuration.
This trend is consistent in both the training and test subsets.
Therefore, similar to the case with system (b) in Table 4, we can identify the number of unobservable energy-dissipating components and their flow types using the training subset.

### E.2 Additional Discussions

Non-Identifiability The characteristic scales of the energy-storing and dissipating components (e.g., spring constants $k$, masses $m$, and damper constants $d$) and the coefficients of the bivector elements among them cancel each other out.
Hence, the overall scale of the system is indeterminate.
Once the coupling pattern has been identified, the coefficients of the bivector elements can be normalized, as shown in Section 4.2.
If external forces are present, they determine the scale of force, which, in turn, uniquely determines the scale of the coupled mass.
This cascading effect determines the system’s overall scale.

When the characteristics of elements are linear, different coupling patterns can result in the same dynamics.
Additionally, a weighted average of these coupling patterns can also be a valid solution.

While non-identifiability may be problematic for system identification, it does not impact prediction accuracy.

Passivity and Convexity In understanding the behavior of dynamical systems, passivity can be a crucial concept (Khalil, ).
Passivity implies that the system always dissipates energy and, without external inputs, converges to a stable invariant set.
In PoDiNNs, general neural networks are used to model energy-dissipating components, which may sometimes learn energy-supplying characteristics, like a negative resistance.
This flexibility is beneficial for learning special diodes used in the FitzHugh-Nagumo model or Chua’s circuit, but it poses a problem when passivity is desired in the system.
A simple solution to this issue is to replace the general neural networks with ones specialized for functions that have zero output for zero input and are monotonic non-decreasing (Wehenkel & Louppe, ).
This ensures the system remains passive.
If convexity is required in the energy function, a neural network enforcing this property can be introduced (Amos et al., ).

PoDiNNs decompose coupled systems into individual components, along with their coupling patterns, and learn the specific characteristics of each component.
This structure allows for the introduction of constraints beforehand and ensure each component to meet the necessary properties.
This flexibility is another key advantage of PoDiNNs.

## Peer Reviews

confidence: 4
contribution: 3
initial_rating: 5
presentation: 3
questions: The models considered all seem to be linear.  Is that correct?  If so, more standard system ID or linear methods should be considered instead of all this sophisticated ML/AI architectures.  

Can this generalize to nonlinear models?  Or simply be applied to nonlinear models with success?  Although you did Fitzhugh-Nagumo and  Chua's model, these are only "slightly" nonlinear and system ID models can work pretty well with those.

What have you not considered non-neural network methods like DMD... or system ID methods like mamba/S4 not been compared?
soundness: 3
strengths: The integration of physics-based principles is always welcome in the dynamical systems world.  The authors have a very nice contribution to make here potentially as the integration of physics principles into neural networks is very important.
summary: The authors present a new architecture, specifically they integrate some clever ideas of Dirac structure in order to unify the port-Hamiltonian and Poisson formulations from geometric mechanics.  Both nice innovations for physics-constrained neural networks.
weaknesses: The models used seem to be all linear: which begs the question about simple linear model regressions such as dynamic mode decomposition, dynamic mode decomposition with control and time-delay embedded DMD which models missing/coupled physics.  These more baseline (non-NN) methods are simply not talked about or considered and I think they should be.

There are statements that are simply not true:  "two key limitations remain in modeling dynamical systems, especially those described by ordinary differential equations (ODEs). The first limitation is the narrow focus on mechanical systems."  This suggests the authors don't know the field well and it is a concern.  People are modeling all kinds of dynamical systems with ML/AI architectures well beyond mechanical systems.

Further: "The second limitation is that most methods treat the system as a single, monolithic entity."  This is also not true.  Many people are working on coupled systems where time-delay embeddings are often used to gather information for missing, coupled and unmeasured variables.  Again, it is odd that these statements exist in the paper which suggests the authors are not aware of the great body of work on model discovery and dynamical systems methods with ML/AI for systems which are not mechanical and which indeed have coupling.

---

confidence: 4
contribution: 3
initial_rating: 5
presentation: 4
questions: Please refer to the bullet points in Weaknesses.
soundness: 3
strengths: •	As the reviewer summarizes above, key limitations are correctly identified such that the contributions of this work are clear. To the best of the reviewer’s knowledge, this work is new.

•	The theoretical analysis is solid, where the definitions and theorems clearly show how the Dirac structure encapsulates internal and external component couplings. Moreover, the corresponding examples of different dynamical systems are well-explained to differentiate the proposed method from existing methods like HNN and NODE.

•	The results on multiple systems look promising, where various experiment scenarios and evaluation metrics are comprehensive to validate PoDiNNs’ capabilities.
summary: This work targets the dynamical system identification using observation data, which is a hot topic and essential application. The key differences with respect to existing works are clearly stated: 1) identifying unknown physics rather than predefined symbolics, 2) considering the coupling behaviors in the system or interactions when the problem covers multiple domains rather than single mechanical systems.
weaknesses: Some technical details, as well as the claimed capabilities, are unclear, which might be because the reviewer is unfamiliar with all kinds of multi-domain dynamical systems. The confusions are listed below.

•	The capability to deal with multi-physics problems is claimed several times in the paper. Specifically, Remark 1 explains the representation of inter-coupling using a bivector element (which is also an NN, right?) Remark 4 with Table 2 briefly demonstrates the scenarios to capture coupled physics in multiple domains. The reviewer would like to know how PoDiNNs represent such interactions.  Is it the same way as the traditional simulation tool, e.g., through iterative refinement of two (or more) simulations of a single domain or system? 

•	The proposed work targets unknown physics/dynamics, which is quite challenging as there are no predefined physical symbolics in PINN-alike works. How to ensure the PoDiNNs capture the correct physics without causing overfitting problem or continuous good performance in extrapolation?

•	Moreover, PoDiNNs focus on behaviors that may not be captured by generic models, e.g., NODE, which seem to need intensive resources. Especially, the couplings in multi-physics usually require heavy computation in traditional simulations. The more fine-grained, the heavier. What is PoDiNNs’ capability in this aspect?

•	For the last paragraph of Sec. 3.4, an example of electric circuit is used. Could the reviewer further explain with more details: why ODEs or using NODE alone cannot capture the current flow and balanced voltage level? Subsequently, how does PoDiNN mitigate the issue of limited representation? A toy example with mathematical derivations or diagram will be helpful.

---

confidence: 1
contribution: 3
initial_rating: 6
presentation: 3
questions: -
soundness: 3
strengths: **S1:** The empirical results show that PoDiNNs achieve lower errors and provide stable predictions for longer than Neural ODEs and HNN variants. As such, the method seems effective for the tasks for which it was designed.
summary: This work studies the application of deep neural networks for data-driven physics simulation. Over the recent years, several neural network simulation models have been proposed among which Hamiltonian Neural Networks (HNNs) and Lagrangian Neural Networks. Although these models are formulated along physics principles and consequently enjoy improved stability of their simulations, two key limitations remain: (1) they focus on mechanical systems, as opposed to for example electric circuits or magnetic fields, and (2) they consider the system as a single monolithic entity. To alleviate these limitations, this work proposes Poisson-Dirac Neural Networks (PoDiNNs), which unifies the port-Hamiltonian and Poisson representations, and explicitly represents the coupling between internal and external components. The empirical evaluation shows that PoDiNNs enjoy stability over longer simulation horizons and in general achieve lower errors than Neural ODE or variants of HNNs.

Unfortunately, I am missing too much of the mathematics and physics background that is required to understand the paper, so I am unable to provide an informative review.
weaknesses: **W1:** The content might be difficult to digest for the audience of ICLR, since a specific mathematical and physics background is required to understand the work. In my impression, most ML researchers will lack this background, and also a substantial part of the AI4Science community might find it difficult to understand the paper.

---

confidence: 2
contribution: 4
initial_rating: 6
presentation: 2
questions: See above.
soundness: 3
strengths: - The paper is well-written and concise. The modelling choices made by the authors are well-motivated; the incorporation of Dirac structure into a DL-based solving method for mechanics has a lot of possible benefits regarding guarantees and interpretability of the learned functions.
- The experimental results are very convincing; in the chosen experimental setups the proposed framework outperforms all baseline methods consistently.
- The experiments on identifying coupling patterns are fascinating, these results show that the proposed framework is indeed able to recover the underlying coupling in a very interpretable manner.
summary: This paper introduces a new method for stable deep-learning based modelling of dynamical systems. Neural ODEs were originally found to have significant limitations in modelling dynamics long-term due to a lack of physical guarantees (e.g. conservation laws) leading to error accumulation. To address this, Hamiltonian NNs incorporate Hamiltonian mechanics (modelling the hamiltonian through a NN). This paper contributes a generalization of Hamiltonian NNs; Poinsson-Dirac Neural Networks (PoDiNNs) to combine both port-hamiltonian and poisson systems for mechanics, enabling learning of coupled dynamical systems across multiple domains, their interactions and degeneracies. This  constitutes a unification of previous works that address/leverage specific properties of dynamical systems like energy dissipation and external inputs. Authors show that this approach moreover allows for the identification of the specific coupling patterns between interacting dynamical systems. The experimental results show clear and consistent improvements over previous approaches in a range of systems with degenerate dynamics (e.g. constraints), dissipation and external inputs.
weaknesses: - The paper is quite dense, making some of the reasoning and motivation hard to follow. Although in places authors elaborate on their arguments with examples, I think it would be good to provide some more visual guidance. For example, figure 1 contains a very high-level overview of the concept of PoDiNNs, but I think it would be good to expand this diagram and clearly indicate what parts of this system the method proposes to replace with learned functions (NNs) and relate it to e.g. eq 5, def 3. I think this would greatly help with readability.
- The experiments chosen by the authors are examples of systems that fit the modelling criteria set by the authors for their framework. In these settings, the proposed framework outperforms previous approaches. However, I think it would also be good to consider/experiment with settings that do not exhibit e.g. degeneracy just to gauge how well this generalized method works in the settings that the originally proposed HNN and LNN are validated on. Would it for example be possible to show performance of your model vs baselines on eg double pendulum (as in HNN or LNN paper)?
- For completeness I think it is important to compare your method against the baselines also in terms of computational complexity / overhead. Currently, no details on the computational complexity or time complexity on either training or inference of the proposed method are given. Please provide these details.
- An arguable weakness of the method is the reliance on specifying problem-specific constraints, i.e. the number of components is a model hyperparameter, but as shown in the last experiment / Appx E it is important to specify this to at least the number of underlying components in the system, or performance will be significantly impacted. Overspecifying is possible of course, and does not seem to have negative impact, but induces additional computational complexity.

---

confidence: 2
contribution: 3
initial_rating: 6
presentation: 2
questions: 1) While the PDNN shows promising results, it would be interesting to compare its performance with state-space model (SSM) based neural operators on the datasets used in this paper.

2) To make the paper more accessible to a broader audience, the authors could consider adding an appendix to provide additional explanations and visualizations of the key concepts and techniques.
soundness: 3
strengths: 1) The paper provides a solid theoretical foundation for the proposed PDNN architecture.

2) The authors conduct extensive experiments to demonstrate the effectiveness of the PDNN in various applications.
summary: This paper introduces a novel neural network architecture, the Poisson-Dirac Neural Network (PDNN), which combines Hamiltonian and Poisson formulations. The PDNN leverages the Dirac structure to model coupled dynamical systems across various domains.
weaknesses: 1) While the paper's theoretical contributions are significant, additional background material, such as differential geometry, etc, in the supplementary information would enhance its accessibility to a broader audience.

Meta-review: The reviews for this submission were leaning towards acceptance of the paper. Among other things, reviewers highlighted as strengths the use of well-explained examples of different dynamical systems that allow to differentiate the proposed method from existing methods, and which demonstrated convincing results in terms of lower errors and more stable predictions over longer time horizons. The presentation of the paper received both praise and requests for improvement, with the accessibility of the paper being the most important improvement area. Several reviewers raised questions about the computational cost of the proposed method compared to baselines, which the authors answered transparently during the rebuttal. The reviewer with the most critical review requested comparisons against other methods such as system ID methods, but has unfortunately not engaged with the authors anymore after their rebuttal, despite requests to do so. Despite one other reviewer agreeing that the paper would be strengthened by a comparison to these types of methods, this reviewer maintained their recommendation that the paper is marginally above the acceptance threshold. Taking all of this into account, the recommendation is to accept this paper. The authors are strongly encouraged to include the additional clarifications that were helpful during the rebuttal in their camera ready version of the manuscript.

Decision: Accept (Poster)