Poisson-Dirac Neural Networks for Modeling Coupled Dynamical Systems across Domains
Authors: Razmik Arman Khosrovian
Takaharu Yaguchi
Hiroaki Yoshimura
Takashi Matsubara
Venue: ICLR 2025 Poster (2025)
Keywords: neural ordinary differential equations
coupled system
Poisson system
Dirac structure
Primary area: learning on time series and dynamical systems
TL;DR: Poisson-Dirac formulation with ports enables neural networks to model various dynamical systems across domains and identify their internal structures.

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.

Published as a conference paper at ICLR 2025
POISSON-DIRAC NEURAL NETWORKS FOR MODELING COUPLED DYNAMICAL SYSTEMS ACROSS DOMAINS
Razmik Arman Khosrovian Takaharu Yaguchi Hiroaki Yoshimura² Takashi Matsubara" OOsaka University "Kobe University "aseed University "Hokkaido University
ABSTRACT
Deep learning has achieved great success in modeling dynamical systems, provid- ing 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 tendeney treat systems mono- lithic. These limitations reduce their applicability to dynamical systems in other domains, such electrical and hydraulic systems, and to coupled systems. To ad- dress these limitations, we propose Poisson-Diis -Drac Neural Networks (PoDiNNs). novel framework based on the Dirae structure that unifies the port-Hamilnonian 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 experi- ments demonstrate that PoDiNNs offer improved accuracy and interpretability modeling unknown coupled dynamical systems from data.
INTRODUCTION
Deep learning has achieved great success in modeling dynamical systems (Chen al. 2018 Anand-
2020). following its successes in image processing and natural language processing (He Vaswani et al., 2017). It provides data-driven approach for predicting the behavior of complex dynamical systems, even when governing equations are unknown (Kasim & Lim 2022 Matsubara & Yaguchi 2023 These models function as computational simulators and show promise in applications acToss diverse fields, including weather forecasting, mechanical design, and system
kumar et al.
al. 2016
control (Lam et al. 2023 Pfaff al., 2020; Horie al. 202I).
However, Greydanus et al. 2019) identified a key limitation of data-driven models: they accumulate modeling errors long-term predictions, leading to rapid failure. Hamiltonian Neural Networks (HNNs) were proposed to overcome this limitation by incorporating Hamiltonian mechanics. In- spired by this approach, many studies have developed neural network models that not only learn the superficial dynamics from data, but also adhere to fundamental laws of physics. Examples include Lagrangian Neural Networks (LNNs) (Cranmer et al.) 2020). Neural Symplectic Forms (NSFs) (Chen al. 2021). Poisson Neural Networks (PNNs) (Jin et al. 2022). Constrained HNNs (CHNNs) (Finzi et al. 2020), and Dissipative SymODENs (Zhong et al., 2020), as shown in Table
Despite recent progress, two key limitations remain in modeling dynamical systems, especially those deseribed by ordinary differential equations (ODEs). The first limitation is the narrow focus on me- chanical systems. Models applied other domains, such as electric circuits magnetic fields, often fail to leverage the governing principles of those systems, such as Kirchhoff's current and voltage laws (Jin al. 2022; Matsubara & Yaguchi 2023). 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 & Mars- den, 2006a). Although the port-Hamiltonian formulation theoretically addresses these interactions, no prior work has fully leveraged its potential. These limitations prevent current methods from ef- fectively handling multiphysics scenarios, where systems from different domains interact, such as those involving DC motors (yan der Schaft & Jeltsema 2014 Gay-Balmas & Yoshimura| 2023).
"Correspondence u1111282 13287ddecca. osaka-u.ac jp and matsubare ra st. hokudaii acc. jp.

---

Published as a conference paper at ICLR 2025
Table 1: Comparison of Methods for Modeling Dynamical Systems
Identiyinng Coupling Patterns
External Iopats
Coordinate Free
Applicable to
Multiphrdes
Degencracy
Disspation
Model
Formulation
HNN Greydanus cl.
2019
canonical Hamiltonian
LNN
2020
Lagrangian
NSF Chea et al. |222 2021
general Hamiltonian
X
X 2
2020
constrained canonical Hamiltonian
CLNN
1azz m 2040
constrained Lagrangian
X
X
Dis. SyMODEN Zhoing al.
2020
portHamiltonian Darboux coordinates
X
PoDiNN (propssed)
Poisson- Dime with ports
*Available only for known holenomic constraints.
f Available only for externa force on mass.
Coupled system
To address these limitations, we propose Poisson-Diie Dirsc Neural Networks (PoDiNNs). which leverage the Dirac structure to unify the port-Hamiltniaan and Poisson formula- tions (Courant 1990; Duindam et al., 2009 van der Schaft 1998; van der Schaft & Jeltsema
energy-storing components 80 -00
external Inputs
Dirac structure (fec) D2
entrgy-storing components
energy-dissipatinn components
flow
fo
effort
(014imosu & Marsden 2006a). The Dirac structure explicitly represents the coupling be- tween internal and external components. A con- ceptual diagram is shown in Fig. [ The advantages of PoDiNNs are summarized below and in Table T and validated through experiments on seven simulation datasets spanning mechanical, rotational, dlectro-magnett and hydraulic domains.
Figure Conceptual diagram of PoDiNNs.
Many real-world systems are coupled systems, consisting of in- teracting components (Yoshimura & Marsden, 2006a). 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.
Identifying Coupling Patterns
Applicable to Various Domains of Systems The Dirac structure allows PoDiNNs to describe various dynamical systems across multiple domains. While most neural network models of ODEs focus mechanical systems, PoDiNNs 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 pects of dynamical systems, such as degenerate dynamics (Finzi et al,, 2020 Jin et al. (2022), energy dissipation, external inputs (Zhong al., 2020). and coordinate transformations (Chen al. 2021 Jin et al.] (022). PoDiNNs offer a unified framework to address all these aspects together by lever- aging the Dirac structure.
BACKGROUND THEORY AND RELATED WORK
Attempts to learn ODEs have a long history Nelles 2001 I see also Appendix A for broader review. Neural Ordinary Differential Equations (Neural ODEs) transformed this field (Chen al. 2018), which train networks to approximate the vector field (the right-hand side of the ODE), which then integrated numerically. However, they are often overly general and fail capture dynamical systems governed by principles such as energy conservation (Greydanus et al.] 2019). 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 pair of n-dimensional generalized coordinates q (q).-.., - Qn) and
2

---

Published as a conference paper at ICLR 2025
momenta p (P1...., PP)) as the state. Given energy function H, the system evolves as:
q
0H 0p:
p. =
aH or equivalently, --1 āy."
(I)
where denotess the time derivative, I and O are the n-dimensional identity and zero matrices, respectively, and VH denotes the gradient of H. Systems describable by Hamilton's 's equations are Hamiltonian systems. See Appendix for an example.
HNNs were introduced to learn such systems, with a neural network HNy serving as the energy function H (Greydanus et al., 2019). HNNs have shown greater long-term robustness than Neural ODEs (Chen et al. 2022; Gruver et al. 2022).
Hamiltonian Systems of Coordinate- Form From the differential geometry perspective, the equations of motion in Eq. (1 are defined on the cotangent bundle M = Q of the configuration space Q R with the standard symplectic form Ω = - dg; A dpt. where q E Q and p E Ty Q, respectively. Here, Ty Q denotes the cotangent space to Q at q, and A denotes the exterior product. The symplectic form 52 leads to a skew-symmetrie linear bundle map St Tu.M + T at each point u .M, which satisfies (S2 (v), w) (u, w) for any u E M and v, E Tu.M, where (+-)) denotes the natural pairing. We will denote by the subscript an assignment at point u and omit the subscripts if an equation holds independently of the point u. Given smooth function H M + IR, its differential is denoted by dH. A vector field X on M assigns a tangent vector X_ E M to each point M. Then, Hamilton's 's equations in the coordinate-free form defines a vector field X called the Hamiltonian vector field by
(2(HH) dH.
Xи defines the time evolution of the state u as ú(() = (Xe) at time t, and u(t) over a certain period is solution of the system. Then, the tuple (H, S2 X²) called a Hamiltonian system. A Hamiltonian system with the standard symplectic form f2 is said to be canonical, and its coordinate system is known as Darboux coordinates.
Hamilton's 's equations also describe the same dynamics using generalized velocities Ta Q instead of generalized momenta p E Ty2. 2, i.c., on the tangent bundle T'Q rather than the cotangent bundle Q. In this case, the symplectic form 12 is replaced by Lagrangian 2-form (Marsden & Ratiu 1999). From another viewpoint, the symplectic form Ω defines the coordinate system. Lagrangian Neural Networks (LNNs) learn the dynamics Neural Symplectic Forms (NSFs) learn the symplectic form directly from data, generalizing HNNS and LNNs for arbitrary coordinate systems (Chen al. 202T). In general, Hamiltonian systems preserve the energy H, as Cx H- = (dH, (dH,HH) = 22(X H, Xx) ) 0, where Lx,, denotes the Lie derivative along XH, and the last equality follows from the kymmetyye of 92.
the tangent bundle TQ (Cranmer al., 2020)
Poisson Systems Because the symplectic form Q is non-degenerate in the sense that the bun- dle map Dh, is non-degenerate, it leads to a 2-tensor B called a Poisson bivector satisfying B(dH, dG) Q(X n, Xa) for any smooth functions H, G on M. B leads linear bundle map Bi T
skew-symmetric M + TuM at each point u E M, which satisfies Bu 3) = (a, 1t, (B)) for any u E MM and a E M.M. Then, it holds that B = (12-)1 Hamilton's equations, Eq. (2), are rewritten as
Xи = (dHH).
The tuple (H, XR) is called a Poisson system. On the Darboux coordinates, B=
A Ju.
Lie-Poisson neural networks were proposed to learn this system with known B EEdredd et al. 2024
2.2 DEGENERACY, DISSIPATION, AND EXTERNAL INPUTS
Degenerate Systems the state u is constrained to submanifold M c M, Hamilton's equations do not directly describe the dynamics. For example, consider a pair of mass-spring systems, indexed by i E (1,2), with a constraint @i = G2 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 's equations (²)) = dH with the standard symplectic form Ω do not describe the dynamics. There are two primary methods for handling such degenerate Hamiltonian systems.
3

---

Published as a conference paper at ICLR 2025
The first approach uses coordinate transformations to reduce the system's 's degrees of freedom and define a submanifold M M, where Hamilton's equations describe the dynamics using the stan- dard symplectic form or Poisson bivector B. The Darboux-Lig theorem ensures the local existence of such transformations. The standard Poisson bivector on M can also be expressed on M while its bundle map is degenerate, and it represents the coordinate transformation. Hence, degenerate Hamiltonian systems are included into Poisson systems. See Appendix [B.2] for an example. [i al.] (2022) introduced Poisson neural networks (PNNs) that combine neural network-basd -hased coor- dinate transformations with SympNets for degenerate systems (Dinh et al., (201)] [Jin et al.) 22200) However, this approach lacks interpretability due to non-unique, nonlinear transformations, and ex- tending 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. (2020) 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-Hamiltonial Systems To incorporate dissipation and external inputs, some studies have ex- plored the port-Hamil Hamilttnonnn formulation with neural networks (Zhong et al. (0020). Although some have mentioned the underlying Dirac structure (Neary & Topeu| |
[Di Persio et al. 2024), they
rely on the following canonical form without fully leveraging its flexibility:
- - (8 Dia))) VH(q. + (c(g)..
(4) where D E Rn*x represents dissipation, f E IR"' is the control input vector, and G' E R= is the control input matrix. This formulation has several disadvantages. The matrix ( () represents the standard symplectic form 52 = Σ, dg: Adp; on R2n, which cannot handle degeneracies. Dissipation from multiple sources, such as dampers or friction, is condensed into single term D, limiting interpretability and modeling performance. Furthermore, this formulation cannot handle externally defined velocities, as seen in models of buildings shaken by the ground. While [idnes et al. (2023) employed non-standard form to address the first limitation, their method still inherits the latter two.
POISSON-DIRAC NEURAL NETWORKS
To overcome the limitations of existing methods, we propose Poisson- Dirae Neural Networks (PoDiNNs), which leverage the Dirac structure to unify Poisson and port-Hamiloniia systems. Detailed background theory and proofs of the following theorems are provided in Appendix
3.1 DIRAC STRUCTURE
Let V be an re-dimensional vector space and V* its dual space, with the natural pairing (+-) between them. Define the symmetric pairing (- >) on V V*
((v,,)) a), (8, @))) (ar, F) + (ā, v) for (v, (8, E VVVV
where ( denotes the direct sum of two vector spaces (Courant 1990).
Definition (Courant (1990); Yoshimura & Marsden (006ab). A Dirac structure on a vector space
V is a vector subspace D C V V- such that D
D'. where D is the orthogonal complement
of D with respect to the pairing ((, .).
(ū, a)) E D, and hence (a, v) = 0 for any (v. a)) E D. A typical Dirac structure can be constructed as follows:
Since D= D-, (ar, 6) + (ā, v) 0 for any (v, (ū,
Theorem 1 (Courant (190)); Yoshimura & Marsden 2006a)). Let V be vector space and A
vector subspace. Define the annihilator A of A as A {a E V* | (a, v) = for all v E A} C
V* Then, Dy A @ A C @ V* is Dirac structure V.
A vector bundle F over a manifold M is a collection of vector spaces Fu, called fibers, smoothly assigned to points u E M, where the manifold M is called the base space. A typical example is the tangent bundle T.M, where the fiber at u is the tangent space TuM, and its dual bundle is the cotangent bundle M. The Whimey sum O defines 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 FF over M and its dual & = Fr, the Dirac structure is defined as a subbundle of J E.
4

---

Published as a conference paper at ICLR 2025
Table 2: Categorization of Components in Different Domains
Domain
Mechanical
Rotational
Eleetro-Magnetic
Hydranlic
Subdomain
Potential
Kinetic
Potential
Kinetic
Electric
Magnetic
Potential
Row (impul)
velocity
force
angular velocity
lorque
current
volluge
volume flow rale
effort (oulput)
force
velocity
torque
angular velocity
vollage
current
pressure
displasement
momentum
angle
angular momentum
electric charge
magnetic flux
volume
energy energy-storing
spring
muss
(potential) friction
inertia
capacitor
inductor
hydmulic
energy-dissipating
damper
resistor
resistor
external imput
external force moving boundary
external torque
vollage source
current source
incoming fluid flow
Definition 2 (Courant 1990); Yoshimura & Marsden (2006a); van der Schaft & Jeltsemal (014
Consider a vector bundle Ff over a manifold M. A distribution A is collection of vector subspaces Au C F each assigned smoothly to M at point u, forming a vector subbundle of F. The annihi- lator A of A is collection of the annihilators Au of Au. also forming subbundle of & = FF. Then, a Dirac structure D is constructed as D A + A which a subbundle of E.
If F + 3 TM # M, the Dirac structure 1) can reformulate Hamiltonian, Poisson, and con- strained Hamiltonian systems (see Appendix DD. Here, we assume that the fibers Fu and E of the vector bundles F and & at u are decomposed as
Fu= OF O and = e c E
(5)
where FS T and 33 = M. A point on the fibers Fu and Ey is denoted by f ( and = (eS eB ee),, respectively. We refer to f E Fu as flows, E & as ef- forts, and both collectively as port variables. Note that our definitions mainly followed those in Duindam et al. (2009); van der Schaft Jeltsemal (2014), while we can find other definitions in Yoshimura & Marsden 2006a).
Theorem 2. Consider vector bundles F and & = F² over M. The collection of
Du e( x Eu f (e)
for the bundle map Bt, Eu + F of bivector B is a Dirac structure F
Here, we define PoDiNNs as special case of the Poisson-Diral -Dircc formulation (Courant 1990
Definition 3 (Poisson-Dirac Neural Network). Let F @ be a vector bundle over a manifold M defined in Eq. 5 which assigns to each M a vector space (F a @F FF.) (ES a &f Ell a &d). Let D C F @ & be a Dirac structure defined in Theorem Let H : :MM- + R be an energy function, ER determined by a mapping Ru FR + E. of the flow f E FR. The effort e' (t) E E., is a time deeeddeendentt function. The functions and Ru, are implemented using neural networks. for each u(t) and E [a, b), it holds that
which determines the effort e dHl E E The effort el
(1F* (()), (') (1), (e (1), en() el),e'(() (1)) E Duct)s the tuple (H, B, R, e
,J ) is called Poisson-Diee Neural Networks (PoDiNNs).
3.2 FLOWS AND EFFORTS FOR COMPONENTS
The point u at the base space M represents the states of the dynamics, which include the displace- ment of a spring, the momentum of a mass, the angle and angular momentum of a rotating rod, the electric charge of capacitor. Intuitively, flows are the inputs to components, while efforts e are the outputs from the components. Components considered in our formulation are summarized in Table 2 with concrete examples in Appendices C and [F]
Similar to Hamiltonian and Poisson systems, the flow fs E Tu.M represents the vector field X on M, defining the time evolution ú of the state u. The effort e E T corresponds to the differential d/f of the Hamiltonian H .M- + R. The superscript $ indicates that these components store energy. In electric circuits, capacitors and inductors are examples, with states as electric charge and magnetic flux, and efforts 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.
5

---

Published as a conference paper at ICLR 2025
The flow f E F& and effort e E E. represent pvergy-dissppating components, such dampers and resistors. A damper's flow f is the velocity (i.e., the rate of extension or compression), and its effort e is the force, which are related as e -df Rt for linear damper with damping coefficient d. Unlike a spring, a damper does not store its own energy, so its flow f is not on Fu but on F&R The superscript R indicates that these components are "resistive." However, components that supply energy, such as special diodes, can also be classified into this category. In any case, each component is implemented PoDiNNs by neural network that approximates its characteristic Ru e
The flow f' E Fll and effort e' E represent external inputs, such as an external force (where force is the effort) or a moving boundary (where velocity is the effort). The superscript indicates that they are "inputs." Their efforts e depend only on time I, not on other components. Their flows f are not required for determining the system's '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
2 For the coordinates di and pæ on M, the basis vectors of the tangent space Tu,M are Dy. and op. respectively. The i-th basis vectors of FRl Fl. Ed and El are denoted by Et E/ <.* and e!- respectively. The bivector assigns to each point u E .M wedge products of the basis vectors of the flow space Fu, such as 2 0 9 ap. 0gs", 0p, AEE, and 0 ^EJ, thereby defining the coupling patterns 0p.
8 among components. For example, Sp. A Ej couples the j-th external input with the i-th mass m; c with the state p.. The effort the j-th external input is expressed as e/5 when the basis is explicit. a_ This is fed to the bivector öp, A = A op. resulting in -cj Up. which forms part of the flow I 0 for mass m., Therefore, we can make the following remarks. öp.
a
Remark 1 (Coupling as Non-zero Elements of Bivector). Coupling between two components is represented by 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 's components.
Remark 2 (Degeneracy of Dynamics Degeneracy of Bundle Map). Constraints between compo- nents that cause degenerate dynamics are reflected in degeneracy of the bundle map Bt. 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), formulatiioon can represent systems across various domains, as sum- marized in 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 C for concrete examples.
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 unified manner.
Previous models, including HNNs, LNNs NSFs, PNNs, and Dissipative SymODENs, approximate the Hamiltonian H, Lagrangian L, 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 from the energy functions H, improving interpretability and gen- eralization performance.

---

Published as a conference paper at ICLR 2025
000
0.9
R
m
m2
m
(f) DC motor
with external force
4-
4
(d) FitzHugh-Nan Nagumo model
10
17
m2
m2
R.
Kg
R
(b) with moving boundary (c) with redundancy
(a)(c) Mass-spring(-dammer) systems
A.P
d2
(e) Chua's circuit
(g)
Hydraulic tank
Figure 2: Diagrams of systems that provide datasets. Detailed definitions are found in Appendix
F
In an absolute coordinate system, the displacements of springs were computed from the absolute po- sitions of springs' edges internally within the energy function. Although PoDiNNs can handle such system, the bivector for energy gergyy-torinn components takes the standard form B = Σ, öp, which hinders the identification of their coupling patterns. PoDiNNs are still able to identify cou- pling patterns involving energy-disse mneggyddssiipating components and external inputs. To identify the coupling patterns between energy-storing components, it is necessary to employ relative coordinate system based on displacements, rather than absolute positions.
a A og;
0
LNNs and NSFs are designed for non-canonical Hamiltonian systems. where the Hamiltonian vector field Xи is implicitly defined as ( dH (Cranmer et al. 2020, Chen et al. 2021). These require matrix inversion to compute XHr, which is computationally expensive, numerically unstable, and are not directly applicable to degenerate dynamics. In contrast, PoDiNNs explicitly define X as part of f (e), offering faster and more stable computations and handling degeneracy.
In mechanical systems, energy-dissipating components such as dampers and friction are character- ized 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 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 electric circuits are generally described by differential-algebriie equations (DAEs). For example, when capacitor is connected 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 wide range of systems and expand the scope of modeling unknown dynamical systems from observations.
direct voltage source in parallel, infinite current
See also Appendix for implementation and further discussions.
4 EXPERIMENTS AND RESULTS
4.1 EXPERIMENTAL SETTINGS
Datasets We evaluated PoDiNNs and related methods to demonstrate their modeling performance using seven simulation datasets, shown Fig. 2 Due page limitations, we briefly overview their characteristics. The full explanations can be found in Appendix
For the mechanical domain, we prepared three mass-spring(-dantee locities of the masses as observations. In the absolute coordinate system, we used the absolute positions of the springs' ends as the springs' states, and the relative coordinate system, their dis- placements. As external inputs, an external force F is applied to system (a), and 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, but PNNs can. In all three systems, the springs and dampers exhibit nonlinear characteristics.
systems (a)-(c), with the ve-
7

---

Published as a conference paper at ICLR 2025
We selected two nonlinear electric circuits, (d) the FitzHugh-Nagumo Nagumo model and (e) Chua's cir- cuit, from the electro-magee megetro-magnetic domain (Izhikevich & FitzHugh 2006 Chua 2007). The capacitor voltage and inductor current were used as the observations.
We also used two multiphysies systems (f) and (g), system (f), a DC motor bridges an electric circuit in the electro- magnettc domain and a pendulum in the rotational domain. In system (g), 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 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 (Paskke et al., 2017). See also Appendix F 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) in the absolute coordinate system, while PNNs were evaluated on system (c). 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 = where m is the mass and p the momentum. Therefore, we employed this form with learnable parameter m, rather than 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ı..... 9n plus the position q of the moving boundary collectively by 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 snergy-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 random state snapshot (*) at n-th step, each model predicted the next state (m++)) after time step AV. Then, all parameters were updated to minimize the squared error between the predicted state a (m++1) and the ground truth u (1+ normalized by state standard deviations. It is known that longer prediction steps can improve robustness against noise (Chen et al. 2020). However,
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 (n) and the ground truth (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 the overall MSE;
MSE u) = / E [MSE(U(") ((("))]
(6)
Lower values indicate better performance. Additionally, we defined the valid prediction time (VPT) the ratio of the number of steps taken before the MSE first exceeds a certain threshold 0 to the total length N of the test trajectory (Botev et al. 2021 Jin et al. 2020; Vlachas et al. 2020);
VPT(ū; u) Nv arg maxn, (n/ (SE(((") (n) 0 for all n
Higher values indicate better performance. The threshold 0 was set 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.
4.2 RESULTS
Numerical Performance The results in Table 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 model ignores oscillations in the
8

---

Published as a conference paper at ICLR 2025
Table 3: Experimental Results.
Mprinn Saans--pinn-D Systems
Dataset
with external force
with moving boundary
(c) with redundancy
Coordinate
(a) relative
(a') absolute
(b) relative
(b") absolute
relative
Model
MSEJ
VPT1
MSEJ
VPTT
MSEL
VPT
MSE.
VPTT
MSEJ
VPT?
Neural
0.92 0.128-0.021
1.68 0.097=0.008
1. 40009
56 0.035
1190.61 0.00
HNN Variant
8.31 20.56
0. 00.007
5.92 12 0.001 40000
200.22
0.000
PoDiNN
4.33-0026 0.622=0.00 7. =0,37
0.26 0.856 8.
0,041
0.11
0.867
x10
e 10
x10
10
X30
6 10
x30
6 = 10
X1U
6 30
Ellcttic Circuits
Multiphysics
(d)
FitzHugh-Nagumoo
(e) Chua's
(D) DC Motor
(g) Hydraulic Tank
Model
MSEL
VPT!
MSEJ
VPTI
MSE
VPT1
MSEL
VPT!
Neural ODE 48.96
11.43 00022
14.744133 0.287 10.06 16.03 #5.15 0.276-0.10
30.62 28 0.045
PoDiNN
1.64
20.022
9.21
0.46910.010
:2,90
0,923 2001
5.42 0,918:00013
X30
8= 10
x10
== 10
x10
= 10
x
0=
Each score represents the median over 10 trials, followed by the d. symbol and the quartile deviation. and (b), and PNNs for system
DDsssipative SymODENs for systems
2.5
ground truth
0.5-
MMMAAM A
couid trach
PAVLROI MAA MAN
0.0
-.55
0.3-
2.5
wMY
nearol ODE
neural DDE
0.00
0.00
5-
neural ODE
0.05-
neural ODE
0.05
cis. SymODEN
0.5-
PNNN
0.00
8.00
PoDINN
PaDiNN
0.00
PoDiNN
0.005
PcDiNN
0.000
20
d0
20
10
20
25
50
75
100
(a') with external force in absolute coordinate system
(c) with redundancy
(e) Chua's 's circuit
(g) Hydraulic tank
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 cemi-transparent color. See also Fig. Al
trajectory and outputs an average path, which keeps the MSE low but results in poor VPT. This
suggests that VPT is more reliable metric than MSE (Botev et al. 2021 Jin et al. 2020 et al. 2020).
Vlachas
We show example trajectories and absolute errors in Figs. and A1
PoDiNNs perform well both in the absolute and relative coordinate systems, showing their adaptabil- ity 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 in- puts using box functions D and
PoDiNNs demonstrate excellent accuracy in system (c) because they successfully simplify the dy- namies by identifying the degeneracy from high-dimensional observations. contrast, PNNs strug- gled to learn the appropriate coordinate transformation. Although the Darboux-Lie theorem guar- antees 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 pattemss in system (e), Chua's circuit. When two coefficients dif- fered by a factor of 1000 or more, we considered the larger one as a detected coupling and the coupling. In all 10 trials, we obtained the bivector à cf ^ σQU
smaller one as effectively zero, indicating B + +225 A a_ with trial-wise positive parameters a, and Eh A 8Q2' were effectively zero. Because R₁ and R3 cannot be directly observed and their indices interchangeable, they were appropriately reordered for analysis.
a A è B = -a 2e
+ lgft ^ aQu C1. and c2 We emphasize that the coefficients of other bivector elements, such as ø5 σश. a
aQz
aQi
9

---

Published as a conference paper at ICLR 2025
We normalized the coefficients of the bivector elements and the characteristics of system com- ponents, as their scales cancel each other out. Then, from the learned bivector B, we can de- rive Kirchhoff's current laws, Lc, -In, + Ix and Ic= = -I,. + In, and Kirchhoff's voltage laws, V, = Vc2 VR; = Vc. - Vcg. and VR, = -Vc,- This allows us to con- struct circuit diagram, which perfectly matches that of Chua's circuit, shown in Fig. 2 (e).
While coupling patterns in more complex circuits may not be unique due to Norton's 's and Thevenin's 's theorems, our formulation can learn one valid realization. Figure Jillustrates the input-output relationship 4
of the neural network representing R2. The results from all 10 trials, shown as colored dashed lines, almost perfectly overlap the true re- lationship, shown as a black dashed line, within the range of +.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.
-
0
3
R2
Figure 4: The identified
characteristies of R2.
Identifying Coupling Patterns of System (g) We also examined the bivector B identified in sys- tem (g), hydraulic tank. In all 10 trials, we consistently obtained the bivector B = kp A 0 av
components, with the coeffi- cients 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 d1 and (22 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 saaplings-betweeee the dampers e and masses, Api A a and a A Bdy and between the external force and the mass, Bpp "E 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.
0.3( öp2 av + Bpp Ban Bpg
a + 0 A à between the
8 ^
adj Bya
Impact of Number of Hidden Components (b) The states of energy-storing components are provided observations, and exter- nal inputs are typically known (or inferred using methods like neural CDE (Kidger al. 2020)). PoDiNNs also require specifying the number and type of energy-dissipaiin ally unknown. To assess the impact of the assumed number of com- ponents, 7dd. we tested system (b) in the relative coordinate system (see Table 44) The correct number of dampers na = 3. When na < 3. performance was extremely poor, indicating incorrect dy- namics due to missing dampers. When Tı > 3, performance was similar to the case for na 3. Redundant dampers provided the extra parameters, which some- times made optimization easier, but were often ignored by learning identical properties to existing dampers, by adopting zero damping coefficient, 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 's internal structure.
Table 4: Impact of # Com- ponents and VPT.
PODINNS
Training
Test
components, which are usu-
===0
0.005
0.000
1
1.009 0.001 10.000
na = 2
0.015
0.001100000 10.000
na 3
0.925 100007
0,58110.040
n. 4
0.932+40010
1.597+00088
Vhe 0.905 0.600:00.35
6 10 4-
See Appendix for additional visualizations and analyses.
CONCLUSION
In this study, we proposed Poisson-Dirac Neural Networks (PoDiNNs), which use a Dirac struc- ture 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 ex- periments with mechanical, rotational, electro- magnetic, and hydraulic systems validate these capa- bilities. Developing methods to address dynamical systems that are described by DAFs and partial differential equations remains for future research.
10

---

Published as a conference paper at ICLR 2025
ACKNOWLEDGMENTS
This study was partly supported by JST PRESTO (JPMUPR24TB, JPMJPR2IC7), CREST (JP- MJCR1914, JPMICR24Q5), ASPIRE (JPMJAP2329), Moonshot R&D (JPMIM12203--14) and ISPS KAKENHI (24K15105), and partly achieved through the use of SQUID at D3 Center, Os- aka University.
ETHICS STATEMENT
This study is purely focused on dynamical systems and modeling, and it is not expected to have any direct negative impact society individuals.
REPRODUCIBILITY STATEMENT
The environment, datasets, methods, evaluation metries, and other experimental settings are given in Section 4.J and Appendix [F] For full reproducibility, it is recommended to run the source code attached as supplementary material.
REFERENCES
Brandon Amos, Lei Xu, and J. Zico Kolter. Input Convex Neural Networks, In International Conference Machine Learning (ICML), 2017.
Anima Anandkumar, Kamyar Azizzadenesheli, Kaushik Bhattacharya, Nikola Kovachki, Zongyi Li, Burigede Liu, and Andrew Stuart. Neural Operator: Graph Kernel Network for Partial Differen- tial Equations. In ICLR 2020 Workshop on Integration of Deep Neural Models and Differential Equations (DeepDiffEq). February 2020.
Philipp Becker, Harit Pandya, Gregor H.W. Gebhardt, Cheng Zhao, and Gerhard Neumann. Re- current Kalman Networks: Factorized Inference in High-Dimensional Deep Feature Spaces. In International Conference on Machine Learning (ICML), 2018.
Petar Bevanda, Stefan Sosnowski, and Sandra Hirche. Koopman Operator Dynamical Models: Learning, Analysis and Control. Annual Reviews in Control, 52:197-212, 2021.
Aleksandar Botev, Andrew Jaegle, Peter Wirnsberger, Daniel Hennes, and Irina Higgins. Which priors matter? Benchmarking models for learning latent dynamics. Advances in Neural Infor- mation Processing Systems (NeurIPS) Track on Datasets and Benchmarks, 2021.
Steven L. Brunton, Joshua L. Proctor, and J. Nathan Kutz. Discovering governing equations from data by sparse identification of nonlinear dynamical systems. Proceedings of the National Academy of Sciences, 113(15):3932-399) 393333977, April 2016.
Ricky T. Q Chen, Yulia Rubanova, Jesse Bettencourt, and David Duvenaud. Neural Ordinary Dif- ferential Equations. In Advances Neural Information Processing Systems (NeurIPS), pp- 1-19, 2018.
Yuhan Chen, Takashi Matsubara, and Takaharu Yaguchi. Neural Symplectic Form
Learning Hamil- tonian Equations on General Coordinate Systems. In Advances in Neural Information Processing Systems (NeurIPS), 2021.
Yuhan Chen, Takashi Matsubara, and Takaharu Yaguchi. KAM Theory Meets Statistical Learning Theory: Hamiltonian Neural Networks with Non-Zero Training Loss. In AAAI Conference on Artificial Intelligence (AAAI), 2022.
Zhengdao Chen, Jianyu Zhang, Martin Arjovsky, and Léon Bottou. Symplectic Recurrent Neural Networks. In International Conference on Learning Representations (ICLR), pp. 1-23, 2020.
Leon O. Chua. Chua circuit. Scholarpedia, 2(10):1488, October 2007.

---

Published as a conference paper at ICLR 2025
Theodore James Courant. Dirac Manifolds. Transactions of the American Mathematical Society, 319(2):631, 1990.
Miles Cranmer, Sam Greydanus, Stephan Hoyer, Peter Battaglia, David Spergel, and Shirley Ho. Lagrangian Neural Networks. In ICLR 2020 Workshop Integration of Deep Neural Models and Differential Equations (DeepDifEgq), pp. 1-9, March 2020.
Luca Di Persio, Matthias Ehrhardt, and Sofia Rizzotto. Integrating Hot-HHmiltonian Systems with Neural Networks: From Deterministic to Stochastic Frameworks. arXiv, March 2024.
Laurent Dinh, Jascha International Conference on Learning Representations (ICLR), pp. 1-32, 2017.
and Samy Bengio. Density estimation using Real NVP. In
R. Dormand and P. Prince, A reconsideration of some embedded Runge-Kutta formulae. Journal of Computational and Applied Mathematics, 15(2):203-211 1986.
Yifan Du and Tamer A. Zaki. Evolutional deep neural network. Physical Review E, 104(4):045303, October 2021.
Vincent Duindam, Alessandro Macchelli, Stefano Stramigioli, and Herman Bruyninckx. Modeling and Control of Complex Physical Systems. Springer Berlin Heidelberg, Berlin, Heidelberg, 2009.
Sølve Eidnes, Alexander J. Stasik, Camilla Sterud, Eivind Bøhn, and Signe Riemer-Sorensee Pseudo-Hamiltonian Neural Networks with State-Dependent External Forces. Physica D: Non- linear Phenomena, 446: 133673, April 2023.
Christopher Eldred, François Gay-Balmaz, Sofia Huraka, and Vakhtang Putkaradze. Lie-Poisson Neural Networks (LPNets): Data-based computing of Hamiltonian systems with symmetries. Neural Networks, 173: 173:106162, May 2024.
Mare Finzi, Ke Alexander Wang, and Andrew Gordon Wilson. Simplifying Hamiltonian and La- grangian Neural Networks via Explicit Constraints. In Advances in Neural Information Process- ing Systems (NeurIPS), 2020.
François Gay-Balmaz and Hiroaki Yoshimura. Systems, variational principles and interconnections in non-equilibrium thermodynamics. Philosophical Transactions of the Royal Society A: Mathe- matical, Physical and Engineering Sciences, 381 (256))2222008, October 2023.
Sam Greydanus, Misko Dzamba, and Jason Yosinski. Hamiltonian Neural Networks. In Advances in Neural Information Processing Systems (NeurIPS), pp. 1-16, 2019.
Nate Gruver, Marc Finzi, Samuel Stanton, and Andrew Gordon Wilson. Deconstructing the Induc- tive Biases of Hamiltonian Neural Networks. In International Conference Learning Represen- tations (ICLR). arXiv, February 2022.
Albert Gu and Tri Dao. Mamba: Linear- Liner-Time Sequence Modeling with Selective State Spaces. In Conference On Language Modeling (COLM), 2024.
Albert Gu, Karan Goel, and Christopher Re. Efficiently Modeling Long Sequences with Structured
State Spaces. In International Conference on Learning Representations (ICLR), 2022.
Ernst Hairer, Christian Lubich, and Gerhard Wanner. Geometric Numerical Integration: Structure-
Preserving Algorithms for Ordinary Differential Equations, volume 31. Springer-Verr Berlin/Heidelberg, 2006.
Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun. Deep Residual Learning for Image Recognition. In IEEE Conference Computer Vision and Pattern Recognition (CVPR), pp. 1-9, December 2016.
Sepp Hochreiter and J'rgen Schmidhuber. Long Short-Term Memory. Neural Computation, 9(8): 1735-1780, November 1997.
Masanobu Horie, Naoki Morita, Toshiaki Hishinuma, Yu Ihara, and Naoto Mitsume. Isometric Transformation Invariant and Equivariant Graph Convolutional Networks. In International Con- ference Learning Representations (ICLR), 2021.
12

---

Published as a conference paper at ICLR 2025
Eugene M. Izhikevich and Richard FitzHugh. hitpr/ddarrrrgaarriieee/ttttahha _mogel, 2006.
FitzHugh-Nagumo model.
Pengzhan Jin, Aiqing Zhu, George Em Karniadakis, and Yifa Tang. Symplectic networks: Intrinsic struecure-preserving networks for identifying Hamiltonian systems. Neural Networks, 132:166- 179, 2020.
Pengzhan Jin, Zhen Zhang, loannis G. Kevrekidis, and George Em Karniadakis. Learning Poisson Systems and Trajectories of Autonomous Systems via Poisson Neural Networks. IEEE Transac- tions on Neural Networks and Learning Systems, pp. 1-13, 2022.
Muhammad Firmansyah Kasim and Yi Heng Lim. Constants of motion network. In Advances in Neural Information Processing Systems (NeurIPS), arXiv, August 2022.
Hassan K. Khalil. Nonlinear Systems, Third Edition. Prentice Hall, 2002.
Patrick Kidger, James Morrill, James Foster, and Terry Lyons. Neural Controlled Differential Equa- tions for Irregular Time Series. In Advances in Neural Information Processing Systems (NeurIPS), volume 33, pp.
Curran Associates, Inc., 2020.
Diederik P. Kingma and Jimmy Ba. Adam: A Method for Stochastic Optimization. In International Conference on Learning Representations (ICLR), pp. 1-15, December 2015.
Rahul G. Krishnan, Uri Shalit, and David Sontag. Structured Inference Networks for Nonlinear State Space Models. In AAAI Conference on Artificial Intelligence (AAAI), pp. 1-13, 2017.
Remi Lam, Alvaro Sanchez-Gonzalee Matthew Willson, Peter Wirnsberger, Meire Fortunato, Fer- ran Alct, Suman Ravuri, Timo Ewalds, Zach Eaton-Rosen, Weihua Hu, Alexander Merose, Stephan Hoyer, George Holland, Oriol Vinyals, Jacklynn Stott, Alexander Pritzel, Shakir Mo- hamed, and Peter Battaglia. Learning skillful medium-range global weather forecasting. Science, 382(6677)):416-142) December 2023.
Ilya Loshchilov and Frank Hutter. SGDR: national Conference on Learning Representations (ICLR), pp. 1-16, 2017.
Stochastic gradient descent with warm restarts. In Inter-
Jerrold E. Marsden and Tudor S. Ratiu. Introduction Mechanics and Symmetry: A Basic Exposi- tion of Classical Mechanical Systems, volume 17 of Texts in Applied Mathematics. Springer New York, New York, NY, 1999.
Takashi Matsubara and Takaharu Yaguchi. FINDE: Neural Differential Equations for Finding and Preserving Invariant Quantities. In International Conference on Learning Representations (ICLR), 2023.
Takashi Matsubara, Ai Ishikawa, and Takaharu Yaguchi. Deep Energy-Based Modeling of Discrete- Time Physies. In Advances in Neural Information Processing Systems (NeurIPS), 2020.
Cyrus Neary and Ufuk Topcu. Compositional Learning of Dynamical System Models Using Port- Hamiltonian Neural Networks. In Proceedings of The 5th Annual Learning for Dynamics and Control Conference, pp. 679-691. PMLR, June 2023.
Oliver Nelles. Nonlinear System Identification. Springer Berlin Heidelberg, Berlin, Heidelberg, 2001.
Adam Paszke, Gregory Chanan, Zeming Lin, Sam Gross, Edward Yang, Luca Antiga, and Zachary Devito. Automatic differentiation in PyTorch. In NeurIPS2017 Workshop on Autodiff. pp. 1-4, 2017.
Tobias Pfaff, Meire Fortunato, Alvaro Sanchez-Gonzaler and Peter W. Battaglia. Learning Mesh- Based Simulation with Graph Networks. pp. 1-18, 2020.
M. Raissi, P. Perdikaris, and G. E. Kariiddakss. Physics-informes neural networks: A deep leamn ing framework for solving forward and inverse problems involving nonlinear partial differential equations. Journal of Computational Physics, 378:686-707, February 2019.
13

---

Published as a conference paper at ICLR 2025
Justin Sirignano and Konstantinos Spiliopoulos. DGM: A deep learning algorithm for solving partial differential equations. Journal of Computational Physics, 375:1339-1130 December 2018.
Arjan van der Schaft. Implicit Hamiltonian systems with symmetry. Reports on Mathematical Physics, 41(2):203-221, 1998.
Arjan van der Schaft and Dimitri Jeltsema. Port-Hamiltonian Systems Theory: An Introductory Overview. Foundations and Trends@ in Systems and Control, 1(2):173-378, 2014.
Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N. Gomez, Lukasz Kaiser, and Illia Polosukhin. Attention is all you need. In Advances Neural Infor- mation Processing Systems (NIPS), 2017.
R. Vlachas, J. Pathak, B. R. Hunt, T. P. Sapsis, M. Girvan, E. Ott, and P. Koumoutsakos. Backprop- agation algorithms and Reservoir Computing in Recurrent Neural Networks for the forecasting of complex spatiotemporal dynamics. Neural Networks, 12::191-217, 2020.
Yi Jen Wang and Chin Teng Lin. Runge-Kutta neural network for identification of dynamical sys- tems in high accuracy. IEEE Transactions on Neural Networks, 9())294-337, 1998.
Antoine Wehenkel and Gilles Louppe. Unconstrained Monotonic Neural Networks. In Advances in Neural Information Processing Systems (NeurIPS). arXiv, 2019.
Hiroaki Yoshimura and Jerrold E. Marsden. Dirac structures in Lagrangian mechanics Part I: Im- plicit Lagrangian systems. Journal of Geometry and Physics, 57(1):133-156, 2006a.
Hiroaki Yoshimura and Jerrold Marsden. Dirac structures in Lagrangian mechanics Part Vari- ational structures. Journal of Geometry and Physics, 57(1):209-22) -5000 2006b.
Yaofeng Desmond Zhong, Biswadip Dey, and Amit Chakraborty. Dissipative SymODEN: Encoding Hamiltonian Dynamics with Dissipation and Control into Deep Learning. In ICLR 2020 Workshop Integration of Deep Neural Models and Differential Equations (DeepDiiffgg) pp. 1-6, 2020.
14

---

Published as a conference paper at ICLR 2025
Appendix
A LESS RELATED WORK
Research on methods for learning dynamical systems has long history (Nelles 2001p. Although the approaches discussed here are not directly related to the primary focus of this paper, we include them for completeness and to clarify key distinctions. Since linear systems are relatively straightforward to handle, our discussion is restricted to nonlinear systems.
classical approach to learning nonlinear dynamics considers discrete-time systems of the form Tn+1 = fzz), where the mapping f is approximated using neural network (Nelles 2001). Even when the underlying system is continuous in time, one can sample the trajectory at fixed intervals and treat it as discrete-time dynamics. A variant of this approach models the dynamics as Zn+11 En + (x which closely resembles the forward Euler method. Building on this idea, Wang & Lin (1998) proposed neural network architecture inspired by the fourth-order Runge- Kutta method. Neural ODEs can be regarded as an extension of this approach although the original work noted connections to ResNet (He et al., (016).
If state x is only partially observable, the dynamics is typically modeled as a combination of latent- state dynamics and an observation model y = ((x), which is commonly referred as a state-space model (SSM), Among neural architectures, long short-term memories (LSTMs) can be interpreted as discrete-time SSMs (Hochreiter & Schmidhuber 1997). More sophisticated models combine neural networks with Bayesian inference to handle uncertainty in latent states (Krishnan et al. [0117] Becker
et al. 2018). Recently, SSMs have been employed in sequence modeling tasks such as language
processing, leading to architectures like S4 and Mamba (Gu et al. [022] Gu & Dao. 2024). These
models are primarily designed for modeling logical reasoning or expressing causal relationships
rather than for modeling physical dynamical systems.
A different viewpoint on dynamical systems leverages the Koopman operator K,, which describes the time evolution of observations rather than the latent states x themselves (Bevanda al. (021). Formally, the operator IC is linear (albeit typically fiinite-dimensionall and satisfies y(t + s) = KC* (²)())).. The Koopman operator provides an alternative perspective for analysis, rather than for modeling.
Another framework that is sometimes confused with our primary focus is that of Physics-Informed Neural Networks (PINNs) (Sirignano & Spiliopoulos) 2018 Raissi et al.) 2019 Du & Zaki (001)) PINNs are designed to solve known symbolic governing equations by employing neural networks as basis functions, analogous to classical numerical methods such as the finite element method. Although PINNs can be used to estimate certain parameters within symbolic equations from data, they are not intended for the discovery of unknown dynamical systems.
For directly learning symbolic governing equations from data, the Sparse Identification of Nonlinear Dynamics (SINDy) framework has been extensively investigated (Brunton et al. 2016). SINDy typi- cally represents the dynamics as a combination of predefined terms in a dictionary, which introduces a trade-off between expressiveness and interpretabilityy
As summarized in T the research trend which our work contributes is focused on data-driven modeling of dynamical systems. We assume that the governing symbolic equations are unknown and not the primary target for estimation, and that the full state of the system observable. While combining the proposed method with the aforementioned approaches could be promising direction for future research, it is not meaningful to directly compare their performance.
Moreover, our study specifically considers dynamical systems as coupled systems across multiple physical domains and focuses on learning and uncovering both the coupling patterns and the char- acteristics of individual components. This distinguishes our work from previous research. Although system identification for coupled linear systems has been extensively studied (Nelles) 2001), to the best of our knowledge, prior work has tackled flexible identification using neural networks.
15

---

Published as a conference paper at ICLR 2025
B EXAMPLES OF FORMULATIONS
Consider the time evolution of a point u on a manifold M. af local coordinate on M denoted by Ig the corresponding basis vector of the tangent space TuM is denoted by ău," M is expressed as X = Σ; X; Uz. rate of change of the state (t) in the z,-direction, i.c., 'ي = X,. For a function F M + R, its differential d/' is given by df Σ, aF dz,, The differential a covector field, that is, a collection dh;
0_ A vector field X on
a_
the curve (t) satisfies ú(t) = Xu(t) X; represents the local
of points on cotangent spaces T.M.. The basis vector of the cotangent space Tw- R is dz,, and the 8_ following relation holds: dr; Øzz = 1 if = and 0 otherwise.
When a function F defines dynamics, 2-tensor field links a vector field X and the covector field dF. Such 2-tensor field can be symplectic form, Poisson bivector, Riemannian metric. For example, link with the negative of a Riemannian metric defines gradient flow. A symplectic form and Poisson bivector are defined using the wedge product A, which satisfies the skew-svmmetry dr; A dxy = - dx, A dr; and the relation (dæ; A dzz) da; dxy. Marsden & Ratiu (1999) and Hairer et al. (206) have thoroughly discussed how to describe dynamical systems using these geometric objects. While theoretical details are left these textbooks, this section introduces a mass-spring system as concrete example.
B.1 MASS- SYSTEM
Consider mass-spring system with spring constant is and mass m. We will write its dynamics by several formulations.
Canonical Hamiltonian Systems In the Darboux coordinates (i.e., the cotangent bundle T Q), the generalized coordinate q is the displacement of spring k. and the generalized momentum is obtained as p me for the velocity of mass m. The manifold M is 2-dimensional Euclidean space R22 The Hamiltonian I is H (g. i.e. Ω dg dp.
2m p² + kkg² and the symplectic form Ω is standard,
Hamilton's 's equations state that (X+) = (dq A dp) B + + = údp jdg equals dH aH dg + aH 3p
aH = p/ m, 1 aH dp, leading to the equations of motion, q = op
-g.
0q
Oq
Non- Non-Canonical Hamiltonian Systems On the tangent bundle A TQ, the velocity e used in place of the momentum p, and the symplectic form Q is the Lagrangian 2-form S m(dg A de).
The Hamiltonian is H (g, v) = =m=² + kkq².
Hamilton's 's equations state that (X))(H+)) = m(dg A du) 8 by
+ ( = midv - modg equals
dHf
BM dg + a1l du, leading to the equations of motion, ÿ = m dt v, ú de
1 aM - ky/ m.
0q
m 0g
Poisson Systems In the Darboux coordinates, the Poisson bivector is B
0.
A 0
Hamilton's 's's
0p
dg
equations state that Xµ 9 ag
Ø + P ap equals (dHH) = ( 0 0
Op
A
0q
(
aH dg +
Bq
bd dp)
Op
08
Öp 0g
8H 8 The equations of motion are q Bp Bq Bp
aH
p/mm m. i
2H
kg.
Oq
On the tangent bundle M TQ, the Poisson bivector is B
( ^ 0g and the Hamiltonian m 0v
is H (q, v) me² + 1 =g Hamilton's equations state that X = 0 + D0
à
0 equals (dHH)
1( 1
A 2. 0H
( 3g
dg + Te du)
aH
1 0H 1 0H 2. The equations of motion are ğ = 1 2H
m De
0y
m de 0g m og dv
de
0, m dq lq/m.
aH
B.2 CONSTRAINED DAAS-NG SYSTEM
Consider pair of mass-spring systems, indexed by E (1, and introduce a constraint such that the two masses are coupled and always have the same displacement and velocity; the dynamics degenerate. This constraint is expressed as qi Q2. The Hamiltonian is given simply by the sum of two systems, H(q.p) _,, of + 2mg + + k + f2242. However, Hamilton's 's equations Ω (X.)) dH with the standard symplectic form Ω cannot describe the dynamics.
16

---

Published as a conference paper at ICLR 2025
Degenerate Systems with Coordinate Transformation Define new coordinates q = Q1 92 and = P2, and consider the submanifold M spanned by P). The Hamiltonian H is unchanged but rewritten as H (q.P) - 2[m "mm-+my D² + (k +k2)q². With the standard symplec-
p 122
mi P1
m2
tic form = dg A dp on M, the equations motion are written as q mi+my P = -- -k: + k-)
Degenerate Systems with Constraint Force The constraint on the coordinates, 91 92 0. is naturally satisfied by the constraint on the velocities, 91 42 0, which is written as the constraint on momenta, P1/ p1/mp P2/m² O. The constants of motion are q1 92 and p1/ma - P2/ P2/m2. A constraint force can be introduced to satisfy these constraints, yielding the same equation as above. Finzi et al. (2020) 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 Poisson Systems The equations of motion on the submanifold M, q
P
P = -Kk + k2)g, can be rewritten in the original coordinate system as (1 = _21 P22 mR
mm++mng 92 =
mi mi++m2? 111 mi mmmm (kiq1 + k2g2), f12 mi+ma (kıq) + k242). Even on the original manifold
P1 PLtP2
M. these equations are obtained from Hamilton's 's equations X = (dHH) with the Poisson bivector 2)) A (8n + 3oz
2) on M, while its bundle map B 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.
B my=ma (m1 Bp1
a + m2 8y2
On the tangent bundle M = 7Q, the Hamiltonian is H(q, 92. 01, U2) = m
+
2 m2 +
kig + 2 and the Poisson bivector B is B equations of motion are á1
+102
Joz
8
A
0
+
84z
8.
Then, the
my 12
02
278 1. +ma01
01
RigI
+
02
Ang H2
+m2g
ma mm
C POISSON SYSTEMS
C.I 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 E (1,2). coupled sequence from fixed wall. Let qi denote the displacement of i-th spring. The equations of motion are given by d1 = P1/m1, Ú2 = P2/ m2 P1 P1/m/.111 P1
- + ke12. and /
/2- = -242- The bivector B leading to these equations is B Øp ^ Tyi pp 0gz 0pz Jaz for the Hamiltonian H (q. p) = 2
8
1 př + 2m, p + k2 + 2242 These terms indicate the couplings
between PL and qı P1 and 92. and p2 and 92- The negative coefficient indicates that the coupling between P1 and 92 is in the opposite direction.
Degenerate Systems as Poisson Systems As shown Appendix B.2 degenerate systems can be expressed as Poisson systems. Recall the case of pair of mass-spring systems, indexed by i E (1,2), constrained so that the two masses always have the same displacement and velocity.
This system is expressed with the Poisson bivector B mitimg (m Øp1
9 + m2 6pp (yy
8 m2 ( + dep
which is degenerate in the sense that its bundle map B is degenerate. This indicates the absence of a corresponding symplectic form. Conversely, by learning the bivector B and examining how degenerates, one can identify the constraints imposed on the system.
The coefficient spring Q1- The effort e âqi
mi for the term
a A a indicates the coupling strength between mass pl and my m1-ma Cgn goes
me tm
Upi aH from spring Q1 distributed to the masses such that
ogi
to mass mi and
my Cy, goes to mass m2.
m22
Coordinate Transformation by Bivector As shown in Appendix B.2 the elements of the Poisson bivector B in Poisson system depend whether the dynamics are defined on the tangent bundle TQ (using generalized velocities as part of the state) or the cotangent bundle T Q (using generalized momenta part of the state). Despite this difference, both coordinate systems represented by Poisson bivectors. Thus, by learning the bivector B to approximate the dynamics a given system, one can implicitly learn the coordinate system employed by that system.
17

---

Published as a conference paper at ICLR 2025
C.2 ELECTRIC CIRCUITS
Our formulation is applicable to electric circuits, as summarized in Table Let Xx and Vx denote the current through and voltage across circuit element X, respectively.
capacitor with capacitance C has the electric charge Q as its state, and stores the energy Hc = 20 Q² Its flow is the current Ic. which leads to the change in the electric charge Q as Q Ic. Its
effort is the voltage Vo because the stored electric charge Q generates the voltage Vo as Vc g aHc. 8Q
inductor with inductance has the magnetic flux yp as its state, and stores the energy H 21 .². Its flow is the voltage VL, which leads to the change in the magnetic flux &p as V1- Its effort is the current l because the magnetic flux vp generates the current [7 IL - ap
0tr,
Consider system composed of an inductor L and capacitor C' coupled in series. The state space M is the space of the magnetic flux and electric charge Q. The total energy is H = HL + Hc, 0 and its differential is dH = ddp + EdQ. Define a bivector B = ai
A 8 oQ. which leads to
10 0
Q à
B$
LaQ
C 0.p
The vector field on M is X equations of motion:
600 + Q σq
Hamilton's equations X (dHH) lead to the
Q Vo and lo Q VE C
However, electrical circuits that can be described as Poisson systems are limited to energy- conservative LC circuits. Our formulation extends this include resistors, diodes, voltage sources, and current sources.
D DIRAC STRUCTURE
D.1 DIRAC STRUCTURE ON A VECTOR SPACE
Proof of Theorem Yoshimura & Marsden 2006a). By Definition I
Dj (3 E V (ax, un) + () 0 for all E A and A.).
Let (Ü, a) E Dv. Then, ö E A and a E A so (a, ū) + (ā, v) = 0 for all (v, E Dv. This implies (Ø. a) E DV Thus, Dy C DV.
Let (w, (3) E Dt. In the above definition of Dj. setting a gives v3 0 for all A. Hence, B E A", Similarly, setting = 0 gives (a, w) = for all a E A which implies w E A. Hence, (w, B) E Dy. Thus, D C Dv.
Therefore, Dy DV-
D.2 DIRAC STRUCTURE ON MANIFOLD
Definition 4 (Courant 1990); Yoshimura & Marsden 2006a)). A distribution A on a manifold M is a collection of vector subspaces Du of langent spaces TuM, each assigned smoothly M at point u, forming a vector subbundle of T.M The annihilator A of A is a collection of the annihilator A of D also forming a subbundle T M. Then, a Dirac structure constructed as D = A + o, which is subbundle of the Pontryagin bundle TM @ M.
A type of Dirac structure on manifold M can be defined using symplectic form Ω on M. Theorem (Yoshimura Marsden (2006a)). Let AM be a distribution on a manifold M. Let P be symplectic form a manifold M. Ω is restricted to AM and denoted by Usw Then, the collection of
(Dax)u ((v,0) E T.M M ll E (AM)us
and (a, w) = (Asn))(( (w,,w) for all w E (AM)u}
18

---

Published as a conference paper at ICLR 2025
is Dirac structure DAM C TTMOT M M, which is said to be induced by Q2
Proof of Theorem Yoshimura & Marsden M is given by
2006a). By definition, the orthogonal of DAM at u E
(D (Dmnu= [(w, B) E Tu x T (ar, w) + (B, v) = 0 for all E (AM)u.
and (a,, w) = (alam)((t) w) for all w E (AM)u).
Let (v, xx; (ū, (Dam)u. (a, + ( v) Sam(m) (u, 8) Ram (Ü. v) 0. The latter equality holds because of the skew-symmetry f. This implies (v, E (DaM)) bm)d. Thus, (Dam)u C (Džm)u
Let (w, (3 (Dam)u. Then, (a, w) + (8, v) 0 for all (8, (Dam)u. First, setting 0 gives (a,w) (S )u (0, w) for any E (AM).. Thus, w (AM) (Am)u- Second, if a satisfies (a. w) = (L,m)u(v, w) for any v E (AM)u. then (a, w) + (3, v) = (l w) + (B.u) = for any v (AM)u. This implies (B,0) (Dam)u (w, v) for any v E (AM)u. Because it has been proved that w E (AM)us 3) E (Dam)u- Thus, (Dim)u (Dam)u.
Therefore, (Dam)u (D
Definition 5 (Hamilton-Diree System (Yoshimura & Marsden 2006b)), Let M be manifold, H M + I be an energy function, and D be a Dirac structure on M. Given a vector field X on M, if it holds for each (1) .M and E [a, b] that
((Xn)ut))
(Dam)u(0)
the tuple (H, Dams XH) H is called Hamilton-Dire system (or implicit Hamiltonian system).
Note that the symplectic form restricted to a distribution AM. denoted by faam. satisfies (Dam)) (w, w) = Pu (v, w) for any v, (AM)u. The symplectic form Ω2 is degenerate in the sense that its bundle map 52 on the tangent space TuM is degenerate, but Ωam is non-degenerate on the distribution (A (Am).. u. If there is no constraint, AM = TMM and fAm f2. When the Dirac structure D is given in Theorem
3 with no constraint, the Hamilton-Dirac system is identical to the Hamiltonian system 22' (Xx) HH in Eq. ( The distribution AM describes constraints on the velocities, and hence constrained Hamiltonian systems can be rewritten as Hamilton- Dirac systems. Therefore, Hamiltonian-Dirs systems are generalizations of Hamiltonian and constrained Hamiltonian systems. For example, the system in Appendix B.2 can be written as a Hamilton- -iirac system follows:
Degenerate Systems as Hamilton-Diral Systems The distribution AM is given by (AM) {(d,, (r, E TuM ④1 - 42 = 0, p1/m1 //// 12/m2 = 0} = span de m2 8ps }. Define the new coordinates q q1 q2 and p m
à
+ Beg mi 0pi + m1 p2, with m2
mi p1
AM s à à ). The symplectic form Ram restricted to the distribution AM should satisfy By 0p
(MaM)) (w, w) = S2 (v, w) for any v, E (AM)u at u. Such form given las = dg dp. Then, the equations of motion are written as I = mitmy" P (ka + k2)g.
D.3 PORT-BASED SYSTEMS
Proof of Theorem
2 By definition, the orthogonal of D at u E M is given by
D = ((f, En, (e. f) + (e, f) = for e E En and f = B (e)).
Let (f, e), (f. E Du. Then, () + (e. f) = (e, (e(()) + (e, Bt = (u(e, +
Bu(e, due to the skew-symmen sommeeyyy of B. This implies (f. e) D., Thus, D. C D
Let (f,c) E D.,- Then, for any e) E Dn (i.., for any e e En and f (e). 0 f)+ (e. (e. f) + (e, Bf
(c. + e, (c. B, C) (e. f) (e, Bi (e, Bu é) for all e E E, implies that f = BE(e), Thus, D C Du.
Therefore, D., D4.
19

---

Published as a conference paper at ICLR 2025
By extending the above Hamilton- Dirac system from the Pontryagin bundle TM + T to the vector bundle F + E, we obtain the so-called Nort-Hamiltoniall systems (Courant) 1990, van der Schaft & Jeltsema| 2014). However, previous neural network-basse methods employed the canoni- cal form in Eq.. 44 and have not attempted to identify coupling patterns or component wise charac-
teristics (Zhong al. 2020 Eidnes et al. 2023
Neary & Topeu 2023
Di Persio et al. 2024
In bond graph representation of 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 al. 2009).
Instead of the symplectic form Ω our formulation employs (Poisson) bivector B to define the Dirac structure D on the tangent and cotangent bundles (T.M J= and T { plus the vector bundles for port variables (FR Fí, E and E¹). As is the case with the symplectic form Ω we can constrain a Poisson bivector on codistribution A which is a subbundle of the cotangent bundle 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 [ Note that, for a bivector to be called Poisson bivector. it must satisfy certain conditions, such as the Jacobi identity (Courant) [990). However, the learned bivector in our formulation does not necessarily satisfy these conditions, nor are these conditions required for the proof' presented earlier. Thus, throughout this paper, we refer to B simply as a bivector, without restricting it to being Poisson bivector.
E IMPLEMENTATION
E. COORDINATE-BASED FORM AND REDUCTION TO EXISTING MODELS
In the main body, we have introduced PoDiNNs in a coordinate-free form for a general and theoret- ical derivation. However, by introducing a specific coordinate system, we can rewrite PoDiNNs in a more intuitive form at the cost of generality. In this section, we will clarify the relationship with previous work through such rewrite.
As stated in Theorem 2 PoDiNNs use a bivector to define the Dirac structure D as subbundle of F @ E. When local coordinates are introduced on M, the corresponding basis vectors for the tangent space T.M (that is, FFS) are a. Additionally, the basis vectors for FR and F are Jz." specified as <. and ç/, respectively. The basis vectors for & are then defined as their duals. The bivector B at u can expressed linear combination of wedge products of pairs of basis vectors for F, F FR $ F!
We assume that the state space M is a Euclidean space. Then, at a specific point u, flows I and efforts e can be represented as vectors, and the linear bundle map Bi, is represented as a skew- symmetric matrix B E Rm where m denotes the total number of dimensions of Fu (Marsden & Ratiu 1999). A skew-symmetric matrix M satisfies the property M -M" which implics that the maximum number of independent degrees of freedom is m - 1). This characteristic arises from the skew-symmen symmetry of the wedge product, particularly the relation 0 A a 0z, 8w: 0w:
0 A Ozy"
0
Each bivector element corresponds specific pair of elements within this matrix.
In this context, the equation f = (e) in Theorem can be expressed with a matrix-vector product:
fs
S
ftt f
Bi
ft e
By substituting the variables as in Definition
we obtain:
[uH(u)
f R
B
Ru Rf)
f
(t)
where the differential d. reduces to the gradient
V, H in Fuclidean space.
20

---

Published as a conference paper at ICLR 2025
In the absence of dissipation (R) and external inputs (I), the equation reduces
d V.HH)
This form aligns with those derived in NSFs and CHNNs, up to certain constraints imposed on the matrix (Chen et al. 2021; Finzi et al. 2020). this sense, PoDiNNs are generalization of these previous models. Since NSFs are generalization of LNNs, PoDiNNs are also generalization of LNNs.
Furthermore, if the Darboux coordinates are introduced on M, the state can be decomposed into the position q and momentum and the bivector B takes the standard form Σ. Bp; A 0g." The equation reduces to:
O - p -I
V P) o VpH(q.
where 2 = [3]- This is identical to Hamilton's equations in Eq. (.. which is the form HNNs employed. Therefore, PoDiNNs are generalization of HNNS. If the target system is known to be learnable by HNNs. it is also learnable by PoDiNNs of the form identical to HNNs. Consequently, a performance comparison between PoDiNNs and HNNs is not meaningful.
Let us consider a specific example from Appendix C72 system where a capacitor C' andnd an inductor are coupled in series. The state of the capacitor C' is its electric charge Q, while the state of the inductor is its magnetic flux sp Their coupling is represented by the bivector B a. ^ aQ
Using local coordinates (Q, :), we can represent B as the matrix B [ o] The Poisson-Dirae
formulation states that applying this matrix to the effort (VqHc, V. yields the flows (Q, :p), that is,
1c Vr
0 -1
Volle WHHH
1r -Vo]
This corresponds to the standard form of Hamilton's equations. We can see that the first line rep- resents Kirchhoff's current law, Ic Jr., and the second line represents Kirchhoff's voltage law, V, = -Vc. The sign inversion follows from the skew-symmetric property of the bivector B. As the number of components increases, the dimensionality the flows, efforts, and matrix also grows, and the bivector may not retain a standard form, as observed in system (c). Nonetheless, the bivector consistently represents Kirchhoff's laws, independent of the component- characteristics, and the Poisson-Dirae formulation ensures the satisfaction of these laws. Similarly, in mechanical systems, bivectors represent the law of action and reaction.
In general, the coordinate systems are not limited to the Darboux coordinates, components can exhibit a variety of coupling patterns, and the matrix Bi can contain elements beyond the restricted values of -1, 0, or 1. The state u is not confined to mechanical position and momentum; it may also be decomposed into various quantitics such as electric charge, voltage, or the volume of fluid in a tank, depending the domain from which these states are derived, shown in Table
E.2 PERFORMANCE AND COMPUTATIONAL COST
In the above case, PoDiNNs use two neural networks to approximate two energy functions Hc and Hr R + R. thereby forming a mapping that defines an ODE, H22 - IR2 (Q, + (Q,+) (Ic, V2). PoDiNNs introduce the assumption that the target system can be decomposed into inter- acting components and use separate neural networks R + R to model n individual components. This approach distinct from NODEs, which directly learn the mapping R² + R without explicit assumptions between variables. In this way. PoDiNNs can mitigate the curse of dimensionality. Moreover, while neural networks used in NODEs are universal approximators capable of learning Kirchhoff's laws and the law of action and reaction, they inherently introduce approximation errors. This means that these laws are only approximately satisfied, and over time, these errors tend to ac- cumulate, becoming significant. In contrast, PoDiNNs explicitly encode these laws using bivector elements, independently of the component-wis characteristics. This explicit separation enhances robustness.
the target system includes n components, we need to compute n neural networks. along with the bivector. This increases the computational cost times that of single NODE, We view
21

---

Published as a conference paper at ICLR 2025
this computational overhead as reasonable trade-off for the improved accuracy and interpretability provided by PoDiNNs.
F DATASETS AND EXPERIMENTAL SETTINGS
F.1 IMPLEMENTATION DETAILS
Following previous studies (Greydanus et al., 2019 Matsubara et al. 2020). 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 hyperbolic tangent activation function. Weight matrices were initialized using Py PyTorch's ' default algorithm.
In PoDiNNs, each element of the bivector B related to energy- dissipating components initial- ized from 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 A = 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- Cormand-Prince method (dopri5) with absolute tolerance atol = 10 and relative tolerance rtol= 10 was used to integrate the ground truth ODEs and neural network models (Dormand & Prince 1986). The Adam optimization algorithm (Kingma & Ba 2015) was applied with parameters (3,, [2) = (0.9, 0.999) and batch size of 100 for updates. The learning rate was initialized at 10-3
and decayed to zero using cosine annealing (Loshchilov & Hutter (017). The number of training
iterations was set to 100,000.
All experiments were conducted on single NVIDIA A100 GPU.
F.2 MECHANICAL SYSTEMS
Overview and Experimental Setting spring generates a force that restores it to its original length when stretched or compressed. Specifically, this force is given by e = -k(g)). where Aq the spring's displacement. In both the Hamiltonian and our formulations, the potential energy U is obtained by integrating this force e with respect to the displacement Aq. Conversely, the force e is obtained as to the partial derivative of the potential energy U with respect to the displacement Ag. The spring's flow f is the rate of extension, Aq.
mass m has velocity e, and its momentum given by p = me. The kinetic energy is 1mm? 2m p. In general, velocity u is more easily observed than momentum as the state of a mass. Therefore, we provided velocity
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 & & A 0
Specifically, the following couplings are possible: mass and spring Sp. zuj" mass and damper op. ^ E,C, mass and external force & Up. A Ej. moving wall and spring s! 0g. and moving wall and damper ç/ E,. These combinations were incorporated as elements of the learnable bivector B, while all other possible pairings were fixed to zero
a
In the absolute coordinate system, the positions y are used to represent the states of the springs, and their displacements Ag; are then calculated within the potential energy function U. Therefore, the potential energy function U depends on all positions q, of springs and moving boundaries collec- tively. In the relative coordinate system, the displacements Aq are used instead, and in this case, the potential energy function U; is defined individually for each displacement Ag:-
In both Dissipative SymODENs and PoDiNNs, since we assume to know the symbolic expression _pp for the kinetic energy, this used, rather than neural network. Also, the velocity was transformed into momentum p using the learnable parameter m with = mu. This approach was originally adopted by Dissipative SyMODENs and is a realistic and practical choice.
22

---

Published as a conference paper at ICLR 2025
For Dissipative SymODENs and PoDiNNs in the absolute coordinate system, a single neural net- work was trained to approximate the potential energy function U using all positions qi of springs and moving boundaries as inputs. In the relative coordinate system, separate neural network U, was used for the potential energy of each spring ki-
For Dissipative SymODENs, the dissipative terms D and input gain G' were modeled using neural networks. In the original paper (Zhong et al.] 2020). 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.
System with External Force This system consists of three springs h; and three masses my arranged sequentially and indexed by i E 21, from a fixed wall. Two dampers d1 and d3 are placed in parallel with springs 11 and k3, respectively. An external force F is applied to mass m3.
(a) Mass- Maas--pring-Dampeer
The masses were set to m, 0.8 + 0.2i. The characteristics of the nonlinear spring were given by (A4:))) = + + 0.1 The characteristics of the nonlinear dampers were defined as d.(:;) = (0.1 0.02i) sgn(1;) (B:))1ux) and sgn the sign function that returns for a positive value and -1 for a negative value. The ini- tial positions of the springs and the initial velocities of the masses were sampled from the uniform distributions (-0,5, 0.5) and U--0.3, 0.3), respectively. The external force F Ft) = e (f) was de- fined as the sum of three sine waves, with each wave's 's amplitude, angular velocity, and initial phase sampled from the uniform distributions U (0.2, 0.5), (0.1 0.2x). and (O, 2x), respectively.
where i denotes the extension velocity of the i-th damper,
System with Moving Boundary This system consists of three springs k and three masses m, arranged sequentially and indexed by E (1,2.3} from a moving wall b. Three dampers d are placed in parallel with springs kg This potentially represents building's response during an earthquake.
(b) Mass--Sring-Damaae
The masses were set to m, = 1.6 - 0.2i. The characteristies of nonlinear spring were given by k; = (0.6 - 0. 0.1i) + 0.1 where Aq: denotes the displacement of the i-th spring. The characteristics of nonlinear damper were defined as (1; = d, sgn(v;)) ee//// for di = 0.10, d2 = 0.05, and d3 0.02. The initial positions of springs k; and the initial velocities of masses my were sampled from the uniform distributions U-,.5, 0.5) and U(-0.3, 0.3), respectively. The position (t) of moving wall b was set as the sum of three sine waves, with each wave's 's amplitude, angular velocity, and initial phase sampled from the uniform distributions U (0.2, 0.4). U (0.05T, 0.2r), and (0, 2Tr., respectively.
In the absolute coordinate system, it is necessary to represent the potential energy of spring k connected to the moving wall b. For both Dissipative SymODENs and PoDiNNs, the position q of the moving wall was included as part of the inputs to the neural network U for the potential energy. along with the positions qi of all springs k. It was also fed Neural ODEs.
The effort el of the moving wall was its velocity thh 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 d1, connected to moving wall b.
(c) Mass-Spring System with Redundancy This system consists of five springs k; and two masses m, arranged in 2-dimensional space, as shown in Fig. (c). Masses my and m2 are con- nected to fixed wall via springs k and h2, respectively. These masses are also connected to each other by spring hz. Additionally, springs ky and Kş diagonally connect masses m2 and mi to the fixed wall, respectively, similar to cross braces. Since there are energy disegg-dissipating com- ponents or external inputs, the total energy is conserved. The masses were set to m. = 5.0 and m2 = 3.0. The natural lengths of spring k and k2 to ( = 3.0, that of Ka to la 4.0, and those of kg and ks to 14 = 5.0. The characteristics of nonlinear spring were ky (Aq) = 2.54g + 3.4223, 3.5 + 2.4Ag " ks (Ag) 2.52q + 1.64ga for the displacement Ag. The initial positions of masses m, and my were sam- pled from uniform distributions U (-0.5, 0.5) in both x and y directions, and their velocities from U.-0.100.1)
kg(Ag) 3.0Ag + 0.5Aq², (g) 2.1Aq + 4. 4.1 (Aq)
23

---

Published as a conference paper at ICLR 2025
Because of two masses in 2-dimensional space, the system has 4 degrees of freedom for configu- ration, 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 M.. Hence, the dynamics is degenerate.
Due to the degeneracy, HNNS not applicable (Greydanus et al. eracy does not come from holonomic constraint of the form that CHNNS assumed, CHNNs are also not applicable (Finzi al. 2020). In PNNs, we used 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- 4-dimensional observations were transformed, and the 8 dimensions were extracted (with the remaining 6 dimensions considered constant) as input to the energy function of HNNs.
2019). Also, since this degen-
F.3 ELECTRICAL SYSTEMS
Overview and Experimental Settings Refer to Appendix C.2 for basic characteristics of capaci- tors and inductors. Since the electric charge is analogous to displacement, capacitor corresponds to spring, while an inductor corresponds
mass in mechanical systems.
The states of a capacitor C' and inductor are the electric charge Q and magnetic flux p. respec- tively, but these are difficult to observe directly. We used the capacitor voltage Vo = Q/C and inductor current I = L as the observations, which linearly correlate with their states. PoDiNNS learned the element characteristies C and L learnable parameters and internally performed the transformations Q = CVc and ip = LI This is the same approach used by Dissipative SymOD- ENs for masses. Strietly speaking, since PoDiNNs can learn coordinate transformations, they can approximate the dynamics even without such transformations.
In a dasssseringg 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 DAE, which cannot be solved explicitly. Additionally, components that share the same effort type, such capacitors and direct 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 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 and dissipating component.
(d) FitzHugh-Nagunn Nagumo Model FitzHugh- model is model for the electrical dynamics of a biological neuron, exhibiting oscillatory behavior when the external current ./ is applied (Izhike- vich & FitzHugh 2006). The governing equations are written as:
V3
V V
W J,
3
(A1)
W = 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 Rz, an inductor L, capacitor C, direct voltage source E, and tunnel diode R1, with an external current J. The characteristics of these elements 1.0, and
are defined as 1/0.08, R₂ 0.8. diode R₁ is characterized by TR, D(vn,) V2 /3- VR, where VR, is the voltage across the diode. The membrane potential V and recovery variable W are represented by the capacitor voltage Ve and the inductor current Iz. In our formulation, resistor R₂ and diode Ri have the current and voltage as their flows, respectively. Also, because the voltage generated by the direct voltage source E is unchanged for any trials, it can be treated as a part of resistor Rg.
-0.7. The current IR, through tunnel
The initial values of V and W were sampled from the uniform distribution external current J was sampled at evenly spaced intervals within the range of 0.1 to 1.5 for each trajectory, and was kept constant within each individual trajectory to evaluate whether each model can learn the oscillatory behavior.
The
24

---

Published as a conference paper at ICLR 2025
Since the external current J/ was kept constant, the variability of the dynamics 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 behav- ior (Chua 2007). It consists of linear elements (a resistor R1, two capacitors C1 and C2. and an inductor L) and a nonlinear element R2, known as Chua's diode. With R₁ = 1, the governing
equations are:
Va, =aVV, - Va - -((RR)).
Vc, Vc, - -Vc, + It,
(A2)
īg -3Vc2:
VR. = Vc..
The parameters 1/C and B L/L. The nonlinear function f(VR,) describes the voltage-
current characteristic of Chua's 's diode:
((VR.) = m, + 0.5(mg - mIV + 1|- IVR, 1(),
(A3)
where VR, is the voltage across the diode, equal to Ve,, and me and my are parameters. our formulation, both resistor R2 and diode R2 have voltages flows.
We set a 15.6, = 28, mg = 8/7, and m4 -5/7. The initial values of Vc, Vo,, and Tr. were sampled from the uniform distribution U(-0..5,..3 0.5).
Due to its chaotic behavior, we set the time step size to At 0.01, step number of test subset to 3,000 and the number of training iterations to 1,000,000.
F,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 and cannot be coupled.
(f) DC Motor In this system, a DC motor bridges the peoctro-magneti and the rotational domains. In the electric circuit, inductor T. and resistor 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 I is attached to the DC motor. The pendulum's angle is measured from the vertical position, and its angular velocity is denoted by w. With a mass m at the rod's end and the gravitational acceleration g. the pendulum has a potential energy of -mgll cos and experiences a torque of -mgl sin f Additionally. friction d occurs at the pivot point of the pendulum.
The DC motor generates torque TDC based on the current DC through the armature. We assume linear relationship with constant K, such that TDC DCp Conversely, when the motor rotates at an angular velocity w, it produces back electromotive force given by VDC - Kw. The governing equations are:
i = w,
ml -mgl sin 0 + KT d(w),
Li - K E R(T).
We set = 2.5, 0.02 (gn(s))((())) [w// for friction d, and R(T) 0.05 for resistor The voltage source F(()) was defined as the sum of five sine waves, where the amplitude, angular velocity, and initial phase of each wave were sampled from the uniform distributions (0.1,, 0.3), (0.05m, 0.3x). and (0, 2m), respectively.
2 1.5, g 1.0, and 0.5, respectively. Also, we set d(w> =
PoDiNNs can model this multiphysics system seamlessly. In the magnetic domain, inductor L has magnetic flux & LI as its state, voltage its flow, and current as its effort. For the pendulum's motion in the rotational domain, the state is angular momentum p = mllw, the flow is torque, and the effort is angular velocity (see Table Because the properties of the armature are represented by
25

---

Published as a conference paper at ICLR 2025
inductor L and resistor R. the remaining function of the DC motor be represented by coupling between inductor L and the pendulum's 's motion p. that is, the bivector element 0i A 2 Both Bp" inductor I and mass m are considered to store kinetic energies, and this type of coupling is referred to as gyrator.
Therefore, the only notable point for the implementation is that magnetic flux yp and angular mo- mentum p. which belong different domains, can be coupled. Since there are only 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 a 0
3p öp learnable.
(g) Hydraulic Tank Consider hydraulic tank storing incompressible fluid, which belongs 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 cross-sectional area and fluid height h, giving the relationship V = Ah. Let p represent the density of the fluid per unit volume, and the gravitational acceleration. The pressure p exerted on the bottom of the tank per unit area given by p = pgh A and the total force acting on the bottom surface pgh.A.
Assuming that fluid is supplied from the bottom of the tank, the potential energy stored in the tank can be expressed as:
Uv = 5v==f L pghA
1 2 p A
Pg y²
2A
Then,
aUv OV = PgV/A p.
Consider two cylinders attached in opposite directions the bottom of the tank, with secssnncctional areas 01 and (22. Each cylinder contains piston, with masses my and m2, respectively. The forces acting on the pistons are -p01 and p02, where the positive direction towards m2. When piston mn moves by displacement 91. volume of fluid 01 flows into the tank. Similarly, when piston m2 moves by Q2. a volume a242 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.] (2009).
Each piston m; moves with velocity u, and is connected fixed wall via spring k; and damper d; for each i. An external force F = e applied to piston m22.
The equations of motion for the system can be written
V = a1 01 021'2,
mi 9j ppa k(q)) di (v1);
m2d2 = pa2 k(q2) d2(02) + F,
Each piston has kinetic energy, and each spring has potential energy. The state of each piston can be expressed as momentum, p; m, and the state of the spring can be described by the displacement for each i. Thus, the total energy of the system is given by:
H
Uilg:) V2
Pg 2A
2m
where U, denotes the potential energy stored in spring k; well the characteristics of dampers d,
The above equations of motion can be expressed using the energy function H with the following bivector B:
0
a
i
B a1
A
02
Op1
@V
Opz
Op1
0q
Op2
Oq2
A
R
^ 32
A NET,
Op
Op2
Op2
where e. and gx represent the basis vectors of the spaces Jt and FI for dampers d, and external force F. respectively.
26

---

Published as a conference paper at ICLR 2025
ground EnEE
2.5
ground truth
Sround
2.5
0.05
ODE
0.06-
ODE
neural ODE
0.00
0.25
da. SymODNN
0.00
PoDiNN
0.1
PoDINN
0.00
PoLxNN
0.00
0.0
0.00
do
(a) with external force
(b) with moving boundary
(b" with moving boundary in absolute coordinate system
relative coordinate system
relative coordinate system
touth
0.01
2.5
ncural ODE
0.02
neural ODE
0.0-
0.00
0.5- PcDINN
PoDINN
200
500
750
100K
50
75
100
(d) Futzluugh NaggNug model
motor
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 cami-transparent color. See also Fig.
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 01 and a2 of the cylinders. This type of coupling is referred to 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, p 10.0, 01 = 1.0, a2 0.3, mn 3.0, and m2 = 1.0. The spring forces were defined as kg:):)) = 0.1q, + +0.01q and the damping forces as (v;); = (0.1 - 0.04i) sgn(v;)]u;;) 11/3 for i = 1,2.
As the system tends to reach equilibrium around V 5, 91 = 10, and 92 = 6, the initial fluid volume and the springs' displacements 91 and 92 were sampled from the uniform distributions U + 0.25), (--10 0.3, -00 + 0.3), and U 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 F()) (L) was defined the sum of three sine waves, where the amplitude, angular velocity, and initial phase of each wave were sampled from the uniform distributions (0.055, 0.2), (0.1 0.3 and (0, 2x), respectively.
In Fig. 3. for clarity, we displayed the changes relative to 5, 10, and 6, rather than the actual values of V. 91, and 92-
G ADDITIONAL RESULTS AND DISCUSSIONS
G.1 ADDITIONAL, RESULLS
Additional Visualization We show example visualizations of trajectories and absolute errors in Fig. Al in addition to Fig. 3 Across all datasets, PoDiNNs demonstrate superior predictive perfor- mance compared to other methods. However, in both systems (b) and (b'), there was one failure. This was due 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.
27

---

Published as a conference paper at ICLR 2025
Identifying Coupling Patterns of System (b) We examined how the bivector B identifies the coupling patterns using system (b), that is, the mass sappsppdaapppe system with moving boundary. When two coefficients differ by factor of 1,000 or more, we considered the larger one detected coupling and the smaller one as effectively zero, indicating no coupling. The basis vector 0 By.
a to mass m,, and 5 to moving boundary corresponds to spring k;, ap.
(1 Ef) + Des + a333") in all 10 trials. Here, the basis vectors Eft E and < for dampers d1, d2, and d3 were appropriately reordered because their indices are interchangeable. 01, 022 and 03 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 2 and a A 2 were effectively zero, indicating that PoDiNNs correctly identified the 5 A Bys Sps absence of non-existent couplings.
In the relative coordinate system, we obtained the bivector B (pp) - 59)
8qg (ag 8 +
B + +a242) + Laps ōpz 1 8 ) Dyz ( (pps Upt
e 2 A (a
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. [A] This suggests that even minor errors the coefficients of the bivector elements can lead to significant prediction failures. Nonetheless, mentioned earlier, such occurrences are extremely rare.
Identifying Coupling Patterns of System (b') In the absolute coordinate system, obtained the bivector B
à A Bep + öps A a Bag + 01(0p 0 - -55) A Ef + az áp: 8p1 a A el in out of 10 trials, where coefficients with a difference of less than 1% Sy2
0 A
D + 8p2 Bi
0 0 A
were considered identical. The relationship between the masses and the springs is represented by the 8p1
EL' + 03 Sps
2. a A 0q. 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. AI
standard Poisson bivector >; Øp.
Identifying Coupling Patterns of System (c) We also examined case of system (c). We denote a
A the displacement of i-th spring in x-direction by qui- The learned bivector B was B 3mal
0gs1
Bqui
8
+
8
A
0qu2
8 +
Janyz Dquz
8_ A
8
Jmmzi
A Jqss
Orags
8 A
Jays
8 +
B
A
+ 0m,2
A 3qzs
+ 0mya 0ags 0m 8gzn 3mys
a
A a
+
a a A 8gys in all 10 trials.
0q.3 3mn 0qps
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; and masses m, 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 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 mi is coupled with springs k1 and ks. while mass m2 is coupled with springs k and kG; and spring ka is coupled with both masses my and m2, but the coupling to mi is in the opposite direction. (c) was fully identified. By reorganizing the a +
This means that the coupling pattern shown in Fig. 2
8 A a r
expression for B, we get B = Jmai
8
8
+ 8
+
Jgza
Jges
Jmy1
Tggy J4gs Jgus
), which indicates that system (c) has eight degrees of freedom. In this way. PoDiNNs can identify the constraints and handle degenerate dynamics.
a A ( 8r + 1h4.s + 89
a
0 + 0m,z2
A (1
à
+ Bass
+ 8aya
8mg2
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 na denote the assumed number of resistors whose flow is current, and ng the assumed number of resistors whose flow is voltage. We tested system (d), the FitzHuug-Naggumo model, to explore the impact of the assumed numbers ng and Rg. and summarized the results in Table A1 The correct numbers for this model are n. = 1 and n, = 1. When fewer components 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
28

---

Published as a conference paper at ICLR 2025
test subsets. Therefore, similar to the with system (b) in Table 4 can identify the number of unobservable dissptiiiii components and their flow types using the training subset.
G.2 ADDITIONAL DISCUSSIONS
Table A1: Impact # Compo- and VPT.
Non-Identifaabiityy 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 ele-
PoDINNs
Training
Test
0.00010000
0.000
nd
n9 =0
0.000
0.000
n4=2, 2.
0.00010:000
1.000
na =3 nh =0
0.000+00000
0,000
my n2
0.000:00000
0.000
ments can be normalized, as shown in Section (42 external
TA =5=5, =0
0.000+00000
0,000+0.000
forces present, they determine the scale of force, which, in tur,, uniquely determines the scale of the coupled mass. This cascading effect determines the system's 's overall scale.
By =0,0, myy
20,000
0.000=0.000
0,
0.007 10.000
0,00010.000
=4=0, m1 =3
0.007
0,00010.000
0.007
40.000
0.000
my= =
0.007 40.000
0.000
When the characteristics elements are linear, different cou- pling patterns result in the same dynamics. Additionally, weighted average of these coupling patterns can also be a valid
nx=1,
0.98S +0.002
0.640
ng= m1 =2
0.630
n,
0.985 10.003
0.681
10,062
na na m1g
0.985
0.636 108
solution.
00904
0.077
ny 3,
0.978 10.000
0.553
While non-identifiability may be problematic for system identi-
14=3,
0.984
0.557
fication, it does not impact prediction accuracy.
0= 10
0 10
Passivity and Convexity In understanding the behavior of dynamical systems, passivity can be crucial concept [Khalil] 2002). Passivity implies that the system always dissipates energy and, without external imputs, converges to stable invariant set. In PoDiNNs, general neural networks are used to model energy-dissipe dissppiissi components, which may sometimes learn energy- supgyys characteristics, like negative resistance. This flexibility is beneficial for learning special diodes used 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 2019). This ensures the system remains passive. If convexity is required in the energy function, a neural network enforcing this property can be introduced (Amos al. 2017).
PoDiNNs decompose coupled systems into individual components, along with their coupling pat- terns, and learn the specific characteristics of each component. This structure allows for the intro- duction of constraints beforehand and ensure each component to meet the necessary properties. This flexibility is another key advantage of PoDiNNs.
29


Reviews:
[{"review_id": "h1Sy6LQKut", "reviewer": "Reviewer_h9yM", "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.", "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.  \n\nCan 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.\n\nWhat have you not considered non-neural network methods like DMD... or system ID methods like mamba/S4 not been compared?", "limitations": "", "rating": 5, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 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.", "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.\n\nThere 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.\n\nFurther: \"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.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "8X4rTrGVUo", "reviewer": "Reviewer_zHG9", "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.", "questions": "Please refer to the bullet points in Weaknesses.", "limitations": "", "rating": 8, "confidence": 4, "soundness": 3, "presentation": 4, "contribution": 3, "strengths": "•\tAs 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.\n\n•\tThe 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.\n\n•\tThe results on multiple systems look promising, where various experiment scenarios and evaluation metrics are comprehensive to validate PoDiNNs’ capabilities.", "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.\n\n•\tThe 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? \n\n•\tThe 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?\n\n•\tMoreover, 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?\n\n•\tFor 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.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "LPh4RdaHXi", "reviewer": "Reviewer_RzPg", "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.\n\nUnfortunately, 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.", "questions": "-", "limitations": "", "rating": 6, "confidence": 1, "soundness": 3, "presentation": 3, "contribution": 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.", "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.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "Sw3XTuSfGU", "reviewer": "Reviewer_F183", "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.", "questions": "See above.", "limitations": "", "rating": 8, "confidence": 2, "soundness": 3, "presentation": 2, "contribution": 4, "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.\n- The experimental results are very convincing; in the chosen experimental setups the proposed framework outperforms all baseline methods consistently.\n- 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.", "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.\n- 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)?\n- 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.\n- 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.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "vt9DWoof7v", "reviewer": "Reviewer_pxB8", "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.", "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.\n\n2) 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.", "limitations": "", "rating": 6, "confidence": 2, "soundness": 3, "presentation": 2, "contribution": 3, "strengths": "1) The paper provides a solid theoretical foundation for the proposed PDNN architecture.\n\n2) The authors conduct extensive experiments to demonstrate the effectiveness of the PDNN in various applications.", "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.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Decision: Accept (Poster)