Structure-Preserving Learning and Prediction in Optimal Control of Collective Motion

1. Problem Setup and General Approach

Understanding and predicting the collective motion of interacting unmanned vehicles is becoming increasingly important for civilian and military uses. Mathematically, the problem we consider can be formulated as follows. Suppose we observe the evolution of a system of N interacting elements, each of them evolving in its own phase space. For the unmanned aerial vehicles (UAVs) or underwater drones, the phase space of the individual particle is the Lie group $\mathrm{SE}{\left(3\right)}^N$, where $\mathrm{SE}\left(3\right)$ is the group of rotations and translations in three-dimensional space. For land-based drones or drones operating on the sea surface, the phase space is also the Lie group $\mathrm{SE}{\left(2\right)}^N$. This, in general, practical applications of unmanned vehicles are occurring on some Lie group G. Suppose that we have observations of each particle of the system, and we know that the objects are connected and controlled, but we have limited information about the actual control procedure and the connection between them. The goal is to predict the collective motion of the system from data, without the need to decipher the actual information about the method of control.

In its most general formulation, without any information about the control, the prediction of motion is likely impossible. However, with some limited assumptions on the nature of control, a solution to this problem can be found. In what follows, we assume that the collective behavior comes from the optimal control procedure. The process of symmetry reduction of problems of collective behavior using optimal control was developed in Justh and Krishnaprasad [1,2]. These papers show that the control Hamiltonians are of the type allowing complete reduction to coupled Lie-Poisson systems. We are going to predict the motion in exactly that case of complete Lie-Poisson reduction as described in [1,2]. Without any knowledge of the system, the actual control Hamiltonian, depending on the directions of control and the control amplitudes, cannot be recovered completely. However, we show that instead of learning the Hamiltonian and the nature of control, we can recreate phase space trajectories and predict dynamics of the system in the future, based on the available data observations only, for an arbitrary control Hamiltonian. In order to do that, we propose a new type of Neural Networks, which we call CO-LPNets (Control Optimal Lie–Poisson Neural Networks), which are a generalization of Lie-Poisson Neural Networks [3] for the problems of optimal control of multiple interacting vehicles.

Data-Driven Discovery of Differential Equations with Structure

Recent advancements in scientific machine learning focus heavily on identifying the governing dynamics of physical systems from observational data. For autonomous systems evolving on a state space manifold, this typically involves reconstructing the underlying phase flow from a discrete sequence of observed states ${\mathbf{u}}_i=\mathbf{u}\left(t_i\right)$.

While generic, unconstrained neural networks excel at interpolating data, they frequently fail in engineering and physics applications where long-term stability hinges on respecting known conservation laws. To address this, frameworks like Physics-Informed Neural Networks (PINNs) [4,5] introduce physical boundaries and structural residuals directly into the loss function. However, because PINNs rely on optimization-based soft constraints, they regularly exhibit optimization failures, spectral bias, and severe error accumulation when applied to highly structured geometric systems, such as canonical Hamiltonian or Lie–Poisson dynamics [6]. For these systems, minor violations of the underlying geometry accumulate over time, destroying the qualitative correctness of long-term trajectories. Consequently, our work shifts the focus from soft-penalty architectures to exact, structure-preserving data-driven representations.

Continuous vs. Discrete Geometric Learning paradigms

Data-driven geometric computing generally splits into two distinct methodological paradigms based on how time is treated:

• Continuous Operator Learning: This paradigm aims to learn the continuous-time vector field or infinitesimal flow of the dynamics, modeling the governing operator as $\dot{\mathbf{u}}=\mathbf{f}\left(\mathbf{u}\right)$. In structured settings, e.g., $\mathbf{f}\left(\mathbf{u}\right)$ describing a Hamiltonian system, approaches like Hamiltonian Neural Networks (HNNs) [7] parameterize the vector field via a learned Hamiltonian function. While effective, a significant drawback of this approach is that generating actual trajectories requires pairing the learned model with an external, secondary numerical scheme (e.g., a structure-preserving integrator [8]). Any discretization error introduced by this secondary solver can still degrade the physical invariants.
• Discrete Flow-Map Learning: Rather than approximating the continuous vector field, this paradigm directly learns the discrete-time flow map ${\mathsf{\Phi }}_h:\mathbf{u}\left(t\right)\mapsto \mathbf{u}(t+h)$ over a fixed time step h. Trajectories are generated strictly through recursive map composition: ${\mathbf{u}}_k={\mathsf{\Phi }}_h\left({\mathbf{u}}_{k-1}\right)$. This architecture inherently bypasses the need for an external ODE solver. Existing literature has successfully leveraged this framework for canonical systems via SympNets [9] and extended it to non-canonical Poisson spaces [3,10,11,12,13].

The methodology developed in this paper uses the discrete data-driven computing paradigm. In other words, we learn discrete mappings and not the Hamiltonian.

By directly identifying the discrete flow map on the reduced space ${\left({\mathfrak{g}}^{*}\right)}^N$, our framework eliminates the accumulation of artificial dissipation typical of continuous-operator approximations. As detailed in the subsequent sections, embedding the Lie–Poisson geometry directly into the network architecture allows our model to preserve non-canonical invariant quantities—specifically, Casimir functions—to machine precision indefinitely, ensuring long-term bounded energy errors and stable rollout trajectories even within chaotic regimes.

2. Literature Review

2.1. Data-Based, Structure-Preserving Computing of Poisson Systems

2.1.1. Continuous Systems with Structure

Since the focus of this paper is the data-based computing of control-Hamiltonian systems, we present a short review of the recent progress in that field.

Various methodologies exist for applying physics-informed data techniques to Hamiltonian systems while respecting their intrinsic mathematical architecture. Much of the existing literature focuses on canonical Hamiltonian systems, where the state vector $\mathbf{u}$ is $2n$-dimensional, composed of coordinates $\mathbf{q}$ and momenta $\mathbf{p}$; thus, $\mathbf{u}=(\mathbf{q},\mathbf{p})$. In these cases, the evolution equation $\dot{\mathbf{u}}=\mathbf{f}\left(\mathbf{u}\right)$ is governed by a Hamiltonian function $H(\mathbf{q},\mathbf{p})$, such that the operator $\mathbf{f}$ takes the form:

$$\mathbf{f}=\mathbb{J}{\nabla }_{\mathbf{u}}H,\mathbb{J}=\left(\begin{array}{cc}0& {\mathbb{I}}_n\\ -{\mathbb{I}}_n& 0\end{array}\right),\iff \dot{\mathbf{q}}={\displaystyle \frac{\partial H}{\partial \mathbf{p}}},\dot{\mathbf{p}}=-{\displaystyle \frac{\partial H}{\partial \mathbf{q}}},$$

(1)

where ${\mathbb{I}}_n$ represents the $n\times n$ identity matrix. If we take an arbitrary smooth function $F(\mathbf{q},\mathbf{p})$, its’ time evolution along a trajectory of (1) is defined by the canonical Poisson bracket:

$${\displaystyle \frac{dF}{dt}}=\{F,H\}={\displaystyle \frac{\partial F}{\partial \mathbf{q}}}·{\displaystyle \frac{\partial H}{\partial \mathbf{p}}}-{\displaystyle \frac{\partial H}{\partial \mathbf{q}}}·{\displaystyle \frac{\partial F}{\partial \mathbf{p}}}.$$

(2)

This bracket acts as a mapping that transforms two smooth functions into another smooth function of the same variables. Mathematically, it is a bilinear, antisymmetric operator acting on differentiable function, which, satisfies the Leibniz rule (acting as a derivative), and fulfills the Jacobi identity:

$$\left\{\right\{F,G\},H\}+\left\{\right\{H,F\},G\}+\left\{\right\{G,H\},F\}=0$$

(3)

for all smooth functions $F,G,H$. Brackets that possess these four properties but differ from the form in (2) are categorized as general or non-canonical Poisson brackets. The resulting equations of motion are called Hamiltonian (usually for canonical brackets) or Poisson systems, for non-canonical Hamiltonian brackets. The non-canonical systems frequently possess a non-trivial null space, leading to the preservation of specific quantities known as Casimir invariants (or Casimirs). These invariants are structural features of the Poisson bracket itself, remaining constant regardless of the specific Hamiltonian chosen. This proposal concentrates on data-driven computational strategies for a significant class of these non-canonical Poisson systems.

One prominent research direction focuses on extracting the system’s Hamiltonian or its underlying Poisson bracket directly from observed data. This approach was pioneered for canonical systems by Hamiltonian Neural Networks (HNN) [7], which approximate $H(\mathbf{q},\mathbf{p})$ by fitting the observed state evolution through the relations in (1). This approach, respecting the structure of the equations, allows for more robust and accurate trajectory approximations than standard neural networks. This approach has been further refined to handle chaotic transitions and adaptive parameter learning [14], with theoretical foundations for the existence of such learned Hamiltonians established in [15].

An alternative strategy is found in Lagrangian Neural Networks (LNNs) [16], which model the Euler-Lagrange equations in the $(\mathbf{q},\dot{\mathbf{q}})$ coordinate-velocity space, bypassing the need for a Legendre transform to $(\mathbf{q},\mathbf{p})$ coordinates. More recently, learning vector fields for non-canonical Poisson brackets was explored in [17], though enforcing the Jacobi identity (3) remains a significant hurdle. Furthermore, recent work [18] utilized kernel ridge regression to learn Hamiltonians from noisy data in Poisson systems. While these methods can accurately reconstruct vector fields, isolating the true Hamiltonian is challenging because it is defined only up to the addition of an arbitrary function of the Casimirs, whereas the ML methods would try to find a particular Hamiltonian.

Crucially, these “continuous” methods assume that once the vector field is learned, it can be solved using standard numerical integration. However, simulating Hamiltonian systems, especially over long time scales, requires caution, as generic solvers are not structure-preserving and thus, for example, violate conservation laws for energy and momenta. To maintain physical validity in long-term simulations, one must employ variational integrators [8,19,20]. While these specialized methods preserve momentum-like quantities to machine precision, they are often significantly more computationally demanding than their non-structure-preserving counterparts.

2.1.2. Discrete Mapping Approach

We utilize a different strategy, namely, determining the phase space transformations that satisfy specific structural properties. For canonical Hamiltonian dynamics, a fundamental theorem by Poincaré asserts that the flow ${\mathsf{\varphi }}_t\left(\mathbf{u}\right)$ of the system (1)—which maps initial conditions $\mathbf{u}$ to a state at time t—is a symplectic transformation [21,22,23]. Formally, for a state $\mathbf{u}=(\mathbf{q},\mathbf{p})$, the flow ${\mathsf{\varphi }}_t\left(\mathbf{u}\right)$ must satisfy:

$${\left({\displaystyle \frac{\partial \mathsf{\varphi }}{\partial \mathbf{u}}}\right)}^T\mathbb{J}\left({\displaystyle \frac{\partial \mathsf{\varphi }}{\partial \mathbf{u}}}\right)=\mathbb{J}.$$

(4)

Rather than identifying the governing equations and subsequently integrating them, several researchers have proposed learning these symplectic mappings directly from empirical data.

Early developments in such learning of such symplectic transformations started with Symplectic Recurring Neural Networks (SRNNs) [24], which approximate transformations using discretized symplectic numerical schemes. That work determined that enforcing canonical structures significantly improves robustness when dealing with noisy datasets. Similarly, Non-Separable Symplectic Neural Networks (NSSNNs) were introduced in [25] to approximate symplectic steps for non-separable Hamiltonians, utilizing the explicit integrator proposed in [26]. From a Lagrangian perspective, [27] explored mappings in the $(\mathbf{q},\dot{\mathbf{q}})$ phase space through structure-preserving discretization, a framework we have recently extended to encompass dissipative forces and thermodynamic systems [28].

Further advancements in learning symplectic maps were achieved through three distinct methodologies [9,29]. The first, SympNets [9], constructs dynamics by composing specific types of symplectic layers. In contrast, Generating Function Neural Networks (GFNNs) [29] derive mappings via a generating function approach, allowing for explicit long-term error estimation. Additionally, [30] introduced HénonNets, which utilize Hénon mappings to learn Poincaré maps while maintaining symplectic integrity. These architectures have successfully simulated the long-term behavior of integrable systems, e.g., pendulums and planetary orbits, and chaotic regimes like the three-body problem, demonstrating the potential of direct transformation learning.

The SympNet framework was adapted for non-canonical Poisson systems in [10] under the name Poisson Neural Networks (PNNs). By invoking the Lie-Darboux theorem, PNNs transform non-canonical coordinates into local canonical ones, assuming a constant-rank Poisson structure. However, this approach faces two major hurdles: the Lie-Darboux transformation is usually (and for most cases) local rather than global, and the model fails to preserve Casimirs exactly. However, Casimir preservation is essential for the reliability of long-term ensemble simulations [31], so we will be interested in methods that preserve Casimirs with machine precision.

An alternative route suggested in [12] involves solving approximate Hamilton-Jacobi equations for symmetry-reduced systems. This allows for the construction of transformations that exactly conserve momentum maps and Casimirs; however, solving the Hamilton-Jacobi equations for general, complex systems remains computationally and analytically prohibitive. More recently, [32] introduced LocSympNets and their symmetric counterparts, SymLocSympNets, which rely on the composition of volume-preserving mappings. These models have shown high accuracy in long-term simulations of rigid body dynamics, linear advection, and particles in magnetic fields. While these methods are highly effective for Poisson-like problems, their broader application has been limited because the mathematical completeness of their proposed activation matrices remains unproven.

2.2. Optimal Control via Neural Networks

Because of the importance of the problem of optimal control in engineering, and flexibility of neural networks, optimal control via neural networks remains an active area of research. For example, [33] introduces the equations of motion in Pontryagin’s principle as a soft constraint, in the framework similar to what that later becomes known as PINNs. The authors of [34] used the input convex neural networks for control tractability. Applications of neural networks to problems such as robotics manipulation [35] and landing problems in aeronautics [36] were considered. For a recent survey of applications of neural networks in optimal control problems we refer the reader to a recent review [37]. While the progress achieved has been substantial, these methods have not used the advances in structure-preserving networks to fully utilize the mathematical structure inherent to Pontryagin’s principle. However, preservation of the structure of the equations in neural networks is essential for reproducing accurate long-term dynamics [3,9], which is exactly the point of our paper.

2.3. Data-Based Predictions for Lie-Poisson Systems

As discussed in more details in Section 4, the problem of optimal control for coupled autonomous systems naturally adopts the Lie-Poisson structure characteristic of many physical systems. Drawing on the work of [1,2], symmetry reduction of controlled Hamiltonian systems results in a complete Lie-Poisson reduction to the direct product of individual Lie algebras, representing the reduced phase space for each vehicle. The symmetry reduction for multi-particle systems is often incomplete, requiring the tracking of relative Lie group elements in addition to motion on the Lie algebra [38]. However, such complexity is bypassed in Lie-Poisson optimal control due to the specific functional form of the Hamiltonian [1,2].

In Lie-Poisson systems, the Poisson bracket is known a priori, as its’ explicit form is derived directly from the underlying Lie algebra commutator. Selecting the specific Hamiltonian is a subsequent task dictated by the physics or control objectives of the problem. To facilitate data-driven computing within this framework, [3,13] developed methods to construct activation maps directly from these brackets. These methods predict dynamics using maps generated by explicit solutions of specific Lie-Poisson systems. As demonstrated in [3], this approach is universally applicable to finite-dimensional Lie-Poisson systems and any Poisson system where explicit integration of the governing equations for suitable transformations is feasible.

The computational advantage of using explicit integration in Lie-Poisson equations to accelerate time-stepping was initially identified in [39], though its application was restricted to specific forms of the Hamiltonian. The methodology we develop in this article, by contrast, accommodates arbitrary Hamiltonians. In fact, the Hamiltonian does not need to be specified to construct the Neural Network; we only require knowledge of the underlying symmetry and the corresponding Lie-Poisson bracket. This is a robust assumption, as symmetries and brackets are fundamental to the problem’s geometric formulation and are known exactly. In contrast, the Hamiltonian is often a modeling choice that may only be known approximately [40].

In Lie-Poisson systems, the configuration manifold is defined on a Lie group G. For a group-invariant system, the dynamics can be reduced from the canonical bracket on the cotangent bundle $T^{*}G$ to the dual of the Lie algebra, ${\mathfrak{g}}^{*}$. This is achieved using the reduced momentum $\mathit{\mu }={\mathbf{g}}^{-1}\mathbf{p}$ (left-invariant) or $\mathit{\mu }=\mathbf{p}{\mathbf{g}}^{-1}$ (right-invariant). For instance, the rotation of a satellite is described by $\mathrm{SO}\left(3\right)$, where $\mathit{\mu }\in \mathfrak{so}{\left(3\right)}^{*}$ can be represented as a vector. The resulting equations of motion and the Poisson bracket for $\mathit{\mu }$ are:

$$\dot{\mathit{\mu }}=\{\mathit{\mu },H\},\iff \dot{\mathit{\mu }}=\mathsf{\Lambda }\left(\mathit{\mu }\right)\nabla H\left(\mathit{\mu }\right),\{F,G\}=\pm 〈\mathit{\mu },\left[{\displaystyle \frac{\partial F}{\partial \mathit{\mu }}},{\displaystyle \frac{\partial G}{\partial \mathit{\mu }}}\right]〉$$

(5)

where $[·,·]$ denotes the Lie algebra commutator. Here, $\mathsf{\Lambda }\left(\mathit{\mu }\right)$ is a matrix-valued function linear in $\mathit{\mu }$, with a structure determined by the Lie group. These brackets often admit Casimir invariants—functions C such that $\{C,H\}=0$ for any H—which are central to our approach.

Recent works [3,13] utilized the fact that the phase flow of any Hamiltonian constitutes a Poisson transformation that exactly preserves the bracket [22]. These works consider the flows generated by the “test” Hamiltonians of the form $F(\mathit{\mu }·\mathbf{a})$, where F is an unknown function and $\mathbf{a}$ is a parameter vector. For these specific forms, the Lie–Poisson equations (5) simplify significantly and can be solved analytically, providing the basis for our Lie–Poisson neural network layers.

By taking a composition of these Poisson transformations, the neural network in [3] fits the parameters of F and $\mathbf{a}$ to approximate the available phase space data. This procedure allows the network to learn the flow with high precision across the entire phase space using very modest amount of data and relatively few parameters. Critically, this architecture ensures that Casimir invariants are conserved to machine precision regardless of the training data.

While efficient, the mathematical completeness of the transformations proposed in [3] remained unproven. In this paper, we demonstrate that while the LPNets approach can be made complete through certain generalizations, an implementation of these methods is impractical. Instead, we propose a more functional generalization of Poisson transformations that enables efficient and practical learning of the evolution across the entire phase space.

3. Contributions of This Paper and Novelty of Results

3.1. Contributions

In this paper, we consider systems arising from the optimal control of the collective motion of a group of particles. The methods we develop allow us to derive explicit Poisson transformations for arbitrary Lie-Poisson dynamics, for an arbitrary number of particles, on arbitrary Lie groups, using a generalization of LPNets. We explicitly develop these techniques for interacting identical particles, with each particle evolving on the groups $\mathrm{SO}\left(3\right)$ of rotations and $\mathrm{SE}\left(3\right)$ of rotations and translations. The first case is relevant to the collective orientation of objects like satellites, or perhaps other objects like antennas, that rotate about their centers. The second case is relevant to the motion of unmanned vehicles in air or water, allowing three-dimensional rotation and translation. For each case, we assume that we know the symmetry group, which leads us to the knowledge of the Poisson bracket for each case. Our goal is to predict the collective motion of these particles from observations. Our algorithm proceeds as follows.

• Step 1 We derive the Poisson transformations approximating the mappings in phase space as the flows from the test Hamiltonians, which are computed in the explicit form, as shown in Section 5. For a given Poisson bracket, a control Hamiltonian generates a phase flow that generates a Poisson map of the phase space into itself. That Poisson map is parameterized by a neural network depending on the initial conditions for each step. To learn the dynamics in phase space, the Poisson map will be approximated as a superposition of the Poisson maps coming from the phase flow of a test Hamiltonian. Since every map approximating the full motion is Poisson, it preserves the Poisson bracket and hence the Casimir invariants (Casimirs) to machine precision; hence, the superposition of the mappings of the phase flow generated by the test Hamiltonian will also preserve the Casimirs with machine precision.
• Step 2 Given a control Hamiltonian coming from the physical system, we generate the phase flow using a high-precision numerical method. That simulation generates data in phase space that are considered the ground truth data. The test Hamiltonians have several parameters; the parameters are optimized in such a way that the superposition of mappings generated by the test Hamiltonians approximates the flow, on average, for all available data.
• Step 3 After the optimization, we compare numerical solutions, obtained via solving differential equations coming from the true control Hamiltonian (the ground truth), with predictions based on the approximation of the phase flow by a sequence of Poisson transformations as flows from test Hamiltonians.

As an example of our technique, we study the flow for two types of control Hamiltonians, which we call ’Dictatorship’ and ’Democracy’. By ’Dictatorship’, we denote the case when only one particle is interacting with all other particles, and there are no other connections. The ’Democracy’ means that all the evolving particles are connected to each other. Following [1,2], we obtained Lie–Poisson reduced dynamics equations for $\mathrm{SO}\left(3\right)$ and $\mathrm{SE}\left(3\right)$ groups in both cases. Interestingly, analytic results for some properties are possible for both of these cases, which we present below in Section 4, extending the results in [1,2]. For reference, we also present the results of Lie–Poisson reduced dynamics equations in the case of one particle, for general reference for these groups, as they have not been presented before in the literature, making the structure of equations for several particles clearer.

3.2. Novelty of Results

This work achieves the following progress compared to the previous literature:

• We present a novel network architecture of Poisson transformations, specifically formulated for the problems of optimal control of collective motion (although in principle it is also applicable to regular Lie-Poisson systems);
• Our method preserves all Casimir invariants with machine precision for all times;
• The network is highly efficient, allowing to accurately learn complete phase space dynamics of (relatively) high dimensional systems with only a modest number of data and network parameters;
• The architecture of our network allows to prove completeness, as opposed to earlier results using Poisson transformations for machine learning [3];
• The network does not need information about the control directions, the number of control dimensions, or anything else beyond the fact that the control Hamiltonian exists and satisfies symmetry reduction principles [1,41];
• We demonstrate the system’s robustness with respect to noise, and report on improvement of prediction quality with small noise in data.

Note that in our method, we never compute either the Hamiltonian H, its gradients, or the equations of motion. Instead, we compute only the composition of Poisson transformations, reproducing the dynamics in phase space of some Poisson system, coming from an unknown control Hamiltonian.

Plan of the Paper

The rest of this paper is organized as follows. Section 4 presents the theoretical background material necessary for this article. Section 5 gives a detailed explanation of how to construct the Poisson transformations out of the test Hamiltonians, which will be the main tool for data-based computing for optimal control problems, including the proof of completeness of our network. Section 6 demonstrates application of CO-LPNets towards two Lie–Poisson systems, one of which evolves on $\mathrm{SO}{\left(3\right)}^N$ Lie group, which is described in Subsection 6.1, and the other one on the group $\mathrm{SE}{\left(3\right)}^N$ Lie group, which is described in Subsection 6.2 (here N is the number of particles). Section 6.3 studies the dependence of the system on the noise in the data. Section 7 provides the Conclusions and directions for future work. In the Appendix, Section A.1 provides detailed derivations of equations of motion; Section B provides Lie–Poisson reduced dynamics equations for the case of “single particle” for reference, in order to present the whole theory more clearly. Section C provides the exact values of the derivatives of the Hamiltonians used for data generations and ground truth computations. Section D of the Appendix provides partial derivatives of the transformations.

4. A Short Review of Symmetry-Reduced Motion of the Particles Under Optimal Control

4.1. The Lie–Poisson Framework for Multi-Agent Optimal Control

To establish a self-contained foundation for our data-driven learning methodology, this section outlines the geometric framework of symmetry-reduced optimal control for interacting particle systems. This short derivation allows us to define mathematical structure, defines main variables and defines how the ground-truth trajectories are generated. While we introduce all essential kinematic relations, graph topologies, and the final reduced dynamical equations here, the exhaustive, step-by-step algebraic derivations are deferred to Appendix A.1 to preserve narrative clarity. Our formulation builds directly upon the geometric control foundations established in [1,2].

System Configuration and Network Topology

We consider an ensemble of N identical particles, where the configuration of each individual agent is represented by an n-dimensional Lie group G with an associated Lie algebra $\mathfrak{g}$. The total configuration space of the collective system is the direct product group $G^N$. Let $\{X_1,\dots ,X_n\}$ be an orthonormal basis for $\mathfrak{g}$ under the standard trace inner product $\langle \langle \xi ,\eta \rangle \rangle =Tr\left({\xi }^{\mathsf{T}}\eta \right)$.

The interaction topology among the N agents is modeled by a connected, undirected graph $\mathcal{G}=(V,E)$ without self-loops, where vertices represent agents and edges denote active communication links. The network architecture is characterized by three standard matrices [42]:

• Adjacency Matrix (A): $a_{ij}=1$ if a connection exists between agents i and j, and 0 otherwise.
• Degree Matrix (D): A diagonal matrix $D=diag(d\left(v_1\right),\dots ,d\left(v_N\right))$, where $d\left(v_i\right)$ is the degree of vertex $v_i$.
• Graph Laplacian (B): Defined as $B=D-A$. By construction, B is symmetric and positive semi-definite, satisfying $B{\mathbf{1}}_N=\mathbf{0}$, meaning the row and column sums are identically zero.

Controlled Dynamics and Coupled Cost Functionals

The left-invariant kinematic evolution of the k-th agent on G is governed by:

$${\dot{g}}_k=g_k{\xi }_k,{\xi }_k\in \mathfrak{g},k=1,\dots ,N$$

(6)

We assume the system is underactuated ($m<n$) with an uncontrolled drift component $X_q$ along a specific algebra direction ($m+1\le q\le n$). The control velocities ${\xi }_k$ are affine in the control inputs $u_k\in {\mathbb{R}}^m$:

$${\xi }_k=X_q+{\displaystyle \sum _{i=1}^m}{u}_{ki}{X}_i$$

(7)

The control objective is to minimize a G-invariant cost functional depending exclusively on the control inputs:

$$\underset{u}{min}{\int }_0^{T}L({\xi }_1,\dots ,{\xi }_N)dt$$

(8)

where the physical coupling between interacting particles is explicitly encoded via a coupling strength parameter $\chi \ge 0$ [1]:

$$L({\xi }_1,\dots ,{\xi }_N)={\displaystyle \frac{1}{2}}{\displaystyle \sum _{k=1}^N}{\parallel {\xi }_k\parallel }^2+\chi {\displaystyle \sum _{k=1}^N}〈〈{\xi }_k,{\displaystyle \sum _{j=1}^N}{b}_{kj}{\xi }_j〉〉$$

(9)

Symmetry Reduction and Reduced Hamiltonian Dynamics

In general mechanical systems exhibiting G-symmetry on a product manifold $G^N$, a complete reduction to the dual algebra ${\left({\mathfrak{g}}^{*}\right)}^N$ is impossible without tracking relative configuration variables ${\alpha }_{ij}=g_i^{-1}{g}_j$, which introduces severe analytical complexity [38]. However, owing to the specific algebraic structure of the control Hamiltonian under Pontryagin’s Maximum Principle [41,43], these optimal control formulations admit a seamless, complete symmetry reduction to a pure Lie–Poisson system.

Let ${\mu }_k={\sum }_{i=1}^n{\mu }_{ki}{X}_i^b\in {\mathfrak{g}}^{*}$ denote the symmetry-reduced co-state (momentum) of the k-th particle, and let ${\tilde{\mu }}_k={({\mu }_{k1},\dots ,{\mu }_{km})}^{\top }\in {\mathbb{R}}^m$ isolate its directly controlled components. By applying the first-order optimality condition $\partial H/\partial u=0$, we find an explicit linear mapping between the optimal control inputs and the controlled momenta components (see Appendix A.1 for the complete inversion step):

$$\mathbf{u}=\mathsf{\Psi }\tilde{\mu },\mathrm{where}\mathsf{\Psi }={\left(({\mathbb{I}}_N+2\chi B)\otimes {\mathbb{I}}_m\right)}^{-1}$$

(10)

Here, $\mathbf{u}={[u_1^{\top },\dots ,u_N^{\top }]}^{\top }$ and $\tilde{\mu }={[{\tilde{\mu }}_1^{\top },\dots ,{\tilde{\mu }}_N^{\top }]}^{\top }$. Substituting these optimal controls back into the pre-Hamiltonian yields the explicitly reduced control Hamiltonian $h:{\left({\mathfrak{g}}^{*}\right)}^N\to \mathbb{R}$:

$$h\left(\breve{\mu }\right)={\displaystyle \sum _{k=1}^N}{\mu }_{kq}+{\displaystyle \frac{1}{2}}{\tilde{\mu }}^{\top }\mathsf{\Psi }\tilde{\mu }$$

(11)

where $\breve{\mu }={[{\mu }_1^{\top },\dots ,{\mu }_N^{\top }]}^{\top }$ represents the aggregate momentum vector of length $Nn$.

Consequently, the global collective trajectories evolve entirely on ${\left({\mathfrak{g}}^{*}\right)}^N$ according to the Lie–Poisson equations:

$$\dot{\breve{\mu }}=\mathsf{\Lambda }\left(\breve{\mu }\right)\nabla h,\mathsf{\Lambda }\left(\breve{\mu }\right)=-{\displaystyle \frac{1}{\sqrt{2}}}diag\left({\widehat{\mu }}_1,{\widehat{\mu }}_2,\dots ,{\widehat{\mu }}_N\right)$$

(12)

where each structural Poisson block ${\widehat{\mu }}_k$ is an $n\times n$ antisymmetric matrix defined via the structure constants ${\mathsf{\Gamma }}_{ij}^s$ of the underlying Lie algebra ($[X_i,X_j]={\sum }_s{\mathsf{\Gamma }}_{ij}^s{X}_s$):

$${\left[{\widehat{\mu }}_k\right]}_{ij}=-{\displaystyle \sum _{s=1}^n}{\mu }_{ks}{\mathsf{\Gamma }}_{ij}^s$$

(13)

Remark 1 

(Geometric Properties). Through left-invariance, equations (12) can be cast in coordinate-free notation via the co-adjoint action: ${\dot{\mu }}_k=-{ad}_{\partial h/\partial {\mu }_k}^{*}{\mu }_k$ [44,45]. Crucially, the structure of the Poisson tensor dictates that the origin $\breve{\mu }=\mathbf{0}$ remains a stationary invariant point for any choice of control Hamiltonian h. This strict structural constraint is heavily leveraged in our data-driven discovery phase.

4.2. Algorithm for Producing the Data in Phase Space

The algorithms presented in this Section are used for producing data for a particular realization of the phase space, which are taken to be N interacting particles in either $\mathrm{SO}\left(3\right)$ or $\mathrm{SE}\left(3\right)$.

We solve the Lie–Poisson reduced dynamics equations (12) for two types of control Hamiltonians described above, which we call ’Dictatorship’ and ’Democracy’. We generate $N_t=40$ short trajectories for $\mathrm{SO}\left(3\right)$ and $N_t=80$ trajectories for $\mathrm{SE}\left(3\right)$ of the length $N_p=51$ points each. We choose $\chi =0.5$ for all simulations. The initial conditions for these trajectories are generated randomly with a uniform distribution in the box ${\mu }_k\in {[-1;1]}^3$ for every particle $k=1,\dots ,N$. The numerical solutions are obtained using a high-precision Lie-Poisson integrator [46], based on an implicit second-order midpoint rule, with second-order accuracy in the time step. The integrator conserves Casimir invariants with machine precision and Hamiltonian with the accuracy of about ${10}^{-14}$. The same integrator is used to construct ground truth for comparison with the reconstructed solutions obtained by the neural network. Such an accuracy may seem excessive as it far exceeds any reasonable accuracy of experimental data. However, it is important to have a structure-preserving integrator to avoid possible drift in solutions for the comparison to the reconstructed solutions over long times.

5. Data-Based Computing of the Coupled Lie-Poisson Control Systems

5.1. General Considerations and Completeness Result

The main goal of this article is to approximate and predict further evolution of the Lie–Poisson system. In what follows, we assume that we do not know the Hamiltonian itself, and thus we do not assume the knowledge of the Lagrangian (9), the form of the controls (7), or a particular form of the Hamiltonian (11). The Lagrangian and controls are assumed to be arbitrary, as long as the non-degeneracy of the Hessian for the pre-Hamiltonian in (see eq.(A5) in the Appendix) is satisfied. Our goal is to predict the evolution of the Lie-Poisson control system based on available data observations only.

To generate the data set, we consider a set of points corresponding to the start and end points of short trajectories obtained from the exact equations of motion as described in Section 4.2. This procedure of data generation gives a set of data points in the phase space connected by the phase space transformation over a short time interval. Our goal is to learn the phase space transformation connecting the data points in some optimal sense that we quantify below.

Remark 2 

(On data observability and applicability of our theory). The ground truth data obtained in Section 4.2 used for learning the system involves information about the co-states of the system μ. Such information is not readily available from the data collected by an external observer. The data can, however, be obtained from the internal control methods of the system itself. The methods for data-based predictions we derive here can then be used to compare the actual evolution of the system to the predicted evolution in case there are inaccuracies in the control application or the analytic expression for the control Lagrangian/Hamiltonian. In these cases, the data-based prediction that we develop here can be a valuable tool for fine-tuning the system and predicting its evolution.

As we mentioned in Section 1, in this work, we do not seek the approximation of the actual control Hamiltonian from data. Instead, we are focusing on the approximation of the mappings of the phase space in the $\mu$ variables between the data points in time using a sequence of structure-preserving Poisson maps. For this purpose, we consider Hamiltonians h for which we can solve the equations of motion analytically. We call such Hamiltonians test Hamiltonians. These test Hamiltonians consist of arbitrary functions of projections of the momenta $\mu$ on a certain vector. We construct Poisson transformations as the flows of these test Hamiltonians, and apply these Poisson transformations sequentially to approximate the flow generated by the real control Hamiltonian. We start our considerations with the following

Lemma 1. 

Suppose $h\left(\mu \right)=f\left(\xi \right)$, where $\xi =〈\alpha ,\mu 〉+\beta$, α is a vector belonging to the Lie algebra ${\mathfrak{g}}^N$, dual to $\breve{\mu }\in {\left({\mathfrak{g}}^{*}\right)}^N$, and β is a constant. Then the Lie–Poisson reduced dynamics (A16) reduces to equations that are linear in μ.

Proof. 

By definition, any Hamiltonian h is a constant, thus, $\xi =〈\alpha ,\mu 〉+\beta =const$ as well. Moreover, since $\nabla h=a{f}^{\prime }\left(\xi \right)=const$. According to (12), $\mathsf{\Lambda }\left(\breve{\mu }\right)$ is a block-diagonal matrix. Then, (12) suggests that the Lie-Poisson reduced dynamics equations for i-th component of k-th particle will look as follows:

$$\begin{array}{c}{\dot{\mu }}_{ki}=-{\displaystyle \sum _{j,s=1}^n}{\mu }_{ks}{\mathsf{\Gamma }}_{ij}^s{\displaystyle \frac{\partial h}{\partial {\mu }_{kj}}}\end{array}$$

that is, for every $k=1,\dots ,N;i=1\dots n$, we have:

$$\begin{array}{c}{\dot{\mu }}_{ki}=-{\displaystyle \sum _{j,s=1}^n}{\mu }_{ks}{\mathsf{\Gamma }}_{ij}^s{a}_{kj}{f}^{\prime }\left(\xi \right).\end{array}$$

(14)

We know that $f^{\prime }\left(\xi \right)=const$ on solutions of (14). Hence, the term multiplying the function $f^{\prime }\left(\xi \right)$ is constant on solutions in (14), and thus the equation (14) can be solved as a linear system of ODEs in the unknown $\mu$. □

Connection to Splitting Methods in Numerical Analysis

Any Hamiltonian $h=f\left(\xi \right)$ described in Lemma 1 generates a phase flow that is a Poisson transformation which can be computed analytically. Suppose the total Hamiltonian h can be separated into a finite number of terms described in Lemma 1, i.e.$h={\sum }_{i=1}^M{f}_i(〈{\alpha }_i,\mu 〉+\beta )$. Then, the idea of splitting methods [47,48] of the traditional numerical analysis is to divide each time step $\Delta t$ into M sub-steps ${\tau }_j$, with ${\tau }_1+{\tau }_2+\dots {\tau }_M=\Delta t$, with each sub-step j yielding the phase flow ${\mathsf{\Phi }}_{{\tau }_j}^j$ according to the solution of (14). The total phase flow after time ${\Delta }_t$ can be approximated as a composition of the phase flows by each sub-Hamiltonian $h_i=f_i(〈{\alpha }_i,\mu 〉+\beta )$:

$${\mathsf{\Phi }}_{\Delta t}={\mathsf{\Phi }}_{{\tau }_M}^M\circ {\mathsf{\Phi }}_{{\tau }_{M-1}}^{M-1}\circ \dots \circ {\mathsf{\Phi }}_{{\tau }_1}^1,{\displaystyle \sum _{j=1}^M}{\tau }_j=\Delta t.$$

(15)

Completeness of transformations defined in Lemma 1

To prove the completeness of transformations derived in Lemma 1, we use a corollary of the Universal approximation theorem for the function and derivatives from [49]:

Theorem 1. 

Take any smooth, monotonically increasing function $\sigma \left(\xi \right)$. Consider any $ϵ>0$ and any compact set U of the phase space. Then, any smooth Hamiltonian $h\left(\mu \right)$ can be approximated in U up to the accuracy ϵ in the second Sobolev norm $W_p^2\left(U\right)$, $1\le p<\infty$ (the function and its first and second derivatives belong to $L_p\left(U\right)$):

$$\begin{array}{c}{∥h-{\displaystyle \sum _{j=1}^M}\sigma \left(〈{\alpha }_j,\breve{\mu }〉+{\beta }_j\right)∥}_{W_p^2\left(U\right)}<ϵ.\end{array}$$

(16)

In fact, the universal approximation theorem in [49] is stronger: it states that one can approximate all derivatives of a function up to order m on a compact domain U with multi-layer networks, not just a single-layer network as in (16). However, multi-layer networks do not satisfy the conditions of Lemma 1 and the resulting Hamiltonians for each sub-step will not be integrable in explicit form. Also, we do not need to approximate the derivatives of the Hamiltonian of order higher than two. Using Theorem (16), we formulate the following

Lemma 2 

(On the approximation of phase flow). Suppose a phase flow of the Poisson system ${\displaystyle \frac{d\breve{\mu }}{dt}}=\mathsf{\Lambda }\left(\breve{\mu }\right)\nabla h$ (A16) is defined by a smooth Hamiltonian $h\left(\breve{\mu }\right)$. Choose a compact set U in the phase space, and some $T>0$. Then, take any $ϵ>0$, and choose the set $U_{ϵ}:\{\mathbf{x}\in U:\mathrm{dist}(\mathbf{x},\partial U>ϵ\}$1 Then, there is an integer $M>0$, and constant vectors ${\alpha }_j$ in ${\mathfrak{g}}^N$, and ${\beta }_j\in \mathbb{R}$ such that the flow is approximated by the phase flow of the Hamiltonian

$$h_{\mathrm{approx}}={\displaystyle \sum _{j=1}^M}\sigma (〈{\alpha }_j,\breve{\mu }〉+{\beta }_j)$$

(17)

up to the accuracy ϵ, for all $0<t<T$ as long as the solution generated by the phase flow of the Hamiltonian h does not leave $U_{ϵ}$.

Proof. 

Let us choose arbitrary $ϵ>0$, and some ${ϵ}_1=ϵ/D$, where $D>0$ is to be determined. We also choose $p>d$, where $d=\dim \breve{\mu }=Nn$. Let us take M and $({\alpha }_j,{\beta }_j)$ such that the true Hamiltonian h and $h_{\mathrm{approx}}$ are ${ϵ}_1$-close in $W_p^2$ norm. (16) is satisfied for $ϵ={ϵ}_1$. Then, the Poisson vector fields defined by ${\displaystyle \frac{d\breve{\mu }}{dt}}=-\mathsf{\Lambda }\left(\breve{\mu }\right)\nabla h$ for $h=h\left(\breve{\mu }\right)$ and $h=h_{\mathrm{appr}}\left(\breve{\mu }\right)$ are ${ϵ}_1$-close in $W_p^2$, and, since $p>d$, by Sobolev’s embedding theorem, they are $C^1$-close uniformly in U.

Suppose L is the Lipshitz constant of the Poisson vector field $\mathbf{f}\left(\breve{\mu }\right)=\mathsf{\Lambda }\left(\breve{\mu }\right)\nabla h$. We define ${\mathbf{f}}_1\left(\breve{\mu }\right)=\mathsf{\Lambda }\left(\breve{\mu }\right)\nabla h_{\mathrm{approx}}$ to be the vector field generated by the approximate Hamiltonian. Then, the solutions $\mathbf{x}\left(t\right)$ and ${\mathbf{x}}_1\left(t\right)$ satisfy the integral equations

$$\mathbf{x}=\mathbf{x}\left(0\right)+{\int }_0^t\mathbf{f}\left(\mathbf{x}\left(s\right)\right)ds,{\mathbf{x}}_1={\mathbf{x}}_1\left(0\right)+{\int }_0^t{\mathbf{f}}_1\left(\mathbf{x}\left(s\right)\right)ds.$$

(18)

Assuming that the initial conditions are the same, ${\mathbf{x}}_1\left(0\right)=\mathbf{x}\left(0\right)$, from (18) we obtain the estimate for the error $e\left(t\right)=\left|\mathbf{x}\left(t\right)-{\mathbf{x}}_1\left(t\right)\right|$, for all t such that $\mathbf{x}\left(t\right)$ and ${\mathbf{x}}_1\left(t\right)$ are in $U_{ϵ}$:

$$\begin{array}{cc}\hfill e\left(t\right)& =\left|\mathbf{x}\left(t\right)-{\mathbf{x}}_1\left(t\right)\right|=\left|{\int }_0^t\mathbf{f}\left(\mathbf{x}\left(s\right)\right)-{\mathbf{f}}_1\left({\mathbf{x}}_1\left(s\right)\right)ds\right|\hfill \\ \hfill & =\left|{\int }_0^t\mathbf{f}\left(\mathbf{x}\left(s\right)\right)-{\mathbf{f}}_1\left(\mathbf{x}\left(s\right)\right)+{\mathbf{f}}_1\left(\mathbf{x}\left(s\right)\right){\mathbf{f}}_1\left({\mathbf{x}}_1\left(s\right)\right)ds\right|\le {ϵ}_1+{\int }_0^{t}Le\left(s\right)ds.\hfill \end{array}$$

(19)

By the integral Gronwall’s inequality [50], equation (19) leads to the estimate for $e\left(t\right)$ as

$$e\left(t\right)\le {\displaystyle \frac{{ϵ}_1}{L}}\left(e^{Lt}-1\right).$$

(20)

Thus, for any $0<t<T$, such that the solution does not leave $U_{ϵ}$, we can take

$${ϵ}_1=ϵ{\displaystyle \frac{L}{e^{LT}-1}},$$

(21)

which guarantees that $e\left(t\right)=|\mathbf{x}\left(t\right)-{\mathbf{x}}_1\left(t\right)|<ϵ$, for $0<t<T$ as long as $\mathbf{x}\left(t\right)$, ${\mathbf{x}}_1\left(t\right)$ staying in U. If the solution leaves U at $T_1<T$, then (21) still works and $e\left(t\right)<ϵ$, but only on the interval $0<t<T_1$. □

Practical Applications of Completeness Results

In principle, Lemma 2, yielding the completeness result, allows us to construct a splitting method based on the approximation of the Hamiltonian by the separable ansatz (17). One would then take M sub-steps on each time step and approximate the flow uniformly in the compact domain U. Thus, Lemma 2 serves as an existence proof of a broad class of splitting methods capable to approximate the phase flow with an arbitrary precision.

In reality, such a straightforward application of the method is not practical. Since only the approximation by single-layer networks is allowed for a separable ansatz (17), the number of terms M needed for accurate approximations of the dynamics can easily become astronomical. Estimates for deep ReLU networks [51], for example, yield the estimates for the function and first derivative with accuracy $ϵ$$M\sim {ϵ}^{-d}$, where d is the number of dimensions, although one-layer networks tend to provide even worse approximations. Taking a modest value $ϵ=0.1$ and the number of dimensions $d=18$ which we will use for three interacting particles in $\mathrm{SE}\left(3\right)$ in Section 6.2, we arrive to the astonishing number of $M\sim {10}^{18}$ substeps per time step. Clearly, writing a machine learning scheme with that number of parameters and having that number of substeps is impractical to say the least. Even if that number is reduced substantially, there will never be enough data in practical applications to construct such a machine learning scheme, which in addition will be extremely slow and memory-demanding. In what follows, we develop a numerical scheme that has the number of substeps growing quadratically with d, the dimension of the phase space, which is highly efficient. Our scheme, exploiting the block-diagonal structure of the Lie-Poisson matrix to radically reduce the parameter count. With these modifications, our scheme will be able to address a high number of dimensions with a relatively small number of data points (only a few thousand), to learn the complete dynamics of the phase space.

Remark 3 

(On polynomial Hamiltonians). A theorem on polynomial expansion [47], Chapter 3, says that if the Hamiltonian $H\left(\breve{\mu }\right)$ is a polynomial, there exist an expansion $h={\sum }_{j=1}^M{f}_j\left(〈{\alpha }_j,\breve{\mu }〉\right)$ where each $f_j$ is a polynomial of one variable. The Hamiltonian given by (11) is quadratic, but that only happened because, following [1,2] we made a very particular choice of the Lagrangian and controls. In general, there is no guarantee that the Hamiltonian is polynomial in components of $\breve{\mu }$. Our goal will be to describe the motion without any restriction on the type of Hamiltonians considered. One may be tempted to use a Taylor expansion of the Hamiltonian $h\left(\breve{\mu }\right)$ to approximate it by a polynomial. Unfortunately, such expansion will fail, for example, if we multiply the quadratic Hamiltonian defined by (11) by a prefactor $(1+\parallel \breve{\mu }{{\parallel }^2)}^{-1}$ and try to predict the motion of solutions in the domain $\parallel \breve{\mu }\parallel \le 2$, as the Taylor expansion will fail beyond $\parallel \breve{\mu }\parallel \ge 1$. We thus proceed without making any assumptions on how well the Hamiltonian can be approximated by polynomials.

5.2. CO-LPNets: General Consideration and Design

The Design of Poisson Transformations

We construct the solutions of the equations (14) with the test Hamiltonian; these solutions will generate the transformations that we will use to learn the dynamics of the phase space. For N interacting particles, each evolving on a Lie group of dimension n, we choose the test Hamiltonians $h_k=a_k{\mu }_{j_k}$, where $j_k$ is the chosen component for the k-th case. This test Hamiltonian generates the phase flow according to the equation

$$h_k=w_k{\mu }_{j_k},{\displaystyle \frac{d\breve{\mu }}{dt}}=\mathsf{\Lambda }\left(\breve{\mu }\right)w_k{\mathbf{e}}_{j_k},\breve{\mu }(t;{\breve{\mu }}_0)={\mathbf{P}}_k(t,{\breve{\mu }}_0;w_k),$$

(22)

where ${\mathbf{e}}_{j_k}$ is the basis vector with 1 at the $j_k$ coordinate and 0 otherwise. The solution of this linear equation depends on the structure of the group, as well as the scalar $a_k$ setting the time scale for the particular equation. Usually, we take K to be a multiple of $Nn$, the total number of particles, and take $j_k=kmod\left(Nn\right)$, so the index $j_k$ runs repeatedly through $1,2,\dots Nn$ a certain number of times.

The crucial observation is that the transformations ${\mathbf{P}}_k$ defined by (22) remain Poisson if we make $w_k$ a function of the initial conditions ${\breve{\mu }}_0$. It is important to note that we set the value of parameter $w_k=w_k\left({\breve{\mu }}_0\right)$ after we have defined the Poisson transformation, which does not break the Poisson structure of the transformation itself.

Definition 1. 

CO-LPNets A composition of Poisson mappings of phase space ${\mathbf{P}}_k$ generated by Hamiltonians $h_k$ defined by (22) over time $\Delta t$ forms a Poisson transformation. We make the parameter $a_k$ to be defined as a function of the initial condition for every step ${\breve{\mu }}_0$ through a neural network $w_k=w_{k,NN}\left({\breve{\mu }}_0\right)$. The total neural network is constructed as the composition of ${\mathbf{P}}_K\circ {\mathbf{P}}_{k-1}\circ \dots {\mathbf{P}}_1$, with the parameter $a_k$ in each transformation depending on an individual neural network. The neural network constructed by such composition is called CO-LPNets . The construction of CO-LPNets is illustrated on Figure 1.

Remark 4 

(On the type of networks for $w_k$). We take a shallow network with the input being all components of ${\breve{\mu }}_k$ and the output being one-dimensional parameter $w_k$, and one hidden layer of width W. There are K identical networks, one predicting each value of $w_k=w_{k,NN}\left(\breve{\mu }\right)$, $k=1,\dots ,N$. The number of parameters in each network is thus $N_p^0=(N·n·W)+2W+1$, and the total number of parameters is $N_{p,tot}=K{N}_p^0$. We will be taking K to be a multiple of the total dimensions $Nn$ (usually equal to $Nn$), so the number of parameters in the network grows with the number of particles and dimensions of the group as $n^2{N}^{2}W$.

In our observations, we take W to be quite small. In particular, in our simulations for $N=3$ particles, we take $W=3$, and $K=Nn$. For particles evolving in $\mathrm{SO}\left(3\right)$ group, the number of parameters in each network is 34, with the total number of parameters being 306. For $\mathrm{SE}\left(3\right)$ group, each network contains 61 parameters, with the total number of parameters equal to 1098.

On the Completeness of CO-LPNets

To clarify the completeness of the CO-LPNets architecture in relation to Lemma 2, we emphasize the geometric structure of the Lie–Poisson equations. The true dynamics are restricted to a symplectic leaf (a coadjoint orbit) within the dual Lie algebra ${\mathfrak{g}}^{*}$. The tangent space to this orbit at any point $\breve{\mu }$ is exactly the image of the Lie–Poisson matrix $\mathsf{\Lambda }\left(\breve{\mu }\right)$.

By choosing the test Hamiltonians $h_k=w_k{\mu }_{j_k}$, the gradient is $\nabla h_k=w_k{\mathbf{e}}_{j_k}$. Consequently, the resulting vector fields $V_k=w_k\mathsf{\Lambda }\left(\breve{\mu }\right){\mathbf{e}}_{j_k}$ correspond to (infinitesimally) flowing along the individual columns of $\mathsf{\Lambda }\left(\breve{\mu }\right)$. As the index $j_k$ cycles through all $Nn$ coordinates, the network explicitly spans every independent basis direction of the coadjoint orbit’s tangent space.

Therefore, learning the state-dependent weights $w_{k,NN}\left(\breve{\mu }\right)$ is equivalent to learning the optimal projection of the true Hamiltonian vector field onto this local basis. Because the composition of flows along a complete set of basis vectors can approximate the flow of any arbitrary vector field in that tangent space, analogous to general splitting methods [47], the finite composition of our transformations ${\mathbf{P}}_k$ is theoretically capable of approximating any dynamically valid phase flow. This geometric approach bridges the theoretical existence of a separable approximation derived in Lemma 2 with a highly efficient practical implementation, explicitly bypassing the potentially astronomical parameter counts required by the generic single-layer separable ansatz.

On Possible Simplifications

Several simplifications of the CO-LPNets can be constructed. One can remove the middle layer and make the output neuron dependent only on one input neuron. Then, each network will contain only two parameters (three in the case of additional scaling), and the total number of parameters will be proportional to K. This was the idea of G-LPNets in [3]. Unfortunately, for Lie-Poisson systems with block-diagonal matrices, such systems may have problems with convergence as Poisson transformations ${\mathbf{P}}_k$ will be isolated to different particles. A solution to this problem is to consider the output neuron depending on a pair of input neurons. Such mapping was implemented for this problem as a test case. However, we believe that the convergence of that algorithm will be slower, since it is not completely clear a priori how many Poisson transformations of this type we will need to achieve convergence.

The Algorithm of CO-LPNets

• Generate several short trajectories of ground truth data. The data will define the ground truth transformations ${\mu }_i^0\to {\mu }_i^f$, $i=1,\dots ,N_p$.
• We construct the transformations ${\mathbf{P}}_k(w_k,{\breve{\mu }}_0)$ and compute the loss function

  $$\begin{array}{cc}\hfill L\left(\overline{\mathbf{W}}\right)& ={\displaystyle \sum _{j=1}^M}{∥{\breve{\mu }}_j^f-{\mathbf{P}}_K\left(w_K\right)\circ {\mathbf{P}}_{K-1}\left(w_{K-1}\right)\circ \dots \circ {\mathbf{P}}_1\left(w_1\right)\left({\breve{\mu }}_j^0\right)∥}^2,\hfill \\ \hfill w_j& =w_{j,NN}\left({\mathbf{W}}_j\right),\overline{\mathbf{W}}:=({\mathbf{W}}_1,{\mathbf{W}}_2,\dots ,{\mathbf{W}}_K),\hfill \end{array}$$

  (23)

  where $\overline{\mathbf{W}}$ stands for weights and biases of all networks defining $w_{k,NN}$, $k=1,\dots K$.
• Find the optimal weights

  $${\overline{\mathbf{W}}}_{*}=\arg \min L\left(\overline{\mathbf{W}}\right)$$

  (24)
• Start with an initial condition ${\breve{\mu }}_0$, and construct $N_r$ reconstruction steps:

  $$\begin{array}{cc}\hfill {\breve{\mu }}_{\alpha +1}& ={\mathbf{P}}_K\left(w_K\right)\circ {\mathbf{P}}_{K-1}\left(w_{K-1}\right)\circ \dots \circ {\mathbf{P}}_1\left(w_1\right)\left({\breve{\mu }}_{\alpha }\right),0\le \alpha \le N_r,\hfill \\ \hfill w_k& =w_{k,NN}\left({\overline{\mathbf{W}}}_{*}\right).\hfill \end{array}$$

  (25)
  Compare the reconstructed solution with the ground truth solution obtained by the Lie-Poisson integrator.

Remark 5 

(On the structure of Poisson transformations). In what follows, we demonstrate that for the cases we consider, namely $G=\mathrm{SO}\left(3\right)$ and $G=\mathrm{SE}\left(3\right)$, the Poisson transformations ${\mathbf{P}}_k(w_k,{\breve{\mu }}_0)$ have the form $\mathbb{A}(a_k,w_k){\breve{\mu }}_0$, where ${\mathbb{A}}_k$ is a matrix of a certain form. The transformation is nonlinear since $w_k=w_{k,NN}\left({\breve{\mu }}_0\right)$. That structure of transformation is a consequence of the solutions to (22). Since every ${\mathbb{A}}_k(\breve{\mu },w_k)$ only depends on its own network through a single parameter $w_k=w_{k,NN}\left(\breve{\mu }\right)$, the optimization of the loss function (23) can be made highly computationally efficient.

6. Specialization to Particular Groups

6.1. $\mathrm{SO}\left(3\right)$ Group

6.1.1. Equations for the Ground Truth Calculations

Let us now apply the general theory to a particular group. We start with the collective motion of objects, with each object evolving on the Lie group of rotations $\mathrm{SO}\left(3\right)$. The $\mathrm{SO}\left(3\right)$ group consists of $3\times 3$ orthogonal matrices $\mathbb{Q}$ such that ${\mathbb{Q}}^{\mathsf{T}}\mathbb{Q}=\mathbb{Q}{\mathbb{Q}}^{\mathsf{T}}={\mathbb{I}}_3$, where ${\mathbb{I}}_3$ is a $3\times 3$ identity matrix, and $det\mathbb{Q}=1$. This example can be applied to a set of rigid objects rotating to achieve the desired orientation, for example, for observations of a particular target in space.

To find the evolution of particles on the N copies of Lie algebra $\mathfrak{so}\left(3\right)$, we defined the scalar product $<X,Y>=Tr\left(X^{\mathsf{T}}Y\right)$ and use the following orthogonal basis:

$$\begin{array}{c}{X}_1={\displaystyle \frac{1}{\sqrt{2}}}\left[\begin{array}{ccc}0& 0& \phantom{-}0\\ 0& 0& -1\\ 0& 1& \phantom{-}0\end{array}\right],X_2={\displaystyle \frac{1}{\sqrt{2}}}\left[\begin{array}{ccc}0& 0& 1\\ 0& 0& 0\\ -1& 0& 0\end{array}\right],X_3={\displaystyle \frac{1}{\sqrt{2}}}\left[\begin{array}{ccc}0& -1& 0\\ 1& 0& 0\\ 0& 0& 0\end{array}\right],\end{array}$$

which is normalized so that $\parallel X_j\parallel =1$, $j=1,2,3$ with respect to the trace norm. The structure constants for the $\mathfrak{so}\left(3\right)$ Lie algebra are ${\mathsf{\Gamma }}_{ij}^k={\displaystyle \frac{1}{\sqrt{2}}}{\epsilon }_{ijk},$ where ${\epsilon }_{ijk}$ is the Levi-Civita symbol.

For each particle, we introduce the vector of ${\mathit{\mu }}_k$ which consists of the components of ${\widehat{\mu }}_k$ in the basis $\left\{X_1,\dots X_n\right\}$:

$$\begin{array}{c}{\mathit{\mu }}_k={({\mu }_{k1},{\mu }_{k2},{\mu }_{k3})}^{\mathsf{T}},k=1,\dots ,N,\iff {\widehat{\mu }}_k={\displaystyle \sum _{j=1}^3}{\mu }_{kj}{X}_j.\end{array}$$

(26)

Thus, according to (12), $\mathsf{\Lambda }\left(\breve{\mu }\right)$ is a block-diagonal matrix of the following form:

$$\begin{array}{c}\mathsf{\Lambda }\left(\breve{\mu }\right)={\displaystyle \frac{1}{\sqrt{2}}}\left[\begin{array}{cccc}{\widehat{\mu }}_1& 0& \dots & 0\\ 0& {\widehat{\mu }}_2& \dots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \dots & {\widehat{\mu }}_N\end{array}\right]\end{array}$$

with k-th block equal to the following $3\times 3$ antisymmetric matrix:

$$\begin{array}{c}{\widehat{\mu }}_k=\left[\begin{array}{ccc}0& -{\mu }_{k3}& {\mu }_{k2}\\ {\mu }_{k3}& 0& -{\mu }_{k1}\\ -{\mu }_{k2}& {\mu }_{k1}& 0\end{array}\right],k=1,\dots ,N,\end{array}$$

(27)

which is just the familiar hat map, mapping vectors in ${\mathbb{R}}^3$ to elements of Lie algebra $\mathfrak{so}\left(3\right)$. As a result of the decoupled Poisson tensor, there is a Casimir $c_k$ for each particle with the index $k=1,\dots N$:

$$\begin{array}{c}\hfill c_k={\mu }_{k1}^2+{\mu }_{k2}^2+{\mu }_{k3}^2=\left|\right|{\mu }_k{\left|\right|}^2,k=1,\dots ,N.\end{array}$$

The Lie–Poisson equations (12) are written in vector form as:

$$\begin{array}{c}{\dot{\mu }}_k={\displaystyle \frac{1}{\sqrt{2}}}{\mu }_k\times {\displaystyle \frac{\partial h}{\partial {\mu }_k}},k=1,\dots ,N,\end{array}$$

which are (formally) N copies of the equations of motion a rigid bodies [21], with the Hamiltonian $h=h({\mu }_1,\dots ,{\mu }_N)$ providing the interaction between the particles.

We assume that the control is applied in the $X_1$ direction and the drift is in the $X_2$ direction. Then, (7) becomes

$${\xi }_k=X_2+u_{k1}{X}_1,k=1,\dots ,N,$$

(28)

After some algebra, the Hamiltonian (A15) for the cases of ’Dictatorship’ yields:

$$\begin{array}{cc}\hfill h=& {\displaystyle \sum _{k=1}^N}{\mu }_{k2}+{\displaystyle \frac{1}{2}}({\displaystyle \frac{1+2\chi }{1+2N\chi }}{\mu }_{11}^2+{\displaystyle \frac{1+2N\chi +4{\chi }^2}{(1+2N\chi )(1+2\chi )}}{\displaystyle \sum _{j=2}^N}{\mu }_{j1}^2\hfill \\ \hfill & +{\displaystyle \frac{4\chi }{1+2N\chi }}{\mu }_{11}{\displaystyle \sum _{j=2}^N}{\mu }_{j1}+{\displaystyle \frac{8{\chi }^2}{(1+2N\chi )(1+2\chi )}}{\displaystyle \sum _{i,j=2,i<j}^N}{\mu }_{i1}{\mu }_{j1}).\hfill \end{array}$$

(29)

Similarly, for the case of ’Democracy’, we obtain

$$\begin{array}{c}\hfill h={\displaystyle \sum _{k=1}^N}{\mu }_{k2}+{\displaystyle \frac{1}{2}}\left({\displaystyle \frac{1+2\chi }{1+2N\chi }}{\displaystyle \sum _{j=1}^N}{\mu }_{j1}^2+{\displaystyle \frac{4\chi }{1+2N\chi }}{\displaystyle \sum _{i,j=1,i<j}^N}{\mu }_{i1}{\mu }_{j1}\right).\end{array}$$

(30)

These two cases (’Dictatorship’ and ’Democracy’) are simulated for $N=3$ particles, computing 40 trajectories of length 51, thus creating 2000 intervals for learning of the phase space. The initial conditions for these trajectories are chosen randomly in the hypercube ${[-1,1]}^d$, where d is the dimension of the phase space, equal to $Nn=9$ in this case. As mentioned in Section 4.2, we use a Lie-Poisson integrator to compute the ground truth data, which preserves the system’s structure and Casimirs with machine precision. We choose the value of the parameter $\chi =0.5$ for all simulations. For reference, we present analytical formulas for the derivatives of the Hamiltonians in Appendix C.1.

Next, we describe how to construct the Poisson transformations that define the transformations of the phase space. We then build CO-LPNets to learn the motion of the system in the whole space where the data is available, and predict the trajectories of the system that are previously unseen by the algorithm.

6.1.2. Poisson Transformations Through Test Hamiltonians

Following Section 5.2, we consider test Hamiltonians which depend only on a single component i of one particle k. Following (22), we start with the Hamiltonian $h_{ki}=w_{ki}{\mu }_{ki}$, where $w_{ki}$ is some constant on that particular time step.

The Lie–Poisson reduced dynamics equations $\dot{\breve{\mu }}=\mathsf{\Lambda }\left(\breve{\mu }\right)\nabla h$ (A16) simplify to give

$${\dot{\mathit{\mu }}}_k={\mathit{\mu }}_k\times {\mathbf{e}}_i{w}_{ki},$$

(31)

where ${\mathbf{e}}_i$ means the i-th standard basis vector in the 3-dimensional space. Assuming that we are given initial conditions for $\breve{\mu }\left(0\right)=({\mathit{\mu }}_1\left(0\right),\dots ,{\mathit{\mu }}_N\left(0\right))$, the solution of such an initial value problem is simply

$${\mathit{\mu }}_k\left(t\right)=\mathbb{R}({\mathbf{e}}_i,-w_{ki}t){\mathit{\mu }}_k\left(0\right),{\mathit{\mu }}_j\left(t\right)={\mathit{\mu }}_j\left(0\right),j\ne k.$$

(32)

where $\mathbb{R}(\mathbf{n},\phi )$ denoting the rotation matrix with respect to the vector $\mathbf{n}$ by the angle $\phi$.

To speed up computation of derivatives of the loss function (23) with respect to the parameters, we have also analytically computed the partial derivatives of matrices $\mathbb{R}$, $i=1,2,3$ with respect to $w_{ki}$. We present these formulas in the Section C.1 of the Appendix.

6.1.3. Application of CO-LPNets for $\mathrm{SO}\left(3\right)$ Group

We generate the data in phase space, as was described in (4.2). We consider $N=3$ coupled particles each evolving in $\mathrm{SO}\left(3\right)$, with the dimension of each component being $n=3$. The total dimension of the space is 9, and the number of generated intervals is $2,000$ with the time step $\Delta t=0.1$. Data is produced in two cases of particle interaction: ’Dictatorship’ and ’Democracy’. We then apply CO-LPNets with $Nn=9$ transformations, the minimum possible number. Each transformation is parameterized by a neural network with 9 input neurons, $W=3$ hidden neurons, and one output neuron. We then optimize the loss defined by (23). We perform an optimization using the Adam algorithm with a learning rate of $0.005$ for 10,000 epochs. During that time, the loss decreases from about $\sim 0.1-0.01$ to $\sim {10}^{-6}$, i.e. by four to five orders of magnitude. After the optimization is finished, we simulate ten trajectories starting with random initial conditions in the cube ${[-1,1]}^{Nn}$, both using the structure-preserving integrator the CO-LPNets with computed parameters.

The results of application of CO-LPNets are shown in Figure 2 for the case of ’Democracy’, with the corresponding summary of errors in Casimirs, energy, and MAE of all components in Figure 3. The same information is shown on Figure 4 for the case of ’Dictatorship’, with the corresponding summary of errors in Casimirs, energy, and MAE of all components in Figure 5. The trajectories are reproduced faithfully over 1000 time steps for every particle, the variations of energy in trajectories by CO-LPNets are bounded, and the Casimirs for each particle are preserved with machine precision.

In all Figures, when we compare the results of the CO-LPNets with the ground truth, we use a consistent color scheme: the ground truth is marked with a solid blue line, and results obtained from CO-LPNets are presented with a solid red line.

6.2. $\mathrm{SE}\left(3\right)$ Group

6.2.1. Equations for Ground Truth Calculations

Definitions

We now proceed to the consideration of equations of motion for the Lie group of rotations and translations. This group gives a highly practical example as it describes the motion of several interacting objects moving in space, exhibiting both rotations and translations. This example is useful for drone applications. An object in space evolving on the group $\mathrm{SE}\left(3\right)$ undergoes rotation $\mathbb{Q}\in \mathrm{SO}\left(3\right)$ and translation $\mathbf{v}\in {\mathbb{R}}^3$. Using these quantities $(\mathbb{Q},\mathbf{v})$, we form a $4\times 4$ matrix $g\in \mathrm{SE}\left(3\right)$, and compute its inverse $g^{-1}$ as follows:

$$g=\left(\begin{array}{cc}\mathbb{Q}& \mathbf{v}\\ {\mathbf{0}}^{\mathsf{T}}& 1\end{array}\right);g_1{g}_2=\left(\begin{array}{cc}{\mathbb{Q}}_1{\mathbb{Q}}_2& {\mathbb{Q}}_1{\mathbf{v}}_2+{\mathbf{v}}_1\\ {\mathbf{0}}^{\mathsf{T}}& 1\end{array}\right);g^{-1}=\left(\begin{array}{cc}{\mathbb{Q}}^{\mathsf{T}}& -{\mathbb{Q}}^{\mathsf{T}}\mathbf{v}\\ {\mathbf{0}}^{\mathsf{T}}& 1\end{array}\right),$$

(33)

where ${\mathbf{0}}^{\mathsf{T}}$ is the row of zeros. The multiplicative structure (33) defines a semidirect product group. This group $SE\left(3\right)$ has a corresponding Lie algebra, which we call $\mathfrak{se}\left(3\right)$. That Lie algebra is six-dimensional; three values are coming from the infinitesimal rotations and three values come from translations.

In order to find an orthonormal basis in the linear space $\mathfrak{se}\left(3\right)$, we describe the elements of $\mathfrak{se}\left(3\right)$ as $4\times 4$ matrices and again employ the trace product of two matrices defined as $<X,Y>=Tr\left(X^{\mathsf{T}}Y\right)$. With this product, we select the following orthogonal basis for the Lie algebra $\mathfrak{se}\left(3\right)$:

$$\begin{array}{c}{G}_1={\displaystyle \frac{1}{\sqrt{2}}}\left[\begin{array}{cccc}0& 0& 1& 0\\ 0& 0& 0& 0\\ -1& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right],G_2={\displaystyle \frac{1}{\sqrt{2}}}\left[\begin{array}{cccc}0& -1& 0& 0\\ 1& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right],G_3={\displaystyle \frac{1}{\sqrt{2}}}\left[\begin{array}{cccc}0& 0& 0& 0\\ 0& 0& -1& 0\\ 0& 1& 0& 0\\ 0& 0& 0& 0\end{array}\right],\\ G_4=\left[\begin{array}{cccc}0& 0& 0& 0\\ 0& 0& 0& 1\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right],G_5=\left[\begin{array}{cccc}0& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 1\\ 0& 0& 0& 0\end{array}\right],G_6=\left[\begin{array}{cccc}0& 0& 0& 1\\ 0& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right],\end{array}$$

which is normalized so that $\parallel G_j\parallel =1$, $j=1,\dots ,6$ with respect to the trace norm. The structure constants for the $\mathfrak{se}\left(3\right)$ Lie algebra are then computed as follows:

$$\begin{array}{c}{\mathsf{\Gamma }}_{61}^5={\displaystyle \frac{1}{\sqrt{2}}},{\mathsf{\Gamma }}_{16}^5=-{\displaystyle \frac{1}{\sqrt{2}}},{\mathsf{\Gamma }}_{62}^4=-{\displaystyle \frac{1}{\sqrt{2}}},{\mathsf{\Gamma }}_{26}^4={\displaystyle \frac{1}{\sqrt{2}}},{\mathsf{\Gamma }}_{43}^5=-{\displaystyle \frac{1}{\sqrt{2}}},{\mathsf{\Gamma }}_{34}^5={\displaystyle \frac{1}{\sqrt{2}}},\\ {\mathsf{\Gamma }}_{42}^6={\displaystyle \frac{1}{\sqrt{2}}},{\mathsf{\Gamma }}_{24}^6=-{\displaystyle \frac{1}{\sqrt{2}}},{\mathsf{\Gamma }}_{53}^4={\displaystyle \frac{1}{\sqrt{2}}},{\mathsf{\Gamma }}_{35}^4=-{\displaystyle \frac{1}{\sqrt{2}}},{\mathsf{\Gamma }}_{51}^6=-{\displaystyle \frac{1}{\sqrt{2}}},{\mathsf{\Gamma }}_{15}^6={\displaystyle \frac{1}{\sqrt{2}}},\\ {\mathsf{\Gamma }}_{31}^2={\displaystyle \frac{1}{\sqrt{2}}},{\mathsf{\Gamma }}_{13}^2=-{\displaystyle \frac{1}{\sqrt{2}}},{\mathsf{\Gamma }}_{32}^1=-{\displaystyle \frac{1}{\sqrt{2}}},{\mathsf{\Gamma }}_{23}^1={\displaystyle \frac{1}{\sqrt{2}}},{\mathsf{\Gamma }}_{12}^3={\displaystyle \frac{1}{\sqrt{2}}},{\mathsf{\Gamma }}_{21}^3=-{\displaystyle \frac{1}{\sqrt{2}}}\end{array}$$

with all other structure constants equal to 0.

Using the mechanical analogy from the mechanical Lie-Poisson systems [44,45], let us separate the full momentum vector for each particle ${\mathit{\mu }}_k$ into two three-dimensional vectors ${\mathsf{\Pi }}_k$ and ${\mathbf{p}}_k$ having the meaning of the angular and linear momenta:

$$\begin{array}{c}\mathrm{vector}\mathrm{of}\mathrm{angular}\mathrm{momenta}:{\mathsf{\Pi }}_k={({\mu }_{k1},{\mu }_{k2},{\mu }_{k3})}^{\mathsf{T}},k=1,\dots ,N\\ \mathrm{vector}\mathrm{of}\mathrm{linear}\mathrm{momenta}:{\mathbf{p}}_k={({\mu }_{k4},{\mu }_{k5},{\mu }_{k6})}^{\mathsf{T}},k=1,\dots ,N.\end{array}$$

(34)

Thus, according to (12), $\mathsf{\Lambda }\left(\breve{\mu }\right)$ is a block-diagonal matrix of the following form:

$$\begin{array}{c}\mathsf{\Lambda }\left(\breve{\mu }\right)={\displaystyle \frac{1}{\sqrt{2}}}\left[\begin{array}{cccc}{\widehat{\mu }}_1& 0& \dots & 0\\ 0& {\widehat{\mu }}_2& \dots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& \dots & 0& {\widehat{\mu }}_N\end{array}\right]\end{array}$$

with k-th block equal to the following $6\times 6$ antisymmetric matrix:

$$\begin{array}{c}{\widehat{\mu }}_k=\left[\begin{array}{cccccc}0& -{\mu }_{k3}& {\mu }_{k2}& 0& -{\mu }_{k6}& {\mu }_{k5}\\ {\mu }_{k3}& 0& -{\mu }_{k1}& {\mu }_{k6}& 0& -{\mu }_{k4}\\ -{\mu }_{k2}& {\mu }_{k1}& 0& -{\mu }_{k5}& {\mu }_{k4}& 0\\ 0& -{\mu }_{k6}& {\mu }_{k5}& 0& 0& 0\\ {\mu }_{k6}& 0& -{\mu }_{k4}& 0& 0& 0\\ -{\mu }_{k5}& {\mu }_{k4}& 0& 0& 0& 0\end{array}\right]=\left[\begin{array}{cc}{\widehat{\mathsf{\Pi }}}_k& {\widehat{p}}_k\\ {\widehat{p}}_k& ⊬\end{array}\right],k=1,\dots ,N,\end{array}$$

(35)

where $⊬$ is a $3\times 3$ matrix of zeros and ${\widehat{\alpha }}_{ij}=-{\epsilon }_{ijk}{\alpha }_k$ is, as usual, the hat map between the components of the vector $\alpha$ and the antisymmetric matrix $\widehat{\alpha }$, just as in (27). The matrix $\widehat{\alpha }$ defined this way posesses the property $\widehat{\alpha }\mathbf{v}=\alpha \times \mathbf{v}$ for any vector $\mathbf{v}\in {\mathbb{R}}^3$. As a result of the decoupled Poisson tensor, there are $2N$ Casimirs, with two Casimirs for each particle:

$$\begin{array}{c}\hfill C_{1,k}={\mu }_{k4}^2+{\mu }_{k5}^2+{\mu }_{k6}^2=\left|\right|{\mathbf{p}}_k{\left|\right|}^2,k=1,\dots ,N,\\ \hfill C_{2,k}={\mu }_{k1}{\mu }_{k4}+{\mu }_{k2}{\mu }_{k5}+{\mu }_{k3}{\mu }_{k6}={\mathsf{\Pi }}_k·{\mathbf{p}}_k,k=1,\dots ,N.\end{array}$$

The Lie–Poisson equations (12) reduce to the well-known form:

$$\begin{array}{c}{\dot{\mathsf{\Pi }}}_k={\displaystyle \frac{1}{\sqrt{2}}}{\mathsf{\Pi }}_{\mathbf{k}}\times {\displaystyle \frac{\partial h}{\partial {\mathsf{\Pi }}_k}}+{\displaystyle \frac{1}{\sqrt{2}}}{\mathbf{p}}_{\mathbf{k}}\times {\displaystyle \frac{\partial h}{\partial {\mathbf{p}}_k}},k=1,\dots ,N\\ {\dot{\mathbf{p}}}_k={\displaystyle \frac{1}{\sqrt{2}}}{\mathbf{p}}_{\mathbf{k}}\times {\displaystyle \frac{\partial h}{\partial {\mathsf{\Pi }}_k}},k=1,\dots ,N\end{array}$$

which are the equations for N underwater vehicles when the centers of mass and buoyancy coincide, extending the classical single-vehicle case developed by Kirchhoff [52,53,54]. We assume that the drift is along the first component of linear momentum (the drone moves forward) and there are two controls in angular momentum in the first and second coordinate. The equation (7) becomes

$${\xi }_k=G_4+u_{k1}{G}_1+u_{k2}{G}_2,k=1,\dots ,N.$$

(36)

Using the general formula (11), after some algebra, we compute the explicit formulas for the Hamiltonian for the case of ’Dictatorship’

$$\begin{array}{cc}\hfill h={\displaystyle \sum _{k=1}^N}{\mu }_{k4}& +{\displaystyle \frac{1}{2}}({\displaystyle \frac{1+2\chi }{1+2N\chi }}\left({\mu }_{11}^2+{\mu }_{12}^2\right)+{\displaystyle \frac{1+2N\chi +4{\chi }^2}{(1+2N\chi )(1+2\chi )}}{\displaystyle \sum _{k=2}^N}\left({\mu }_{k1}^2+{\mu }_{k2}^2\right)\hfill \\ \hfill & +{\displaystyle \frac{4\chi }{1+2N\chi }}\left({\mu }_{11}{\displaystyle \sum _{k=2}^N}{\mu }_{k1}+{\mu }_{12}{\displaystyle \sum _{k=2}^N}{\mu }_{k2}\right)\hfill \\ \hfill & +{\displaystyle \frac{8{\chi }^2}{(1+2N\chi )(1+2\chi )}}{\displaystyle \sum _{i,j=2,i<j}^N}\left({\mu }_{i1}{\mu }_{j1}+{\mu }_{i2}{\mu }_{j2}\right))\hfill \end{array}$$

(37)

For the case of ’Democracy’, the Hamiltonian is given by

$$\begin{array}{cc}\hfill h={\displaystyle \sum _{k=1}^N}{\mu }_{k4}& +{\displaystyle \frac{1}{2}}({\displaystyle \frac{1+2\chi }{1+2N\chi }}{\displaystyle \sum _{k=1}^N}\left({\mu }_{k1}^2+{\mu }_{k2}^2\right)\hfill \\ \hfill & +{\displaystyle \frac{4\chi }{1+2N\chi }}{\displaystyle \sum _{i,j=1,i<j}^N}\left({\mu }_{i1}{\mu }_{j1}+{\mu }_{i2}{\mu }_{j2}\right)).\hfill \end{array}$$

(38)

The forms for the Hamiltonian (37) and (38) are based on the explicit form of the matrix $\mathsf{\Psi }$, which we present in Section C.2.

Data Generation for Learning

We generate data for two cases, ’Dictatorship’ and ’Democracy’, by simulating the system with $N=3$ particles, creating 80 trajectories of length 51, thus creating 4000 intervals for learning of the phase space. The initial conditions for these trajectories are also chosen randomly in the hypercube ${[-1,1]}^d$, where d, the dimension of the phase space, is now equal to $Nn=18$ in this case. Again, as mentioned in Section 4.2, we use a Lie-Poisson integrator to compute the ground truth data, which preserves the system’s structure and Casimirs with machine precision. We choose the value of the parameter $\chi =0.5$ for all simulations. We also present analytical formulas for the matrices $\mathsf{\Psi }$ and derivatives of the Hamiltonians in Appendix C.2.

Next, we construct Poisson transformations as flows generated by test Hamiltonians in order to create Poisson transformations for CO-LPNets.

6.2.2. Poisson Transformations Through Test Hamiltonians

To construct the Poisson transformations for CO-LPNets in the case of $\mathrm{SE}\left(3\right)$ group we consider Hamiltonians which depend only on a single component i, $1\le i\le 6$ of one particle k, $1\le k\le N$; more precisely, we take $h_{ki}=w_{ki}{\mu }_{ki}$ To construct the Poisson transformations explicitly, we consider two possibilities: $1\le i\le 3$ (test Hamiltonian corresponding to the angular momentum) and $4\le i\le 6$ (test Hamiltonian corresponding to the linear momentum).

Test Hamiltonians for Angular Momenta

This is the case $h_{ki}=w_{ki}{\mu }_{ki}$ with $1\le i\le 3$. Equations (12) simplify as follows:

$$\begin{array}{cc}\hfill & {\dot{\mathsf{\Pi }}}_k={\mathsf{\Pi }}_k\times {\mathbf{e}}_i{w}_{ki},{\dot{\mathit{p}}}_k={\mathit{p}}_k\times {\mathbf{e}}_i{w}_{ki},\hfill \\ \hfill & {\dot{\mathsf{\Pi }}}_j=\mathbf{0},{\dot{\mathbf{p}}}_j=\mathbf{0},j\ne k.\hfill \end{array}$$

(39)

with ${\mathsf{\Pi }}_k$ and ${\mathbf{p}}_k$ defined as in (34). Again, ${\mathbf{e}}_i$ means the i-th standard basis vector in the 3-dimensional space; $w_{ki}$ is a constant on every time step that will need to be approximated by the neural network. Suppose we are given initial conditions for $\mathsf{\Pi }\left(t\right)$ and $\mathbf{p}\left(t\right)$ at the point $t=0$, which we denote as ${\mathsf{\Pi }}_k\left(0\right)$ and ${\mathit{p}}_k\left(0\right)$. The solution of the initial value problem (39) is then given by

$$\begin{array}{cc}\hfill {\mathsf{\Pi }}_k\left(t\right)& =\mathbb{R}({\mathbf{e}}_i,-w_{ki}t){\mathsf{\Pi }}_k\left(0\right),{\mathbf{p}}_k\left(t\right)=\mathbb{R}({\mathbf{e}}_i,-w_{ki}t){\mathit{p}}_k\left(0\right),\hfill \\ \hfill {\mathsf{\Pi }}_j\left(t\right)& ={\mathsf{\Pi }}_j\left(0\right),{\mathbf{p}}_j\left(t\right)={\mathit{p}}_j\left(0\right),j\ne k.\hfill \end{array}$$

(40)

Test Hamiltonians for Linear Momenta

This is the case $h_{ki}=w_{ki}{\mu }_{ki}$ with $4\le i\le 6$. Equations (12) are now written as follows: for the kth particle, we get the system of equations:

$$\begin{array}{cc}\hfill & {\dot{\mathsf{\Pi }}}_k={\mathbf{p}}_k\times {\mathbf{e}}_{i-3}{w}_{ki},{\dot{\mathbf{p}}}_k=\mathbf{0},\hfill \\ \hfill & {\dot{\mathsf{\Pi }}}_j=\mathbf{0},{\dot{\mathbf{p}}}_j=\mathbf{0},j\ne k.\hfill \end{array}$$

(41)

with the solution

$$\begin{array}{cc}\hfill {\mathsf{\Pi }}_k\left(t\right)& ={\mathbf{p}}_k\times {\mathbf{e}}_{i-3}{w}_{ki}t+{\mathsf{\Pi }}_k\left(0\right),{\mathbf{p}}_k\left(t\right)={\mathit{p}}_k\left(0\right),\hfill \\ \hfill {\mathsf{\Pi }}_j\left(t\right)& ={\mathsf{\Pi }}_j\left(0\right),{\mathbf{p}}_j\left(t\right)={\mathit{p}}_j\left(0\right),j\ne k.\hfill \end{array}$$

(42)

The transformations (40) and (42) can be written as block-diagonal matrices with $j\ne k$ blocks being identities, and k-th block being a combination of rotations and translations. For brevity, we do not present them here. For the computations of the loss function (23), we present partial derivatives of these transformations with respect to parameters $w_{ki}$ in Section C.2 of the Appendix.

6.2.3. Results from the Application of CO-LPNets for Particles Evolving on $\mathrm{SE}\left(3\right)$

We generate the data in phase space, as was described in (4.2). We consider $N=3$ coupled particles each evolving in $\mathrm{SE}\left(3\right)$, having dimension of $n=6$ with the total dimension of the phase space being $Nn=18$. We generate $4,000$ points, to match the number of points per dimension with the $\mathrm{SO}\left(3\right)$ case described in Section 6.1.3, with $\sim 222$ points per dimension of the phase space. Again, the time step is chosen to be $\Delta t=0.1$. We produce data for learning in two cases: ’Dictatorship’ and ’Democracy’. We then apply CO-LPNets with 18 rotations, the minimum number of rotations allowed to cover all possible motions in space, and optimize the loss defined by (23). We perform the optimization using the Adam algorithm with a learning rate of $0.005$ for 10,000 epochs. During that time, the loss decreases from about $\sim 0.1$ to $\sim {10}^{-5}$, i.e. by four orders of magnitude. After the optimization is finished, we simulate ten trajectories starting with random initial conditions in the cube ${[-1,1]}^{Nn}$, both using the structure-preserving integrator, the CO-LPNets with computed parameters.

The results of application of CO-LPNets are shown in Figure 6 for the case of ’Democracy’, with the corresponding summary of errors in Casimirs, energy, and MAE of all components in Figure 7. The same information is shown on Figure 8 for the case of , with the corresponding summary of errors in Casimirs, energy and MAE of all components in Figure 9. We see that CO-LPNets faithfully reproduce the trajectory over an intermediate time scale (200 steps), preserving Casimirs with machine precision. The Hamiltonian (energy) is not conserved by CO-LPNets, but the deviations of energy are bounded.

6.3. Studies of the Effects of Noise in Data on Method Accuracy

Since one cannot expect the data forming the input to the neural networks to be perfect, studying the robustness of the system with respect to noise is of ultimate importance for practical applications. We study the case where the underlying system is fully deterministic, but the collected data representing the system’s evolution are only known up to a certain accuracy.

More precisely, we perturb the initial data by a small amount of additive noise. We assume that the noise follows a zero-mean normal distribution $\mathcal{N}(0,{\eta }^2)$, where the noise amplitude (standard deviation) is chosen from the discrete set $\eta \in \{0.001,0.002,0.003,0.004,0.005\}$. We train the network on this noisy data, reconstruct $N=10$ test trajectories, and compare them against the unperturbed ground truth. To ensure statistical robustness, we evaluate five independent noise realizations for every value of $\eta$ and average the final metrics.

For brevity, we present the results of this study exclusively for the $SO\left(3\right)$ group with three interacting particles under the democracy topology, as previously detailed in Sec. 6.1. In Figure 10, we present the reconstructed trajectories from the noisy data for all nine momentum components. Visually, the trajectories generated under small $\eta$ values remain nearly indistinguishable from the exact ground truth.

In order to rigorously analyze the system’s performance, we evaluate the integration errors and conservation properties as a function of the noise level (Figure 11). As expected by the geometric nature of the architecture, all Casimir invariants are preserved to machine precision irrespective of the noise level (left panel). However, examining the trajectory accuracy reveals an interesting phenomenon: the reconstruction actually improves for small amounts of injected noise ($\eta =0.001$ and $0.002$) compared to the noiseless baseline ($\eta =0$), before eventually degrading for $\eta \ge 0.003$.

This counter-intuitive result is a manifestation of noise-induced regularization. When trained on perfectly clean data, the network can overfit to specific data points present in the training set. By injecting a small amount of noise, the learning procedure is forced to smooth out these isolated errors and learn the true, underlying continuous vector field representing the Hamiltonian gradients. Consequently, the network trades exact precision on scattered training points for improved global generalization across the continuous phase space. Thus, introducing a small amount of noise during training acts as a regularizer and is actively advantageous to the long-term predictive accuracy of our structure-preserving system.

7. Conclusions and Future Work

The main contribution of this paper is the presentation of a novel method for learning the full phase-space behavior of general Lie-Poisson systems, particularly those relevant to optimal control. Our models are highly efficient, utilizing a relatively small number of trainable parameters and requiring a minimal amount of data. We discussed the completeness of the learned transformations and established a formal connection between our work and the established field of structure-preserving splitting methods. We call our method CO-LPNets. The illustrative examples of predicting the evolution of interacting particles on the $\mathrm{SO}\left(3\right)$ and $\mathrm{SE}\left(3\right)$ groups demonstrate the potential of our method and its applicability to higher-dimensional systems. By design, CO-LPNets strictly preserve the underlying Poisson structure, thereby conserving Casimir invariants to machine precision. Furthermore, we demonstrated the robustness of the method with respect to noise, highlighting the advantage that a small amount of noise in the training data can actively improve the accuracy of long-term predictions.

In the future, we aim to extend this progress to the optimal control of interacting particles on the $\mathrm{SU}\left(2\right)$ group, which is essential for quantum computing. Another compelling problem is the interaction of several agents evolving on different Lie groups, where one group is a subgroup of another. For example, a flying drone evolves on the $\mathrm{SE}\left(3\right)$ group, and a land-based or sea-surface drone evolves on the $\mathrm{SE}\left(2\right)$ group, which is a subgroup of $\mathrm{SE}\left(3\right)$. Thus, we can consider the evolution of such a control system on the direct group product $\mathrm{SE}\left(3\right)\times \mathrm{SE}\left(2\right)$. Theoretical studies for such work have been undertaken in [55]; it would be interesting to extend our machine learning methods to these more general and practically important cases.

Finally, our methods can be applied to the rigorous, structure-preserving Reduced Order Modeling (ROM) of Lie-Poisson systems. While recent literature allows for the rigorous construction of ROMs for canonical systems [56,57,58], there is no apparent mathematical mechanism to reduce a high-dimensional Lie-Poisson bracket to a lower-dimensional space. Our method allows us to compute the dominant motions directly in the full Lie-Poisson space using the most ‘essential’ Poisson mappings for the dynamics, potentially leading to a radical reduction in computational complexity. This promising approach will be developed in future work.

Appendix A.1. The Lie-Poisson Theory of Optimal Control of Coupled Interacting Particles

This section outlines the theoretical framework for symmetry-reduced optimal control of N interacting particles. This chapter presented for reference to define the system and explain how the ground truth data for the application of our learning method was computed. Our derivation follows [1,2] quite closely, including the notation used, although we try to be concise here and skip a few details. Readers interested in the comprehensive derivations are encouraged to consult those primary sources.

We consider a system of N identical particles, where the configuration of each individual agent is represented by a Lie group G. Consequently, the total configuration manifold for the collective system is given by the direct product group $G^N$.

To reflect on the connection of this theory to mechanics, we first address symmetry reduction for mechanical systems before extending it to optimal control. A mechanical system on $G^N$ is traditionally characterized by the canonical Poisson bracket on the cotangent bundle ${\left(T^{*}G\right)}^N$. In this representation, the system state is defined by the coordinates $(g_1,\dots ,g_N)$ and their conjugate momenta ${\mathbf{p}}_1,\dots ,{\mathbf{p}}_N$.

As introduced in Section 1, the full symmetry reduction of such canonical systems yields Lie-Poisson systems (5) (see [45], Chapter 10). For any system defined on a Lie group G, there is an associated Lie algebra $\mathfrak{g}$ with a bracket $[\alpha ,\beta ]$ for elements $\alpha ,\beta \in \mathfrak{g}$. We define $\langle \mu ,\alpha \rangle$ as the duality pairing between the algebra $\mathfrak{g}$ and its dual space ${\mathfrak{g}}^{*}$, where the momentum co-vector is $\mu =g^{-1}p$. For functions $F,H:{\mathfrak{g}}^{*}\to \mathbb{R}$, the partial derivatives reside in $\mathfrak{g}$, forming the Lie–Poisson bracket structure seen in (5). Since our focus is on left-invariant dynamics, we adopt the negative sign convention in the bracket definition through the paper.

Pontryagin’s Maximum Principle: Linking Optimal Control and Hamiltonian Dynamics

The Pontryagin’s Maximum Principle provides the formal bridge between optimal control and Hamiltonian mechanics. Consider a system governed by the differential equation $\dot{x}=a(x\left(t\right),u\left(t\right),t)$ subject to control inputs $u\left(t\right)$. Here, x, and u may be scalars of vectors; if x is a vector, then a is a vector-valued function with the same number of components as x. Suppose we also have fixed boundary conditions $x\left(t_0\right)=x_0^{*}$ and $x\left(t_1\right)=x_1^{*}$. Then, the objective is to minimize the following cost functional:

$$\eta ={\int }_{t_0}^{t_1}L(x\left(\tau \right),\dot{x}\left(\tau \right),u\left(\tau \right),\tau )d\tau ,$$

(A1)

where the control Lagrangian L encodes optimality criteria such as minimum time, energy, or effort. The Pontryagin’s Maximum Principle [59,60] introduces a co-state $\lambda$ and defines the control Hamiltonian as $H(x,u,t)=\lambda \left(t\right)·a(x,u,t)-L(x,u,t)$. The necessary conditions for an extremum are:

$$\left\{\begin{array}{cc}\hfill 0& ={\displaystyle \frac{\partial H}{\partial u}}(\mathrm{Control}\mathrm{Optimality})\hfill \\ \hfill \dot{x}& ={\displaystyle \frac{\partial H}{\partial \lambda }}(\mathrm{State}\mathrm{Equation})\hfill \\ \hfill \dot{\lambda }& =-{\displaystyle \frac{\partial H}{\partial x}}(\text{Co-state}/\mathrm{Adjoint}\mathrm{Equation})\hfill \end{array}\right.$$

(A2)

The last two equations are identical in their form to the canonical Hamilton’s equations, where x serves as the coordinate and $\lambda$ as the momentum.

In the present study, we assume the control Lagrangian L is unknown, rendering the exact Hamiltonian in (A2) inaccessible. Furthermore, we assume no prior knowledge regarding the specific nature of the control $u\left(t\right)$—including its dimensionality, application axes, or functional form. However, we do assume that the symmetry of the configuration manifold is known, and the resulting Pontryagin’s Hamiltonian possesses the required symmetry criteria. This assumption is physically and mathematically justified, as the space in which the particles evolve is usually geometrically well-defined. Our objective is to forecast the system’s trajectory by synthesizing this known geometric symmetry with observational data of its past evolution. To this end, we apply the symmetry reduction theory for optimal control as formulated in [1,2].

Symmetry Reduction Applied to Optimal Control

The Lie–Poisson Reduction Theorem [45] (see the proof for $\mathrm{GL}\left(n\right)$ presented in Section 13.2 of that paper), derives the Lie–Poisson bracket from the general canonical bracket through a process known as symmetry reduction. In classical mechanics, this procedure is utilized to map the evolution of a system from the full phase space, which is typically the cotangent bundle of a configuration manifold, to the Lie algebra. Since that the phase space is the cotangent bundle of the full Lie group (all coordinates and momenta), the reduction to the Lie algebra of that group, which is a linear space, is a drastic simplification. A practical example is the simulation of satellite dynamics via Euler’s equations; the motion is calculated within the space of angular momenta rather than the more complex space of orientation matrices and their associated momenta.

Symmetry reduction in optimal control mirrors the Lie-Poisson reduction for in canonical Hamiltonian systems [61,62,63]. As described in [41] and [43] (pp. 227-267), this process reduces the optimality conditions in (A2) to a Lie–Poisson system when Lie group symmetries are present. For an optimal control problem involving a single agent, the state variable x in (A2) is assumed to reside on a Lie group G. In this geometric context, the co-state $\lambda \in T^{*}G$ acts as the mathematical analogue of mechanical momentum. If the Pontryagin Lagrangian is determined solely by the control u, and that control directly influences velocities, as is the case in (A10) below, the resulting Hamiltonian H permits a complete symmetry reduction to a Lie–Poisson system.

This theoretical framework extends to the collective dynamics of multiple interacting controlled agents. Drawing again on the mechanical analogy, we begin with a canonical Hamiltonian $H(g_1,\dots ,g_N,p_1,\dots ,p_N)$, where each configuration $g_i$ belongs to G and each momentum $p_i$ resides in the cotangent space $T_{g_i}^{*}G$. If the system exhibits symmetry only with respect to the group G, rather than the full product manifold $G^N$, a complete Lie-Poisson reduction is generally not possible. Instead, the reduced space must include additional Lie-group-valued variables representing relative configurations, defined as ${\alpha }_{ij}=g_i^{-1}{g}_j$. The resulting governing equations then consist of Lie-Poisson components coupled with the evolution of these relative elements ${\alpha }_{ij}$ [38], introducing significant analytical complexity.

To illustrate this point, consider $G=SO\left(3\right)$ describing the evolution of coupled rigid bodies. The global rigid body symmetry in a coupled two-body system [61,64] allows for a reduction from the initial phase space which is a cotangent bundle of $G\times G$, i.e.,

$$\begin{array}{c}{T}^{*}G\times T^{*}G\cong G\times G\times \mathfrak{g}\times \mathfrak{g}\end{array}$$

to the reduced manifold $G\times \mathfrak{g}\times \mathfrak{g}$. Physically, this reduction captures the coupled evolution of the two momenta alongside the relative orientation of the objects. Consequently, in mechanical problems involving $G^N$ manifolds with G-symmetry, complete reduction to ${\left({\mathfrak{g}}^{*}\right)}^N$ is rare unless the Hamiltonian possesses a highly specific structure.

Similar constraints would, in priniple, apply to the optimal control of interacting particles [1,2]. However, due to the unique functional form of the control Hamiltonian, these systems actually achieve a complete symmetry reduction to the pure Lie-Poisson case. This useful fact bypasses the need to explicitly track relative orientations.

The following section provides the mathematical formulation for the collective dynamics of N controlled particles. This foundation is necessary for examining the ’Dictatorship’ and ’Democracy’ control Hamiltonian models discussed in Section A.2. Although the primary focus of this paper is a data-driven methodology, a rigorous exposition of this general theory is essential to illustrate how the ground-truth data used for learning is generated.

The General Theory of Controlled Dynamics of N Interacting Particles

In the following derivation, we adhere to the notation established in [1,2]. We model the connectivity between particles using a connected, undirected graph where the vertex set is $V=\{v_1,v_2,\dots ,v_N\}$. The graph contains no self-loops, and $d\left(v_i\right)$ represents the degree of vertex $v_i$. We begin with several foundational definitions.

Definition A1 

(Definition of graph quantities).

1. 

Let A be the adjacency matrix with elements $a_{ij}$ for $i,j=1,\dots ,N$. The entries are defined such that $a_{ij}=1$ if a connection exists between vertices $v_i$ and $v_j$, and $a_{ij}=0$ otherwise.

2. 

The degree matrix D is a diagonal matrix defined as

$D=diag(d\left(v_1\right),\dots ,d\left(v_N\right))$.

3. 

The graph Laplacian is given by $B=D-A$, with elements denoted by $b_{ij}$ for $i,j=1,\dots ,N$.

By construction, the matrices A, D, and B are symmetric [42]. Furthermore, the Laplacian satisfies $B{\mathbf{1}}_N=\mathbf{0}$, where ${\mathbf{1}}_N={[11\dots 1]}^{\mathsf{T}}$ is the $N\times 1$ vector of ones. Consequently, ${\mathbf{1}}_N$ serves as an eigenvector of B corresponding to the zero eigenvalue; physically, this implies that the row and column sums of B are identically zero.

Consider an ensemble of N agents interacting via a control protocol modeled by this graph. In the framework used here, the particles (vehicles) correspond to vertices, and edges represent active communication or physical coupling. We assume each agent evolves on a finite-dimensional Lie group G; therefore, the collective system state resides in $G^N$, the direct product of N copies of G. Let n denote the dimension of the Lie group G and its associated Lie algebra $\mathfrak{g}$. We define the inner product $\langle \langle \xi ,\eta \rangle \rangle =Tr\left({\xi }^{\mathsf{T}}\eta \right)$ for $\xi ,\eta \in \mathfrak{g}$. This inner product defines the trace norm ${\parallel \xi \parallel }^2=\langle \langle \xi ,\xi \rangle \rangle =Tr\left({\xi }^{\mathsf{T}}\xi \right)$. We select a basis $\{X_1,X_2,\dots ,X_n\}$ for the Lie algebra, assumed to be orthonormal under the trace inner product. Additionally, let $\langle ·,·\rangle$ represent the duality pairing between the primal and dual spaces [2].

The derivation of symmetry reduction for Pontryagin’s maximum principle in this context was originally presented in [1,2,41], which we reproduce succinctly here. We focus on the specific case where the state $x=g\in G$. We examine the optimal control problem (A1) where the Lagrangian $L\left(u\right)$ depends exclusively on the controls, leading to the minimization objective:

$$\begin{array}{c}{\mathrm{Min}}_u{\int }_0^{T}L\left(u\right)dt.\end{array}$$

(A3)

Since $L\left(u\right)$ is independent of the group configuration g, it is a G-invariant Lagrangian. Each variable $u\left(t\right)$ generates the control ${\xi }_{u\left(t\right)}=:\xi$ within the Lie algebra $\mathfrak{g}$, governed by the kinematic equation:

$$\dot{g}=g\xi ,$$

(A4)

subject to fixed boundary conditions $g\left(0\right)=g_0$ and $g\left(T\right)=g_1$. The canonical state and co-state equations (the last two equations in (A2)) are superseded by equations derived from the Lie–Poisson bracket via the Lie–Poisson Reduction Theorem (see, e.g., [45], Section 13.2). The link between controls and symmetry-reduced momenta corresponds to the first condition of (A2). We define the pre-Hamiltonian as:

$$\begin{array}{c}H(p,g,u)=\langle p,g\xi \rangle -L\left(u\right)=\langle \mu ,\xi \rangle -L\left(u\right),\end{array}$$

(A5)

where $p\in T_g^{*}{G}^N$ is the momentum in the cotangent space, and $\mu \in {\left({\mathfrak{g}}^{*}\right)}^N$ resides in the dual of the collective Lie algebra. While $\mu$ should be formally referred to as the ’symmetry-reduced co-state’, we will adopt the mechanical analogy and refer to it simply as ’momentum’.

For simplicity of exposition, we restrict our analysis to Lie groups admitting a matrix representation. This is not really a limitation of the method since all Lie groups likely to be encountered in practical problems can be represented in terms of matrices. Consequently, both the group elements g and the algebra elements $\xi$ are expressed as block-diagonal matrices:

$$\begin{array}{c}g=diag(g_1,g_2,\dots ,g_N),\xi =diag({\xi }_1,{\xi }_2,\dots ,{\xi }_N),\end{array}$$

where the aggregate control vector is $u=(u_1,u_2,\dots ,u_N)$ with total dimension $mN$. In this formulation, the pre-Hamiltonian (A5) becomes:

$$\begin{array}{c}H(p,g,u)={\displaystyle \sum _{k=1}^N}\langle p_k,g_k{\xi }_k\rangle -L\left(u\right)={\displaystyle \sum _{k=1}^N}\langle {\mu }_k,{\xi }_k\rangle -L\left(u\right).\end{array}$$

(A6)

The condition for optimality with respect to the control variables u i.e., the first extremal condition in (A2), becomes:

$${\displaystyle \frac{\partial H(\mu ,u)}{\partial u}}=0\iff 〈\mu ,{\displaystyle \frac{\partial \xi }{\partial u}}〉-{\displaystyle \frac{\partial L}{\partial u}}=0,$$

(A7)

where $\mu =({\mu }_1,\dots ,{\mu }_N)$ is the composite momentum vector. By the implicit function theorem [65], one can solve (A7) for $u\left(\mu \right)$ provided the Hessian ${Hess}_{u}H(\mu ,u)$ is non-degenerate. If $\xi$ is linear in u, this non-degeneracy condition reduces to the requirement that the Hessian of the Lagrangian, $det{Hess}_{u}L\left(u\right)\ne 0$, which we assume holds true.

Adopting the framework of [1,2], we consider the Lagrangian specific to interacting particles:

$$\begin{array}{c}L=L({\xi }_1,\dots ,{\xi }_N)={\displaystyle \frac{1}{2}}\left({\displaystyle \sum _{k=1}^N}|{\xi }_k{|}^2+\chi {\displaystyle \sum _{k=1}^N}{\displaystyle \sum _{j=1}^N}{a}_{kj}{|{\xi }_k-{\xi }_j|}^2\right),\end{array}$$

(A8)

where $\chi \ge 0$ is a coupling constant. Utilizing the relationship between the adjacency matrix A and the Laplacian B (where diagonal elements $b_{ii}=\sum a_{ij}$ and off-diagonals $b_{ij}=-a_{ij}$), we can restructure the Lagrangian (A8) as:

$$\begin{array}{c}L=L({\xi }_1,\dots ,{\xi }_N)={\displaystyle \frac{1}{2}}{\displaystyle \sum _{k=1}^N}{\left|{\xi }_k\right|}^2+\chi {\displaystyle \sum _{k=1}^N}〈〈{\xi }_k,{\displaystyle \sum _{j=1}^N}{b}_{kj}{\xi }_j〉〉.\end{array}$$

(A9)

We assume the control velocities ${\xi }_k$ are affine functions of the control inputs $u_k$:

$$\begin{array}{c}\hfill {\xi }_k=X_q+{\displaystyle \sum _{i=1}^m}{u}_{ki}{X}_i,k=1,\dots ,N,\end{array}$$

(A10)

Here, $u_k\in {\mathbb{R}}^m$, and the basis $\{X_1,\dots ,X_n\}$ is orthonormal. The system is underactuated ($m<n$) with a drift component $X_q$ aligned with an uncontrolled direction ($m+1\le q\le n$). While the precise form of (A10) is treated as unknown during the discovery phase, the knowledge of its’ particular form is essential for generating the ground truth data.

We introduce the dual basis for the coalgebra ${\mathfrak{g}}^{*}$, denoted $\{X_1^b,\dots ,X_n^b\}$. Letting ${\mu }_{ki}$ be the expansion coefficients such that ${\mu }_k={\sum }_{i=1}^n{\mu }_{ki}{X}_i^b$, the optimality condition (A7) simplifies to:

$$\begin{array}{c}〈{\mu }_k,{\displaystyle \frac{\partial {\xi }_k}{\partial u_{ki}}}〉-{\displaystyle \frac{\partial L}{\partial u_{ki}}}=0\iff 〈{\mu }_k,X_i〉-{\displaystyle \frac{\partial L}{\partial u_{ki}}}=0\iff {\mu }_{ki}-{\displaystyle \frac{\partial L}{\partial u_{ki}}}=0\end{array}$$

(A11)

for $k=1,\dots ,N$ and $i=1,\dots ,m$. Differentiating the Lagrangian (A8) establishes the relationship between momenta and controls:

$${\displaystyle \frac{\partial L}{\partial u_{ki}}}=u_{ki}+\chi {\displaystyle \sum _{j=1}^N}{b}_{jk}{u}_{ji}={\mu }_{ki},k=1,\dots ,N,i=1,\dots ,m.$$

(A12)

To vectorize this, let $\breve{\mu }={\left[\begin{array}{ccc}{\mu }_1^{\mathsf{T}}& \dots & {\mu }_N^{\mathsf{T}}\end{array}\right]}^{\mathsf{T}}$ be the full momentum vector of length $Nn$. We denote the vector of the m controlled components for the k-th agent as ${\tilde{\mu }}_k=({\mu }_{k1},\dots ,{\mu }_{km})$. The relationship between the controlled momentum components and control inputs is then:

$${\tilde{\mu }}_k=u_k+2\chi {\displaystyle \sum _{j=1}^N}{b}_{jk}{u}_{ji}\iff \left[\begin{array}{c}{\tilde{\mu }}_1\\ ⋮\\ {\tilde{\mu }}_N\end{array}\right]=(({\mathbb{I}}_N+2\chi B)\otimes {\mathbb{I}}_m)\left[\begin{array}{c}{u}_1\\ ⋮\\ u_N\end{array}\right],$$

(A13)

where ⊗ is the Kronecker product. Since the eigenvalues of B are non-negative, the matrix inverse exists for $\chi \ge 0$. We define the coupling matrix $\mathsf{\Psi }$ as:

$$\begin{array}{c}\mathsf{\Psi }={(({\mathbb{I}}_N+2\chi B)\otimes {\mathbb{I}}_m)}^{-1}={({\mathbb{I}}_N+2\chi B)}^{-1}\otimes {\mathbb{I}}_m.\end{array}$$

(A14)

Solving for the controls yields $u=\mathsf{\Psi }\tilde{\mu }$. Substituting these optimal controls back into the pre-Hamiltonian (A6) yields the reduced control Hamiltonian (up to a constant) [1,2]:

$$h={\displaystyle \sum _{k=1}^N}{\mu }_{kq}+{\displaystyle \frac{1}{2}}\left[\begin{array}{ccc}{\tilde{\mu }}_1& \dots & {\tilde{\mu }}_N\end{array}\right]\mathsf{\Psi }\left[\begin{array}{c}{\tilde{\mu }}_1\\ ⋮\\ {\tilde{\mu }}_N\end{array}\right].$$

(A15)

Crucially, this reduced system evolves entirely on ${\left({\mathfrak{g}}^{*}\right)}^N$, avoiding the need for variables describing relative orientation. The collective dynamics for $\breve{\mu }$ are governed by the Lie–Poisson equations:

$$\dot{\breve{\mu }}=\mathsf{\Lambda }\left(\breve{\mu }\right)\nabla h,\mathsf{\Lambda }\left(\breve{\mu }\right)=-{\displaystyle \frac{1}{\sqrt{2}}}\left[\begin{array}{cccc}{\widehat{\mu }}_1& 0& \dots & 0\\ 0& {\widehat{\mu }}_2& \dots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \dots & {\widehat{\mu }}_N\end{array}\right],$$

(A16)

where the Poisson tensor block ${\widehat{\mu }}_k$ is an $n\times n$ antisymmetric matrix defined by:

$$\begin{array}{c}{\widehat{\mu }}_k=-{\displaystyle \sum _{s=1}^n}{\mu }_{ks}\left[\begin{array}{cccc}{\mathsf{\Gamma }}_{11}^s& {\mathsf{\Gamma }}_{12}^s& \dots & {\mathsf{\Gamma }}_{1n}^s\\ {\mathsf{\Gamma }}_{21}^s& {\mathsf{\Gamma }}_{22}^s& \dots & {\mathsf{\Gamma }}_{2n}^s\\ ⋮& ⋮& \ddots & ⋮\\ {\mathsf{\Gamma }}_{n1}^s& {\mathsf{\Gamma }}_{n2}^s& \dots & {\mathsf{\Gamma }}_{nn}^s\end{array}\right].\end{array}$$

(A17)

Here, ${\mathsf{\Gamma }}_{ij}^s$ are the structure constants satisfying $[X_i,X_j]={\sum }_{s=1}^n{\mathsf{\Gamma }}_{ij}^s{X}_s$.

Remark A1 

(Geometric structure). As a result of the complete symmetry reduction, equations (A16) can be expressed in coordinate-free notation as ${\dot{\mu }}_k=-{\mathrm{ad}}_{\partial h/\partial {\mu }_k}^{*}{\mu }_k$, where ${\mathrm{ad}}^{*}$ denotes the co-adjoint action of the Lie algebra $\mathfrak{g}$ (assuming left invariance [44,45]).

Remark A2. 

The structure of the Poisson tensor implies that $\breve{\mu }=0$ remains a stationary point for any control Hamiltonian h.

To further detail the reduced Hamiltonian $h\left(\breve{\mu }\right)$, we must next consider specific graph Laplacians B and particle dynamics on specific Lie groups.

Appendix A.2. Two Different Cases of Adjancency Matrix

In what follows, to be concrete, we consider two cases of an adjacency matrix. The first case refers to the adjacency matrix A when one ’leading’ particle is connected to all other particles, which we call ’Dictatorship’. Another case happens when the adjacency matrix A is such that all particles are connected with all other particles, which we will denote as ’Democracy’. In these two cases, we can explicitly compute the inverse matrix ${({\mathbb{I}}_N+2\chi B)}^{-1}$ in the corresponding part of the matrix $\mathsf{\Psi }$ defined by (A14) in its explicit form. The matrix $({\mathbb{I}}_N+2\chi B)$ is always invertible, but the explicit formula in (A18) is useful for analytic considerations.

The Case of ’Dictatorship’ Governance

At first, we consider a ’Dictatorship’ case, that is, when particle 1 is connected to all other particles. There are no other connections between particles (edges of the graph). In this case the adjacency matrix A is given by the $N\times N$ matrix

$$\begin{array}{c}A=\left[\begin{array}{cccccc}0& 1& 1& \dots & 1& 1\\ 1& 0& 0& \dots & 0& 0\\ 1& 0& 0& \dots & 0& 0\\ ⋮& ⋮& ⋮& \ddots & ⋮& ⋮\\ 1& 0& 0& \dots & 0& 0\\ 1& 0& 0& \dots & 0& 0\end{array}\right],\end{array}$$

The degree matrix is the $N\times N$ matrix $D=diag(N-1,1,1,\dots ,1,1)$, and the graph Laplacian is the $N\times N$ matrix

$$\begin{array}{c}B=D-A=\left[\begin{array}{cccccc}N-1& -1& -1& \dots & -1& -1\\ -1& 1& 0& \dots & 0& 0\\ -1& 0& 1& \dots & 0& 0\\ ⋮& ⋮& ⋮& \ddots & ⋮& ⋮\\ -1& 0& 0& \dots & 1& 0\\ -1& 0& 0& \dots & 0& 1\end{array}\right].\end{array}$$

Moreover, for this particular case of ’Dictatorship’ control, we can, after some algebra, obtain the exact formula for the inverse of $({\mathbb{I}}_N+2\chi B)$ as:

$$\begin{array}{c}{({\mathbb{I}}_N+2\chi B)}^{-1}={\left(\left[\begin{array}{cccccc}1+2\chi (N-1)& -2\chi & -2\chi & \dots & -2\chi & -2\chi \\ -2\chi & 1+2\chi & 0& \dots & 0& 0\\ -2\chi & 0& 1+2\chi & \dots & 0& 0\\ ⋮& ⋮& ⋮& \ddots & ⋮& ⋮& \\ -2\chi & 0& 0& \dots & 1+2\chi & 0\\ -2\chi & 0& 0& \dots & 0& 1+2\chi \end{array}\right]\right)}^{-1}\\ =\left[\begin{array}{cccccc}{\displaystyle \frac{1+2\chi }{1+2N\chi }}& {\displaystyle \frac{2\chi }{1+2N\chi }}& {\displaystyle \frac{2\chi }{1+2N\chi }}& \dots & {\displaystyle \frac{2\chi }{1+2N\chi }}& {\displaystyle \frac{2\chi }{1+2N\chi }}\\ {\displaystyle \frac{2\chi }{1+2N\chi }}& {\displaystyle \frac{1+2N\chi +4{\chi }^2}{(1+2N\chi )(1+2\chi )}}& {\displaystyle \frac{4{\chi }^2}{(1+2N\chi )(1+2\chi )}}& \dots & {\displaystyle \frac{4{\chi }^2}{(1+2N\chi )(1+2\chi )}}& {\displaystyle \frac{4{\chi }^2}{(1+2N\chi )(1+2\chi )}}\\ {\displaystyle \frac{2\chi }{1+2N\chi }}& {\displaystyle \frac{4{\chi }^2}{(1+2N\chi )(1+2\chi )}}& {\displaystyle \frac{1+2N\chi +4{\chi }^2}{(1+2N\chi )(1+2\chi )}}& \dots & {\displaystyle \frac{4{\chi }^2}{(1+2N\chi )(1+2\chi )}}& {\displaystyle \frac{4{\chi }^2}{(1+2N\chi )(1+2\chi )}}\\ ⋮& ⋮& ⋮& \ddots & ⋮& ⋮\\ {\displaystyle \frac{2\chi }{1+2N\chi }}& {\displaystyle \frac{4{\chi }^2}{(1+2N\chi )(1+2\chi )}}& {\displaystyle \frac{4{\chi }^2}{(1+2N\chi )(1+2\chi )}}& \dots & {\displaystyle \frac{1+2N\chi +4{\chi }^2}{(1+2N\chi )(1+2\chi )}}& {\displaystyle \frac{4{\chi }^2}{(1+2N\chi )(1+2\chi )}}\\ {\displaystyle \frac{2\chi }{1+2N\chi }}& {\displaystyle \frac{4{\chi }^2}{(1+2N\chi )(1+2\chi )}}& {\displaystyle \frac{4{\chi }^2}{(1+2N\chi )(1+2\chi )}}& \dots & {\displaystyle \frac{4{\chi }^2}{(1+2N\chi )(1+2\chi )}}& {\displaystyle \frac{1+2N\chi +4{\chi }^2}{(1+2N\chi )(1+2\chi )}}\end{array}\right].\end{array}$$

(A18)

The Case of ‘Democratic’ Control

Next, we consider a ‘Democracy’ case, that is, when each particle is connected to all of the others. In other words, we have a complete graph. In this case, the adjacency matrix A is given by the $N\times N$ matrix

$$\begin{array}{c}A=\left[\begin{array}{cccccc}0& 1& 1& \dots & 1& 1\\ 1& 0& 1& \dots & 1& 1\\ 1& 1& 0& \dots & 1& 1\\ ⋮& ⋮& ⋮& \ddots & ⋮& ⋮\\ 1& 1& 1& \dots & 0& 1\\ 1& 1& 1& \dots & 1& 0\end{array}\right],\end{array}$$

the degree matrix is the $N\times N$ matrix $D=diag(N-1,N-1,N-1,\dots ,N-1,N-1)$, and the graph Laplacian is the $N\times N$ matrix

$$\begin{array}{c}B=D-A=\left[\begin{array}{cccccc}N-1& -1& -1& \dots & -1& -1\\ -1& N-1& -1& \dots & -1& -1\\ -1& -1& N-1& \dots & -1& -1\\ ⋮& ⋮& ⋮& \ddots & ⋮& ⋮\\ -1& -1& -1& \dots & N-1& -1\\ -1& -1& -1& \dots & -1& N-1\end{array}\right].\end{array}$$

In this particular case, the inverse of the matrix $({\mathbb{I}}_N+2\chi B)$ is also given explicitly as

$$\begin{array}{c}{({\mathbb{I}}_N+2\chi B)}^{-1}=\left[\begin{array}{cccccc}{\displaystyle \frac{1+2\chi }{1+2N\chi }}& {\displaystyle \frac{2\chi }{1+2N\chi }}& {\displaystyle \frac{2\chi }{1+2N\chi }}& \dots & {\displaystyle \frac{2\chi }{1+2N\chi }}& {\displaystyle \frac{2\chi }{1+2N\chi }}\\ {\displaystyle \frac{2\chi }{1+2N\chi }}& {\displaystyle \frac{1+2\chi }{1+2N\chi }}& {\displaystyle \frac{2\chi }{1+2N\chi }}& \dots & {\displaystyle \frac{2\chi }{1+2N\chi }}& {\displaystyle \frac{2\chi }{1+2N\chi }}\\ {\displaystyle \frac{2\chi }{1+2N\chi }}& {\displaystyle \frac{2\chi }{1+2N\chi }}& {\displaystyle \frac{1+2\chi }{1+2N\chi }}& \dots & {\displaystyle \frac{2\chi }{1+2N\chi }}& {\displaystyle \frac{2\chi }{1+2N\chi }}\\ ⋮& ⋮& ⋮& \ddots & ⋮& ⋮\\ {\displaystyle \frac{2\chi }{1+2N\chi }}& {\displaystyle \frac{2\chi }{1+2N\chi }}& {\displaystyle \frac{2\chi }{1+2N\chi }}& \dots & {\displaystyle \frac{1+2\chi }{1+2N\chi }}& {\displaystyle \frac{2\chi }{1+2N\chi }}\\ {\displaystyle \frac{2\chi }{1+2N\chi }}& {\displaystyle \frac{2\chi }{1+2N\chi }}& {\displaystyle \frac{2\chi }{1+2N\chi }}& \dots & {\displaystyle \frac{2\chi }{1+2N\chi }}& {\displaystyle \frac{1+2\chi }{1+2N\chi }}\end{array}\right].\end{array}$$

(A19)

The inversion of the matrix $({\mathbb{I}}_N+2\chi B)$ does not allow explicit formulas as in the case of ’Dictatorship’ (A18) or ’Democracy’ (A19). Thus, it is advantageous to consider the dynamics for these particular cases of interactions between the particles, as they allow additional analytical insights into the problem.

Analytical solutions may be obtained for a single particle dynamics for particular groups in terms of elliptic functions. Such solution for a particle evolving on $\mathrm{SE}\left(2\right)$ group was obtained in [2]. In the Appendix, we derive these analytical solutions for $\mathrm{SO}\left(3\right)$ group in Section B.2 and a particular realization of the control in $\mathrm{SE}\left(3\right)$ group in Section B.3.0.1. While the single-particle solution does not address the collective behavior, we present these solutions for reference since they allow analytical solutions, and, to our knowledge, have not been described before.

Appendix B.1. Single Particle Equations in the General Case

This section computes the equations of motion for the single particle, $N=1$. It is clearly not a coupled case, but it also has its own pedagogical merit. In this case, the Poisson tensor will comprise only one block and, thus, will simplify (A17) for a single particle:

$$\begin{array}{c}\hfill \mathsf{\Lambda }=\widehat{\mu }=-{\displaystyle \sum _{s=1}^n}{\mu }_{1s}\left[\begin{array}{cccc}{\mathsf{\Gamma }}_{11}^s& {\mathsf{\Gamma }}_{12}^s& \dots & {\mathsf{\Gamma }}_{1n}^s\\ {\mathsf{\Gamma }}_{21}^m& {\mathsf{\Gamma }}_{22}^s& \dots & {\mathsf{\Gamma }}_{2n}^s\\ ⋮& ⋮& \ddots & ⋮\\ {\mathsf{\Gamma }}_{n1}^s& {\mathsf{\Gamma }}_{n2}^s& \dots & {\mathsf{\Gamma }}_{nn}^s\end{array}\right],\end{array}$$

(A20)

where ${\mathsf{\Gamma }}_{ij}^s$, again, denote structure constants for the corresponding Lie algebra $\mathfrak{g}$. Moreover, now the matrix ${\mathbb{I}}_N={\mathbb{I}}_1=1$, and we can take $B=0$. Then $({\mathbb{I}}_N+2\chi B)=1$, ${({\mathbb{I}}_N+2\chi B)}^{-1}=1$. Finally,

$$\begin{array}{c}\mathsf{\Psi }={({\mathbb{I}}_N+2\chi B)}^{-1}\otimes {\mathbb{I}}_m=1\otimes {\mathbb{I}}_m={\mathbb{I}}_m.\end{array}$$

(A21)

According to (A15), the Hamiltonian looks as follows, since the particle index always equals $k=N=1$:

$$\begin{array}{c}\hfill h={\mu }_{1q}+{\displaystyle \frac{1}{2}}{\displaystyle \sum _{j=1}^m}{\mu }_{1j}^2,\end{array}$$

(A22)

where m is the number of controls. Reference [2] considered the case of $\mathrm{SE}\left(2\right)$ group; we extend these considerations for $\mathrm{SO}\left(3\right)$ and $\mathrm{SE}\left(3\right)$ groups.

Appendix B.2. The Case of Single Particle Dynamics Evolving on SO (3)

In this Section, we consider the Lie–Poisson reduced dynamics for a single particle evolving on $\mathrm{SO}\left(3\right)$ group. We use $m=1$ (one control in the first element) and $q=2$ (drift in the second element). With the use of (A17) and in terms of notation (26) for k=1, the Poisson tensor is

$$\begin{array}{c}\mathsf{\Lambda }={\displaystyle \frac{1}{\sqrt{2}}}\left[\begin{array}{ccc}0& -{\mu }_3& {\mu }_2\\ {\mu }_3& 0& -{\mu }_1\\ -{\mu }_2& {\mu }_1& 0\end{array}\right]={\displaystyle \frac{1}{\sqrt{2}}}\widehat{\mu }.\end{array}$$

According to (A22), the Hamiltonian becomes

$$h={\mu }_2+{\displaystyle \frac{1}{2}}{\mu }_1^2.$$

(A23)

The Casimir function is $c={\mu }_1^2+{\mu }_2^2+{\mu }_3^2={\left|\right|\mu \left|\right|}^2$ which clearly satisfies $\mathsf{\Lambda }\nabla c=0$.

For this single particle problem, the Lie–Poisson reduced dynamics equations $\dot{\breve{\mu }}=\mathsf{\Lambda }\left(\breve{\mu }\right)\nabla h$ (A16), where $\nabla h={({\mu }_1,1,0)}^T$, are given by

$$\begin{array}{c}\hfill {\displaystyle \frac{d}{dt}}\left[\begin{array}{c}{\mu }_1\\ {\mu }_2\\ {\mu }_3\end{array}\right]={\displaystyle \frac{1}{\sqrt{2}}}\left[\begin{array}{ccc}0& -{\mu }_3& {\mu }_2\\ {\mu }_3& 0& -{\mu }_1\\ -{\mu }_2& {\mu }_1& 0\end{array}\right]\left[\begin{array}{c}{\mu }_1\\ 1\\ 0\end{array}\right]\end{array}$$

or, equivalently,

$$\begin{array}{c}\hfill \left\{\begin{array}{c}{\dot{\mu }}_1=-{\displaystyle \frac{1}{\sqrt{2}}}{\mu }_3,\hfill \\ {\dot{\mu }}_2={\displaystyle \frac{1}{\sqrt{2}}}{\mu }_1{\mu }_3,\hfill \\ {\dot{\mu }}_3=-{\displaystyle \frac{1}{\sqrt{2}}}{\mu }_1{\mu }_2+{\displaystyle \frac{1}{\sqrt{2}}}{\mu }_1.\hfill \end{array}\right.\end{array}$$

(A24)

One verifies that h is a conserved quantity, as we have

$$\begin{array}{c}\hfill \dot{h}={\dot{\mu }}_2+{\mu }_1{\dot{\mu }}_1={\displaystyle \frac{1}{\sqrt{2}}}{\mu }_1{\mu }_3-{\mu }_1{\displaystyle \frac{1}{\sqrt{2}}}{\mu }_3=0.\end{array}$$

Hence, the Hamiltonian and the Casimir function are conserved quantities. From (A24) we obtain the second-order equation

$$\begin{array}{c}\hfill {\ddot{\mu }}_1=-{\displaystyle \frac{1}{\sqrt{2}}}{\dot{\mu }}_3={\displaystyle \frac{1}{2}}{\mu }_1({\mu }_2-1)={\displaystyle \frac{1}{2}}{\mu }_1\left((h-1)-{\displaystyle \frac{1}{2}}{\mu }_1^2\right),\end{array}$$

that is,

$${\ddot{\mu }}_1-{\displaystyle \frac{(h-1)}{2}}{\mu }_1+{\displaystyle \frac{1}{4}}{\mu }_1^3=0.$$

(A25)

In fact, we observe that (A25) is the same as Equation (36) in [1] and eq. (3.24) in [2] (which, we remind the reader, was written for $\mathrm{SE}\left(2\right)$ group), with the only difference $h\to (h-1)$. Equation (A25) can be solved in terms of Jacobi elliptic functions. Interestingly enough, (A25) also appears in the theory of elasticity [66].

Appendix B.3. The Case of SE (3) with a Particular Hamiltonian Allowing Analytical Solution

Lie–Poisson Reduced Dynamics for a Single Particle

In this Section, we will introduce the Lie–Poisson reduced dynamics for a single particle, evolving on $\mathrm{SE}\left(3\right)$ group. We use controls in the first two coordinates (rotations) and drift in 4-th coordinate (first linear momentum). Using (A17) and notation (34) for $k=1$, the Poisson tensor is

$$\begin{array}{c}\mathsf{\Lambda }={\displaystyle \frac{1}{\sqrt{2}}}\left[\begin{array}{cccccc}0& -{\mu }_3& {\mu }_2& 0& -{\mu }_6& {\mu }_5\\ {\mu }_3& 0& -{\mu }_1& {\mu }_6& 0& -{\mu }_4\\ -{\mu }_2& {\mu }_1& 0& -{\mu }_5& {\mu }_4& 0\\ 0& -{\mu }_6& {\mu }_5& 0& 0& 0\\ {\mu }_6& 0& -{\mu }_4& 0& 0& 0\\ -{\mu }_5& {\mu }_4& 0& 0& 0& 0\end{array}\right]={\displaystyle \frac{1}{\sqrt{2}}}\left[\begin{array}{cc}\widehat{\mathsf{\Pi }}& \widehat{p}\\ \widehat{p}& ⊬\end{array}\right].\end{array}$$

With the use of (A21):

$$\begin{array}{c}\mathsf{\Psi }={\mathbb{I}}_2=\left[\begin{array}{cc}1& 0\\ 0& 1\end{array}\right].\end{array}$$

According to (A22), the Hamiltonian looks like

$$h={\mu }_4+{\displaystyle \frac{1}{2}}{\mu }_1^2+{\displaystyle \frac{1}{2}}{\mu }_2^2.$$

(A26)

There are two Casimir functions: $c_1={\mu }_4^2+{\mu }_5^2+{\mu }_6^2={\left|\right|\mathit{p}\left|\right|}^2$ and $c_2={\mu }_1{\mu }_4+{\mu }_2{\mu }_5+{\mu }_3{\mu }_6=\mathsf{\Pi }·\mathit{p}$ which clearly satisfy $\mathsf{\Lambda }\nabla c_1=0$ and $\mathsf{\Lambda }\nabla c_2=0$ .

For this single particle problem the Lie–Poisson reduced dynamics equations $\dot{\breve{\mu }}=\mathsf{\Lambda }\left(\breve{\mu }\right)\nabla h$ (A16), where $\nabla h={({\mu }_1,{\mu }_2,0,1,0,0)}^T$, are given by

$$\begin{array}{c}\hfill {\displaystyle \frac{d}{dt}}\left[\begin{array}{c}{\mu }_1\\ {\mu }_2\\ {\mu }_3\\ {\mu }_4\\ {\mu }_5\\ {\mu }_6\end{array}\right]={\displaystyle \frac{1}{\sqrt{2}}}\left[\begin{array}{cccccc}0& -{\mu }_3& {\mu }_2& 0& -{\mu }_6& {\mu }_5\\ {\mu }_3& 0& -{\mu }_1& {\mu }_6& 0& -{\mu }_4\\ -{\mu }_2& {\mu }_1& 0& -{\mu }_5& {\mu }_4& 0\\ 0& -{\mu }_6& {\mu }_5& 0& 0& 0\\ {\mu }_6& 0& -{\mu }_4& 0& 0& 0\\ -{\mu }_5& {\mu }_4& 0& 0& 0& 0\end{array}\right]\left[\begin{array}{c}{\mu }_1\\ {\mu }_2\\ 0\\ 1\\ 0\\ 0\end{array}\right]\end{array}$$

or, equivalently,

$$\begin{array}{c}\hfill \left\{\begin{array}{c}{\dot{\mu }}_1=-{\displaystyle \frac{1}{\sqrt{2}}}{\mu }_2{\mu }_3,\hfill \\ {\dot{\mu }}_2={\displaystyle \frac{1}{\sqrt{2}}}{\mu }_1{\mu }_3+{\displaystyle \frac{1}{\sqrt{2}}}{\mu }_6,\hfill \\ {\dot{\mu }}_3=-{\displaystyle \frac{1}{\sqrt{2}}}{\mu }_5\hfill \\ {\dot{\mu }}_4=-{\displaystyle \frac{1}{\sqrt{2}}}{\mu }_2{\mu }_6,\hfill \\ {\dot{\mu }}_5={\displaystyle \frac{1}{\sqrt{2}}}{\mu }_1{\mu }_6,\hfill \\ {\dot{\mu }}_6=-{\displaystyle \frac{1}{\sqrt{2}}}{\mu }_1{\mu }_5+{\displaystyle \frac{1}{\sqrt{2}}}{\mu }_2{\mu }_4.\hfill \end{array}\right.\end{array}$$

(A27)

One verifies that h is a conserved quantity, as we have

$$\begin{array}{cc}\hfill \dot{h}=& {\dot{\mu }}_4+{\mu }_1{\dot{\mu }}_1+{\mu }_2{\dot{\mu }}_2=\hfill \\ \hfill & =-{\displaystyle \frac{1}{\sqrt{2}}}{\mu }_2{\mu }_6+{\mu }_1\left(-{\displaystyle \frac{1}{\sqrt{2}}}{\mu }_2{\mu }_3\right)+{\mu }_2\left({\displaystyle \frac{1}{\sqrt{2}}}{\mu }_1{\mu }_3+{\displaystyle \frac{1}{\sqrt{2}}}{\mu }_6\right)=0.\hfill \end{array}$$

Hence, the Hamiltonian and both Casimir functions are conserved quantities. No further simplification of equations is possible.

However, if we were to consider a slightly different Hamiltonian with $q=6$ (drift in the 3rd linear momentum coordinate), so that $h={\mu }_6+{\displaystyle \frac{1}{2}}{\mu }_1^2+{\displaystyle \frac{1}{2}}{\mu }_2^2$, then Lie-Poisson equations (A16) can be simplified and solved exactly for particular initial conditions. For this single particle problem the Lie–Poisson reduced dynamics equations $\dot{\breve{\mu }}=\mathsf{\Lambda }\left(\breve{\mu }\right)\nabla h$ (A16), where $\nabla h={({\mu }_1,{\mu }_2,0,0,0,1)}^T$ are given by

$$\begin{array}{c}\hfill {\displaystyle \frac{d}{dt}}\left[\begin{array}{c}{\mu }_1\\ {\mu }_2\\ {\mu }_3\\ {\mu }_4\\ {\mu }_5\\ {\mu }_6\end{array}\right]={\displaystyle \frac{1}{\sqrt{2}}}\left[\begin{array}{cccccc}0& -{\mu }_3& {\mu }_2& 0& -{\mu }_6& {\mu }_5\\ {\mu }_3& 0& -{\mu }_1& {\mu }_6& 0& -{\mu }_4\\ -{\mu }_2& {\mu }_1& 0& -{\mu }_5& {\mu }_4& 0\\ 0& -{\mu }_6& {\mu }_5& 0& 0& 0\\ {\mu }_6& 0& -{\mu }_4& 0& 0& 0\\ -{\mu }_5& {\mu }_4& 0& 0& 0& 0\end{array}\right]\left[\begin{array}{c}{\mu }_1\\ {\mu }_2\\ 0\\ 0\\ 0\\ 1\end{array}\right]\end{array}$$

or, equivalently,

$$\begin{array}{c}\hfill \left\{\begin{array}{c}{\dot{\mu }}_1=-{\displaystyle \frac{1}{\sqrt{2}}}{\mu }_2{\mu }_3+{\displaystyle \frac{1}{\sqrt{2}}}{\mu }_5,\hfill \\ {\dot{\mu }}_2={\displaystyle \frac{1}{\sqrt{2}}}{\mu }_1{\mu }_3-{\displaystyle \frac{1}{\sqrt{2}}}{\mu }_4,\hfill \\ {\dot{\mu }}_3=0\hfill \\ {\dot{\mu }}_4=-{\displaystyle \frac{1}{\sqrt{2}}}{\mu }_2{\mu }_6,\hfill \\ {\dot{\mu }}_5={\displaystyle \frac{1}{\sqrt{2}}}{\mu }_1{\mu }_6,\hfill \\ {\dot{\mu }}_6=-{\displaystyle \frac{1}{\sqrt{2}}}{\mu }_1{\mu }_5+{\displaystyle \frac{1}{\sqrt{2}}}{\mu }_2{\mu }_4.\hfill \end{array}\right.\end{array}$$

(A28)

One verifies that h is a conserved quantity, as we have

$$\begin{array}{cc}\hfill \dot{h}=& {\dot{\mu }}_6+{\mu }_1{\dot{\mu }}_1+{\mu }_2{\dot{\mu }}_2=\hfill \\ \hfill & =-{\displaystyle \frac{1}{\sqrt{2}}}{\mu }_1{\mu }_5+{\displaystyle \frac{1}{\sqrt{2}}}{\mu }_2{\mu }_4+{\mu }_1\left(-{\displaystyle \frac{1}{\sqrt{2}}}{\mu }_2{\mu }_3+{\displaystyle \frac{1}{\sqrt{2}}}{\mu }_5\right)+{\mu }_2\left({\displaystyle \frac{1}{\sqrt{2}}}{\mu }_1{\mu }_3-{\displaystyle \frac{1}{\sqrt{2}}}{\mu }_4\right)=0.\hfill \end{array}$$

From the third equation of (A28) we obtain that ${\mu }_3=const$. Let us suppose the initial condition ${\mu }_3\left(0\right)=0$, that is, ${\mu }_3\equiv 0$. In this case the system of equations (A27) simpifies as follows:

$$\begin{array}{c}\hfill \left\{\begin{array}{c}{\dot{\mu }}_1={\displaystyle \frac{1}{\sqrt{2}}}{\mu }_5,\hfill \\ {\dot{\mu }}_2=-{\displaystyle \frac{1}{\sqrt{2}}}{\mu }_4,\hfill \\ {\dot{\mu }}_3=0({\mu }_3=0),\hfill \\ {\dot{\mu }}_4=-{\displaystyle \frac{1}{\sqrt{2}}}{\mu }_2{\mu }_6,\hfill \\ {\dot{\mu }}_5={\displaystyle \frac{1}{\sqrt{2}}}{\mu }_1{\mu }_6,\hfill \\ {\dot{\mu }}_6=-{\displaystyle \frac{1}{\sqrt{2}}}{\mu }_2{\mu }_5+{\displaystyle \frac{1}{\sqrt{2}}}{\mu }_2{\mu }_4.\hfill \end{array}\right.\end{array}$$

(A29)

Differentiating the fourth equation of (A29) we obtain the second-order equation

$$\begin{array}{c}\hfill {\ddot{\mu }}_6=-{\displaystyle \frac{1}{2}}{\mu }_6\left({\mu }_1^2+{\mu }_2^2\right)-{\displaystyle \frac{1}{2}}\left({\mu }_4^2+{\mu }_5^2\right).\end{array}$$

Substituting expressions for the Hamiltonian h and the first Casimir $c_1$, we obtain:

$$\begin{array}{c}{\ddot{\mu }}_6=-h{\mu }_6+{\displaystyle \frac{3}{2}}{\mu }_6^2-{\displaystyle \frac{c_1}{2}}.\end{array}$$

(A30)

Notice that ${\displaystyle \frac{1}{2}}{\left({\dot{\mu }}_6^2\right)}^{\prime }={\dot{\mu }}_6{\ddot{\mu }}_6$. Integration of this equality yields

$$\begin{array}{c}{\displaystyle \frac{1}{2}}\left({\dot{\mu }}_6^2\right)=\int {\ddot{\mu }}_6{\dot{\mu }}_{6}dt=\int {\ddot{\mu }}_{6}d{\mu }_6.\end{array}$$

Substituting ${\ddot{\mu }}_6$ from (A30) and performing the integration explicitly, we derive

$$\begin{array}{c}{\displaystyle \frac{1}{2}}{\dot{\mu }}_6^2=-{\displaystyle \frac{h{\mu }_6^2}{2}}+{\displaystyle \frac{{\mu }_6^3}{2}}-{\displaystyle \frac{c_1}{2}}{\mu }_6+D,\end{array}$$

(A31)

where D is an arbitrary constant which appears due to integration.

Our next goal is to perform a linear exchange of variables ${\mu }_6=Au+B$ with such constants A and B (to be determined below) that the equation (A31) will turn into

$$\begin{array}{c}{\dot{u}}^2=4u^3-g_{2}u-g_3\end{array}$$

(A32)

which can be solved using Weierstrass elliptic functions ($g_2$ and $g_3$ are some constants). Performing the above mentioned exchange of variables, the equation (A31) is rewritten as

$$\begin{array}{c}{\dot{u}}^2=A{u}^3+(-h+3B)u^2+\left(-{\displaystyle \frac{2hB}{A}}+{\displaystyle \frac{3B^2}{A}}-{\displaystyle \frac{c_1}{A}}\right)u\end{array}$$

(A33)

$$\begin{array}{c}+\left(-{\displaystyle \frac{h{B}^2}{A^2}}+{\displaystyle \frac{B^3}{A^2}}-{\displaystyle \frac{c_{1}B}{A^2}}+{\displaystyle \frac{2D}{A^2}}\right).\end{array}$$

(A34)

Comparing (A32) and (A34) we see that we must have $A=4$ and $B={\displaystyle \frac{h}{3}}$. That is, the desired exchange of variables is ${\mu }_6=4u+{\displaystyle \frac{h}{3}}$. Then the equation (A34) takes the form

$$\begin{array}{c}{\dot{u}}^2=4u^3-\left({\displaystyle \frac{h^2}{12}}+{\displaystyle \frac{c_1}{4}}\right)u-\left({\displaystyle \frac{h^3}{216}}+{\displaystyle \frac{c_{1}h}{48}}-{\displaystyle \frac{D}{8}}\right).\end{array}$$

(A35)

That is, comparing (A32) and (A35), we obtain constants $g_2={\displaystyle \frac{h^2}{12}}+{\displaystyle \frac{c_1}{4}}$ and $g_3={\displaystyle \frac{h^3}{216}}+{\displaystyle \frac{c_{1}h}{48}}-{\displaystyle \frac{D}{8}}$.

Appendix C.1. SO(3) Group

Dictatorship

In the case of ‘Dictatorship’, the Hamiltonian (A15) becomes

$$\begin{array}{cc}\hfill h=& {\displaystyle \sum _{k=1}^N}{\mu }_{k2}+{\displaystyle \frac{1}{2}}({\displaystyle \frac{1+2\chi }{1+2N\chi }}{\mu }_{11}^2+{\displaystyle \frac{1+2N\chi +4{\chi }^2}{(1+2N\chi )(1+2\chi )}}({\mu }_{21}^2+\dots +{\mu }_{N1}^2)\hfill \\ \hfill & +2{\displaystyle \frac{2\chi }{1+2N\chi }}{\mu }_{11}({\mu }_{21}+\dots +{\mu }_{N1})+2{\displaystyle \frac{4{\chi }^2}{(1+2N\chi )(1+2\chi )}}{\displaystyle \sum _{i,j=2,i<j}^N}{\mu }_{i1}{\mu }_{j1}).\hfill \end{array}$$

We remind the reader that the total gradient is split into the gradients with respect to the momenta ${\mathit{\mu }}_k$ of each particle:

$$\begin{array}{c}\nabla h={\displaystyle \frac{\partial h}{\partial \breve{\mu }}}=\left[\begin{array}{c}\partial h/\partial {\mathit{\mu }}_1\\ \partial h/\partial {\mathit{\mu }}_2\\ ⋮\\ \partial h/\partial {\mathit{\mu }}_N\end{array}\right]=\left[\begin{array}{c}{(\nabla h)}_1\\ {(\nabla h)}_2\\ ⋮\\ {(\nabla h)}_N\end{array}\right]\end{array}$$

with the gradients expressed as

$$\begin{array}{c}{(\nabla h)}_k={\displaystyle \frac{\partial h}{\partial {\mathit{\mu }}_k}}=\left[\begin{array}{c}\partial h/\partial {\mu }_{k1}\\ \partial h/\partial {\mu }_{k2}\\ \partial h/\partial {\mu }_{k3}\end{array}\right]=\left[\begin{array}{c}\partial h/\partial {\mu }_{k1}\\ 1\\ 0\end{array}\right],k=1,\dots ,N.\end{array}$$

The partial derivatives are computed explicitly as

$$\begin{array}{cc}\hfill \partial h/\partial {\mu }_{11}& ={\displaystyle \frac{1+2\chi }{1+2N\chi }}{\mu }_{11}+{\displaystyle \frac{2\chi }{1+2N\chi }}{\displaystyle \sum _{j=2}^N}{\mu }_{j1},\hfill \\ \hfill \partial h/\partial {\mu }_{k1}& ={\displaystyle \frac{2\chi }{1+2N\chi }}{\mu }_{11}+{\displaystyle \frac{1+2N\chi +4{\chi }^2}{(1+2N\chi )(1+2\chi )}}{\mu }_{k1}+{\displaystyle \frac{4{\chi }^2}{(1+2N\chi )(1+2\chi )}}{\displaystyle \sum _{\begin{array}{c}i=2\\ i\ne k\end{array}}^N}{\mu }_{i1}\hfill \end{array}$$

with $k=2,\dots ,N$.

Democracy

Similarly, we compute gradients for each particle as

$$\begin{array}{c}{(\nabla h)}_k={\displaystyle \frac{\partial h}{\partial {\mu }_k}}=\left[\begin{array}{c}\partial h/\partial {\mu }_{k1}\\ \partial h/\partial {\mu }_{k2}\\ \partial h/\partial {\mu }_{k3}\end{array}\right]=\left[\begin{array}{c}\partial h/\partial {\mu }_{k1}\\ 1\\ 0\end{array}\right],k=1,\dots ,N.\end{array}$$

We see that the derivatives with respect to ${\mu }_{k1}$ are computed as

$$\begin{array}{c}\partial h/\partial {\mu }_{k1}={\displaystyle \frac{1+2\chi }{1+2N\chi }}{\mu }_{k1}+{\displaystyle \frac{2\chi }{1+2N\chi }}{\displaystyle \sum _{i=1,i\ne k}^N}{\mu }_{i1}\forall k=1,\dots ,N.\end{array}$$

Appendix C.2. SE(3) Group

Dictatorship

The matrix $\mathsf{\Psi }$ in (A15) can be computed explicitly as

$$\begin{array}{c}\mathsf{\Psi }={({\mathbb{I}}_N+2\chi B)}^{-1}\otimes {\mathbb{I}}_2\\ =\left[\begin{array}{ccccccc}{\displaystyle \frac{1+2\chi }{1+2N\chi }}& 0& {\displaystyle \frac{2\chi }{1+2N\chi }}& 0& \dots & {\displaystyle \frac{2\chi }{1+2N\chi }}& 0\\ 0& {\displaystyle \frac{1+2\chi }{1+2N\chi }}& 0& {\displaystyle \frac{2\chi }{1+2N\chi }}& \dots & 0& {\displaystyle \frac{2\chi }{1+2N\chi }}\\ {\displaystyle \frac{2\chi }{1+2N\chi }}& 0& {\displaystyle \frac{1+2N\chi +4{\chi }^2}{(1+2N\chi )(1+2\chi )}}& 0& \dots & {\displaystyle \frac{4{\chi }^2}{(1+2N\chi )(1+2\chi )}}& 0\\ 0& {\displaystyle \frac{2\chi }{1+2N\chi }}& 0& {\displaystyle \frac{1+2N\chi +4{\chi }^2}{(1+2N\chi )(1+2\chi )}}& \dots & 0& {\displaystyle \frac{4{\chi }^2}{(1+2N\chi )(1+2\chi )}}\\ \dots & \dots & \dots & \dots & \dots & \dots & \dots \\ {\displaystyle \frac{2\chi }{1+2N\chi }}& 0& {\displaystyle \frac{4{\chi }^2}{(1+2N\chi )(1+2\chi )}}& 0& \dots & {\displaystyle \frac{4{\chi }^2}{(1+2N\chi )(1+2\chi )}}& 0\\ 0& {\displaystyle \frac{2\chi }{1+2N\chi }}& 0& {\displaystyle \frac{4{\chi }^2}{(1+2N\chi )(1+2\chi )}}& \dots & 0& {\displaystyle \frac{4{\chi }^2}{(1+2N\chi )(1+2\chi )}}\\ {\displaystyle \frac{2\chi }{1+2N\chi }}& 0& {\displaystyle \frac{4{\chi }^2}{(1+2N\chi )(1+2\chi )}}& 0& \dots & {\displaystyle \frac{1+2N\chi +4{\chi }^2}{(1+2N\chi )(1+2\chi )}}& 0\\ 0& {\displaystyle \frac{2\chi }{1+2N\chi }}& 0& {\displaystyle \frac{4{\chi }^2}{(1+2N\chi )(1+2\chi )}}& \dots & 0& {\displaystyle \frac{1+2N\chi +4{\chi }^2}{(1+2N\chi )(1+2\chi )}}\end{array}\right].\end{array}$$

Using that expression, we arrive at (37). Taking partial derivatives of (37) leads to

$$\begin{array}{c}\nabla h={\displaystyle \frac{\partial h}{\partial \breve{\mu }}}=\left[\begin{array}{c}\partial h/\partial {\mu }_1\\ \partial h/\partial {\mu }_2\\ \dots \\ \partial h/\partial {\mu }_N\end{array}\right]=\left[\begin{array}{c}{(\nabla h)}_1\\ {(\nabla h)}_2\\ \dots \\ {(\nabla h)}_N\end{array}\right]\end{array}$$

with

$$\begin{array}{c}{(\nabla h)}_k={\displaystyle \frac{\partial h}{\partial {\mathit{\mu }}_k}}=\left[\begin{array}{c}\partial h/\partial {\mu }_{k1}\\ \partial h/\partial {\mu }_{k2}\\ \partial h/\partial {\mu }_{k3}\\ \partial h/\partial {\mu }_{k4}\\ \partial h/\partial {\mu }_{k5}\\ \partial h/\partial {\mu }_{k6}\end{array}\right]=\left[\begin{array}{c}\partial h/\partial {\mu }_{k1}\\ \partial h/\partial {\mu }_{k2}\\ 0\\ 1\\ 0\\ 0\end{array}\right],k=1,\dots ,N.\end{array}$$

We need to find the derivatives of h with respect to ${\mu }_{k1}$ and ${\mu }_{k2}$. We start with ${\mu }_{k1}$. We see that

$$\begin{array}{cc}\hfill \partial h/\partial {\mu }_{11}& ={\displaystyle \frac{1+2\chi }{1+2N\chi }}{\mu }_{11}+{\displaystyle \frac{2\chi }{1+2N\chi }}\left({\mu }_{21}+\dots +{\mu }_{N1}\right);\hfill \\ \hfill \partial h/\partial {\mu }_{k1}& ={\displaystyle \frac{2\chi }{1+2N\chi }}{\mu }_{11}+{\displaystyle \frac{1+2N\chi +4{\chi }^2}{(1+2N\chi )(1+2\chi )}}{\mu }_{k1}+{\displaystyle \frac{4{\chi }^2}{(1+2N\chi )(1+2\chi )}}{\displaystyle \sum _{i=2,i\ne k}^N}{\mu }_{i1}\hfill \\ \hfill & k=2,\dots ,N.\hfill \end{array}$$

Democracy

The matrix $\mathsf{\Psi }$ can be again computed explicitly as

$$\begin{array}{c}\mathsf{\Psi }={({\mathbb{I}}_N+2\chi B)}^{-1}\otimes {\mathbb{I}}_2\\ =\left[\begin{array}{ccccccccc}{\displaystyle \frac{1+2\chi }{1+2N\chi }}& 0& {\displaystyle \frac{2\chi }{1+2N\chi }}& 0& \dots & {\displaystyle \frac{2\chi }{1+2N\chi }}& 0& {\displaystyle \frac{2\chi }{1+2N\chi }}& 0\\ 0& {\displaystyle \frac{1+2\chi }{1+2N\chi }}& 0& {\displaystyle \frac{2\chi }{1+2N\chi }}& \dots & 0& {\displaystyle \frac{2\chi }{1+2N\chi }}& 0& {\displaystyle \frac{2\chi }{1+2N\chi }}\\ {\displaystyle \frac{2\chi }{1+2N\chi }}& 0& {\displaystyle \frac{1+2\chi }{1+2N\chi }}& 0& \dots & {\displaystyle \frac{2\chi }{1+2N\chi }}& 0& {\displaystyle \frac{2\chi }{1+2N\chi }}& 0\\ 0& {\displaystyle \frac{2\chi }{1+2N\chi }}& 0& {\displaystyle \frac{1+2\chi }{1+2N\chi }}& \dots & 0& {\displaystyle \frac{2\chi }{1+2N\chi }}& 0& {\displaystyle \frac{2\chi }{1+2N\chi }}\\ \dots & \dots & \dots & \dots & \dots & \dots & \dots & \dots & \dots \\ {\displaystyle \frac{2\chi }{1+2N\chi }}& 0& {\displaystyle \frac{2\chi }{1+2N\chi }}& 0& \dots & {\displaystyle \frac{1+2\chi }{1+2N\chi }}& 0& {\displaystyle \frac{2\chi }{1+2N\chi }}& 0\\ 0& {\displaystyle \frac{2\chi }{1+2N\chi }}& 0& {\displaystyle \frac{2\chi }{1+2N\chi }}& \dots & 0& {\displaystyle \frac{1+2\chi }{1+2N\chi }}& 0& {\displaystyle \frac{2\chi }{1+2N\chi }}\\ {\displaystyle \frac{2\chi }{1+2N\chi }}& 0& {\displaystyle \frac{2\chi }{1+2N\chi }}& 0& \dots & {\displaystyle \frac{2\chi }{1+2N\chi }}& 0& {\displaystyle \frac{1+2\chi }{1+2N\chi }}& 0\\ 0& {\displaystyle \frac{2\chi }{1+2N\chi }}& 0& {\displaystyle \frac{2\chi }{1+2N\chi }}& \dots & 0& {\displaystyle \frac{2\chi }{1+2N\chi }}& 0& {\displaystyle \frac{1+2\chi }{1+2N\chi }}\end{array}\right].\end{array}$$

Using that matrix $\mathsf{\Psi }$ and the expression (A15), we arrive to the expression for the Hamiltonian given by (38). We then compute the gradients of that Hamiltonian as follows:

$$\begin{array}{c}\nabla h={\displaystyle \frac{\partial h}{\partial \breve{\mu }}}=\left[\begin{array}{c}\partial h/\partial {\mu }_1\\ \partial h/\partial {\mu }_2\\ \dots \\ \partial h/\partial {\mu }_N\end{array}\right]=\left[\begin{array}{c}{(\nabla h)}_1\\ {(\nabla h)}_2\\ \dots \\ {(\nabla h)}_N\end{array}\right]\end{array}$$

$$\begin{array}{c}h={\displaystyle \sum _{k=1}^N}{\mu }_{k2}+{\displaystyle \frac{1}{2}}\left({\displaystyle \frac{1+2\chi }{1+2N\chi }}{\mu }_{11}^2+\dots +{\displaystyle \frac{1+2\chi }{1+2N\chi }}{\mu }_{N1}^2+2{\displaystyle \frac{2\chi }{1+2N\chi }}{\displaystyle \sum _{i,j=1,i<j}^N}{\mu }_{i1}{\mu }_{j1}\right).\end{array}$$

Hence

$$\begin{array}{c}\nabla h={\displaystyle \frac{\partial h}{\partial \breve{\mu }}}=\left[\begin{array}{c}\partial h/\partial {\mu }_1\\ \partial h/\partial {\mu }_2\\ \dots \\ \partial h/\partial {\mu }_N\end{array}\right]=\left[\begin{array}{c}{(\nabla h)}_1\\ {(\nabla h)}_2\\ \dots \\ {(\nabla h)}_N\end{array}\right],\end{array}$$

with the Hamiltonian given by the following explicit formulas:

$$\begin{array}{c}h={\displaystyle \sum _{k=1}^N}{\mu }_{k4}+{\displaystyle \frac{1}{2}}({\displaystyle \frac{1+2\chi }{1+2N\chi }}{\mu }_{11}^2+{\displaystyle \frac{1+2\chi }{1+2N\chi }}{\mu }_{12}^2+\dots +{\displaystyle \frac{1+2\chi }{1+2N\chi }}{\mu }_{N1}^2+{\displaystyle \frac{1+2\chi }{1+2N\chi }}{\mu }_{N2}^2\\ +2{\displaystyle \frac{2\chi }{1+2N\chi }}{\displaystyle \sum _{i,j=1,i<j}^N}{\mu }_{i1}{\mu }_{j1}+2{\displaystyle \frac{2\chi }{1+2N\chi }}{\displaystyle \sum _{i,j=1,i<j}^N}{\mu }_{i2}{\mu }_{j2}).\end{array}$$

Hence

$$\begin{array}{c}\nabla h={\displaystyle \frac{\partial h}{\partial \breve{\mu }}}=\left[\begin{array}{c}\partial h/\partial {\mu }_1\\ \partial h/\partial {\mu }_2\\ \dots \\ \partial h/\partial {\mu }_N\end{array}\right]=\left[\begin{array}{c}{(\nabla h)}_1\\ {(\nabla h)}_2\\ \dots \\ {(\nabla h)}_N\end{array}\right]\end{array}$$

We compute the gradients of h with respect to ${\mu }_{k1}$ for each particle as

$$\begin{array}{c}{(\nabla h)}_k={\displaystyle \frac{\partial h}{\partial {\mathit{\mu }}_k}}=\left[\begin{array}{c}\partial h/\partial {\mu }_{k1}\\ \partial h/\partial {\mu }_{k2}\\ \partial h/\partial {\mu }_{k3}\\ \partial h/\partial {\mu }_{k4}\\ \partial h/\partial {\mu }_{k5}\\ \partial h/\partial {\mu }_{k6}\end{array}\right]=\left[\begin{array}{c}\partial h/\partial {\mu }_{k1}\\ \partial h/\partial {\mu }_{k2}\\ 0\\ 1\\ 0\\ 0\end{array}\right],k=1,\dots ,N.\end{array}$$

We see that the derivative with respect to first component ${\mu }_{k1}$ for each particle is given by

$$\begin{array}{c}\partial h/\partial {\mu }_{k1}={\displaystyle \frac{1+2\chi }{1+2N\chi }}{\mu }_{k1}+{\displaystyle \frac{2\chi }{1+2N\chi }}{\displaystyle \sum _{i=1,i\ne k}^N}{\mu }_{i1}\forall k=1,\dots ,N.\end{array}$$

The derivative with respect to second component ${\mu }_{k2}$ is given by

$$\begin{array}{c}\partial h/\partial {\mu }_{k2}={\displaystyle \frac{1+2\chi }{1+2N\chi }}{\mu }_{k2}+{\displaystyle \frac{2\chi }{1+2N\chi }}{\displaystyle \sum _{i=1,i\ne k}^N}{\mu }_{i2}\forall k=1,\dots ,N.\end{array}$$

All in all, we see that the Hamiltonian h, and, hence, partial derivatives $\partial h/\partial {\mu }_{k1}$ and $\partial h/\partial {\mu }_{k2}$ are symmetric to each other with respect to the first and second components of all particles. This fact makes sense intuitively as well, since we suppose that control is applied equivalently in both components.