Evaluation I: Isolation Study on Flexibility Axes sec:eval

exp_T0.75Tile Axis Isolation. (a) Performance comparisons of accelerators with
different level of Tile flexibility running Mnasnet model. (b) Venn diagram of
Tile map space on Mnasnet. We use on-chip buffer size=4K in this experiments.
Scale: Total Data points in $W^{\omega}_{T}$ = $\pi(40)^{2}$.
Runtime/Energy/EDP are normalized against values of -1000. exp_T_sweep0.7The
runtime performance and workload Flexion of fully-Tile-flexible accelerators
with different buffer sizes. (W-F: Workload-dependent Flexion.) exp_O0.75Order
Axis Isolation. (a) Performance comparisons of accelerators with different
level of Tile flexibility running Mnasnet model. (b) Venn diagram of Tile map
space on Mnasnet. Runtime/Energy/EDP are normalized against values of -0100.
exp_P0.75Parallelism Axis Isolation. Performance comparisons of accelerators
with different level of Parallelism flexibility running Mnasnet model with (a)
16x64 PE array and (b) 32x32 PE array. (c) Venn diagram of Parallelism map
space on Mnasnet. PRuntime/Energy/EDP are normalized against values of -0010.
exp_S0.75Shape Axis Isolation. (a) Performance comparisons of accelerators
with different level of Tile flexibility running Mnasnet model. (b) Venn
diagram of Tile map space on Mnasnet. Runtime/Energy/EDP are normalized
against values of -0001. exp_S_sweep0.7The runtime performance and workload
Flexion of fully-Shape-flexible accelerators with different size of PE arrays.
(W-F: Workload-dependent Flexion) exp_flex0.755ptRuntime performance of
accelerators with different flexibility. -0000-X-Opt means class-0000 (fixed)
accelerator that is optimized for DNN model X (X= Alexnet, Mnasnet, etc.). The
runtime in each row is normalized by the runtime of the values of
-0000-Alexnet-Opt. In this section, we demonstrate how architects could
utilize the flexibility-aware toolflow to tractably explore different map-
space constraints on candidate designs. Simultaneously, we leverage the
toolflow constraints to target the DSE in order to isolate and quantify the
contributions of the T/O/P/S axes.

###  Methodology

Workloads. We use MnasNet [tan2019mnasnet], which is a high speed and high
accuracy model found by Neural Architecture Search (NAS), for our evaluation,
since networks found by NAS have been recently achieving state-of-the-art
performance.We also look into three types of models: Vision (Alexnet
[alexnet], Resnet50 [resnet], MobilenetV2 [sandler2018mobilenetv2]), Language
(BERT [bert]), and Recommendation (NCF [ncf], DLRM [dlrm]). Target
Accelerators. For the targeted isolation studies, we formulate three
accelerators with different levels of flexibility: i) an inflexible
accelerator (), using a fixed and unchangeable TOPS configuration, ii) a
partially flexible accelerator () with some constrained capabilities that we
describe in the relevant sections, and iii) a fully flexible accelerator (),
which can support the full map-space, as was described earlier in
sec:part_flex. These accelerators are built with the HW resources and baseline
configuration listed in table:baseline_dataflow, and the area overhead for
different flexibility is shown in table:area_cost. For each study, we
characterized the full neural network across all layers and present the
overall runtime, energy, and energy-delay product (EDP) for the best-found
mappings. Furthermore, we then present detailed analysis of three layers with
significantly different aspect ratios from MnasNet to illustrate the
correlation between performance. Optimization Objective. Across our
experiments, we set the optimization objective to minimum runtime. However,
other objectives such as energy, area, energy-delay-product (EDP),
performance-per-watt, and so on, are all feasible. We can also use multi-
objective formulation to optimize two or more objectives simultaneously. Map-
space Visualization. Finally, for each of the targeted isolation studies, we
present a Venn diagram in the style of fig:design_space_def drawn to scale
(provided in each figure) with the areas capturing the number of mappings.
Each Venn diagram plots the average of the metrics across all layers of the
model. Axis Isolation. We investigate the impact of standalone flexibility of
Tile (class-1000), Order (class-0100), Parallelism (class-0010) and Shape
(class-0001) with the target accelerators described by varying the constraints
given to the tool. This aspect of the evaluation is not intended to reflect
the typical usage of the tool, but rather as a limit study that characterizes
the most extreme points in the overall design space.

###  Isolation: Tile Flexibility

We study the impact of tile flexibility with three kinds of accelerators:
-1000 (or -0000111-X (X=0000 to 1111) indicate the same constraint in all 16
classes.), uses the fixed baseline tile sizes; -1000 uses hard-partitioned
buffer with partition ratio of 1:1:1 for weight, input, and output buffer, and
uses flexible tile sizes; -1000 uses soft-partitioned buffer and flexible tile
sizes. Note that the hardware trade-off between hard-partitioning and soft-
partitioning significantly increases the energy per on-chip memory access, but
can also result in fewer expensive off-chip accesses [parashar2019timeloop].
Map Space. Hardware-dependent Flexion (H-F) $A_{X}/C_{X}$ is decided by the
partitioning strategy which is around 0.22 for 1:1:1 partition. Value 0.22 is
the area ratio of $A_{X}$ over $C_{X}$ in fig:exp_T(b). $A^{\omega}_{X}$ is
the intersection of $W^{\omega}_{X}$ and $A_{X}$, showing how many different
tiles that workload $\omega$ can explore and are supported by the
accelerators. $A^{\omega}_{X}$ of -1000 is 1, since it leverages only 1
choice, while $A^{\omega}_{X}$ of -1000 will be smaller than the one for
-1000, since hard-partitioning is a stricter constraint. These can be observed
by the increase of overlapped area of $W^{\omega}_{X}$ and $A_{X}$ from -1000
to -1000. Note that the workload-dependent Flexion (W-F)
$A^{\omega}_{X}/W^{\omega}_{X}$ is directly related to the ratio of the
supported map space to the entire map space, where larger values implies the
given flexibility in the investigating accelerator can support the current
workload better. Performance. -1000 allows the tile sizes to be flexible,
which reduce runtime by 98%, as shown in Layer-1 (L1) in fig:exp_T(a). It
chooses smaller K, C sizes and larger Y, X sizes to tailor for Layer-1. While
enabling soft-partitioning (-1000), with more tile choices, the runtime
improves by another $6\times$. In Layer-16, the C tile sizes are chosen to
divide the C dimension size perfectly. The tile configurations are chosen
under the criterion of best runtime, and thus energy does not show significant
difference. However, if we used least energy as criterion, different points
would be picked, demonstrating a larger energy difference between -1000 and
-1000/-1000. Overall, with this criterion, in class-1000, -1000 can achieve
speedup of $4.8\times$ and $1.9\times$ over -1000 and -1000 in end-to-end
Mnasnet. Sensitivity Analysis: Buffer Size. The space of supported tile shapes
is directly affected by the buffer size, which is the intersection area of
$S_{X}$ and $W^{\omega}_{S}$ in fig:exp_T(b). As shown in fig:exp_T_sweep,
with increasing buffer sizes, the workload-dependent Flexion increases and the
runtime improves. We can observe that the runtime improvement saturates at a
buffer size of around 6.4KB, which is sufficient for most of the layers in
MnasNet. Conversely, when the buffer size is small, the flexibility to execute
on soft-partitioned buffer (-1000) becomes more valuable to achieve a good
runtime. Takeaways. For a model like Mnasnet with large diversity in layer
dimensions, we find that tile flexibility is crucial for the accelerator to
capture the map-space. For the most part, the hard-partitioned -1000 provides
significant runtime benefits over an -1000 but is strictly worse than -1000.
We also observe that tile flexibility is more crucial for accelerators with
smaller buffer sizes.

###  Isolation: Order Flexibility

The three types of accelerators are configured as follows: -0100, using the
output-stationary order YXKCRS, -0100, with three order choices including
output, input, and weight stationary order, and -0100, with flexible order.
Map Space. H-F $A_{X}/C_{X}$ is affected by the number of supported order
choices, which depends on the available HW support described in
fig:hw_cost_pic. For order, the map space does not depend on the number of
compute or memory resources such as PEs and buffer sizes. Therefore, $C_{X}$
can encompasses all design points in $W^{\omega}_{X}$, as shown in
fig:exp_O(b). Performance. From Layer-16 in fig:exp_O(a), we observe that
simply adding 3 out of 6! order choices—and consequently increasing workload
Flexion by a small amount—improves runtime to a near-optimal value. This
demonstrates that for accelerator design, partially supporting order
flexibility may expose a better cost-performance trade-off, if the supported
order choices expose distinct behaviors. In our evaluation, we also found that
many different orders will end up with the similar runtime. This observation
motivates the design of –0100, a simpler and lower-cost accelerator that
achieves similar performance as –0100. Overall, in class-0100, full
flexibility support can still achieve $1.12\times$ and $1.01\times$ speedup
over the two baselines. Takeaways. We find that the output stationary loop-
order in -0100 works reasonably well across all layers. Adding support for
three loop-orders (output, weight and input) in -0100 gets nearly the same
performance as supporting all 6! loop orders in -0100.

###  Isolation: Parallelism Flexibility

The three types of accelerators are as follows: -0010, using default K-C
parallelism, -0010, with a choice of K-C or Y-X parallelism, and -0010. Map
Space. Since we investigate CONV-accel (6-dim) and consider two-way
parallelism. We have $C_{X}$=6x5=30 different choices, while the $A_{X}$ for
-0010 and -0010 is 1 and 2, respectively. Performance. fig:exp_P(a) shows that
different workloads can have different optimum parallelism choices, and
supporting only a subset of choices may not be sufficient to approach optimal
performance (Layer-16, 19). Layer-29 is a depth-wise CONV, where there is no
cross-input channel computation. -0010’s restriction of parallelism only
across K and C dimensions deprives it of other parallelism opportunities such
as YX- or RS-parallelism, which could be a better choice for this layer.
Overall in class-0010, -0010 consistently achieves around $1.6\times$ and
$1.3\times$ better runtime relative to the fixed and partially flexible
accelerators, respectively. Takeaways. We observe the flexibility in
parallelism is highly-sensitive to the shape of the physical array. We also
see that non-conventional parallelism choices (such as XK or KS) are found by
the mapper, suggesting that supporting full flexibility in parallelism choices
is valuable.

###  Isolation: Shape Flexibility

The accelerators compared are as follows: -0001, using a square PE array,
-0001-A, a flexible PE array composed of a modular $16\times 16$ PE-array
building block, -0001-B, a flexible PE array with a $4\times 4$ building
block, and -0001. Map Space. $C_{X}$ is decided by the size of the PE array.
$A_{X}$ for the -0001 is the number of different array shapes that can be
composed by the smaller building block. A smaller building block can enable
more fine-grained PE array shape exploration and thus exposes a larger
configuration space, which can be observed by the different size of $A_{X}$
for the two partially flexible accelerators in fig:exp_S(b). Performance. When
the size of parallelism dimension is larger than the size of the building
block ($32\times 32$, $16\times 16$, $4\times 4$, or flexible), we need to
fold the computation. When they are not perfectly divisible, the last fold of
computation will introduce compute under-utilization, as shown in Layer-15 and
Layer-16 of -0001-A in fig:exp_S(a). Interestingly, -0001-B can mostly reach
optimum performance, since most of the layers have K, C dimension sizes that
are multiples of 4. Note that we use the default K-C parallelism strategy in
this experiment. -0001-B reaches almost the same performance as the fully-
flexible one with only $6\%$ flexibility. Overall, we observe that full
flexibility achieve speedup of $1.05\times$, $1.02\times$, and $1.001\times$
over the three baselines, respectively. Sensitivity Analysis: Array Size. A
larger number of PEs can allow more array configuration options and hence has
larger hardware Flexion, as shown in fig:exp_S_sweep. We can also observe the
diminishing return around $45\times 45$ to $64\times 64$ points when the
provided number of PEs is starting to become larger than the maximum K-C
parallelism that some of the layers can leverage, leading to under utilization
of the PEs. Takeaways. We observe the flexibility in array shape is highly
crucial for utilization across layers. However, it is highly sensitive to the
dimensions being mapped spatially. This indicates that flexibility in P and S
go hand-in-hand (more discussion on this in sec:future_proof). One key
observation is that a partially flexible accelerator using multiple fixed-size
small arrays (e.g., 4x4) reaches almost the same performance as the fully-
flexible one. This observation is consistent with some design choices in the
existing accelerator designs, e.g., the $4\times 4$ matrix multiplier in CUDA
tensor cores [cuda_guide]. The partially flexible accelerators are also
sufficient for many manual-design CNNs, e.g., ResNet50 [resnet], VGG16 [vgg],
whose dimension sizes are mostly multiples of 64 in K, C dimensions. However,
shape flexibility could become increasingly important because of two use
cases, as follows. (1) More and more neural network designs are relying on
NAS, which could generate more diverse layer shapes. MnasNet is one such NAS
network, but it is highly modularized, while more diverse networks could be
developed in future to target different applications. (2) Pruning techniques
can change layer shapes in diverse ways.

##  1 Evaluation II: Accelerator Future-Proofing

In this experiment, we demonstrate how this work enables a first-of-its-kind
experiment: assessing the advantage of adding flexibility to an accelerator at
design time in terms of future-proofing, quantifying the trade-off between
hardware cost and the performance gain on “future” workloads. To start, we
design an inflexible accelerator tailored for Alexnet, -0000-Alexnet-Opt, by
finding a fixed configuration of TOPS to optimize its runtime. This represents
an accelerator such as Eyeriss [chen2016eyeriss_jssc] which was a optimized
for a fixed workload in a certain year ($\sim$2014). Fixed Alexnet-Opt Accel
on Future Models. We now progress our accelerator into the future. When new
neural network models (model X) are introduced. We design -0000-X-Opts,
representing a new accelerator optimized just for this new workload—though for
comparison purposes we always use roughly the same number of hardware
resources as -0000-Alexnet-Opt. First, we use the -0000-Alexnet-Opt to run the
future models and compare the performance witht the future model (X) optimized
-0000-X-Opt. We observe that -0000-Alexnet-Opt can achieve similar performance
to -0000-X-Opt when X is a CNN model, which follows the intuition since
Alexnet is also a CNN, which might share similar characteristics. However, for
fully-connected (FC) layer dominated models such as BERT, DLRM, and NCF, we
see 1/0.25 = 4x to 1/0.07=14x slowdown. This shows that a fixed configuration
can hurt performance when new models start to show notably different
characteristics (e.g., CONV- vs FC-dominated). Single-Axes Flexible Alexnet-
Opt Accel on Future Models. Next, we go back to the “past” to explore the
effect of adding single-axes flexibility when designing Alexnet-Opt
accelerator, as shown in rows 3-6 of fig:exp_flex. We see that in MnasNet,
-1000-Alexnet-Opt is $1.00/0.57=1.75\times$ better than -0000-Alexnet-Opt. We
found tile flexibility to be crucial for the three CNNs and BERT. However, it
is not as effective in DLRM and NCF. The reason is as follows. In this
experiment, -0000-Alexnet-Opt found Y-K parallelism to be optimal for Alexnet,
and this stays valid for other CNNs, too. BERT is composed of matrix-matrix
multiplications, the GEMM (M,N,K) dimensions got mapped to the (K_conv,C,Y)
dimensions of the CONV accelerator respectively. The Y-K parallelism in the
accelerator helped BERT as well. However, DLRM and NCF use matrix-vector
multiplications, where the K dimension is 1 for the vector. Mapping them on
the -0000-Alexnet-Opt accelerator leads to under-utilization along the
Y-parallelism dimension of the accelerator, leading to loss in performance. In
order to run DLRM and NCF efficiently, the ability to switch to other
parallelism strategies become crucial. Multiple-Axes Flexible Alexnet-Opt
Accel on Future Models. There are 10 other variants of combinatorial flexible
accelerators, as shown in row 7-16 of fig:exp_flex. We see that (P)+(S)
together (0011) work well and show significant performance gain compared to
the separated counterpart (0010, 0001). (T)+(P) with only two flexibility axes
enabled can achieve comparable performance with -1111-Alexnet-Opt in multiple
cases. Interestingly, we see that (T)+(O) works much worse than (T)+(P). This
shows that (T), which individually works well still needs to find the right
other flexibility axes to have another level of performance boost. Fully
Flexible Alexnet-Opt Accel on Future Models. Finally, we investigate two kinds
of flexible accelerators: -1111-Alexnet-Opt, and -1111-Alexnet-Opt.
-1111-Alexnet-Opt uses the configurations described in the evaluations in
sec:eval (using variant A for Shape). While more flexibility support
introduces larger area cost, it is noteworthy that it brings with potential to
leverage better runtime _even in the original AlexNet_. This is because
AlexNet itself has a lot of variation across its layers. _A fixed accelerator
that optimizes for the average-case across all layers may miss the potential
to reach the best case for each individual layer._ For e.g., -1111-Alexnet-Opt
speeds up the baseline -0000-Alexnet-Opt by 1/0.06=16.7x times in Alexnet.
Furthermore, -1111-Alexnet-Opt constantly achieves better performance than
-0000-X-Opts. Takeaways. Accelerators optimized for certain layers (such as
CONV2D) are more future proof, even if the layer dimensions of future models
change. This is because most models show abundant parallelism across the input
and output channels. Often tile flexibility seems sufficient to capture most
benefits. However, for FC-dominated layer types, flexibility in parallelism
and shape starts reaping more benefits since there are fewer dimensions and
GEMMs can be irregular. As models with several other layer types emerge (e.g.,
FC/DWCONV/Attention/…), investing in flexibility is valuable.


# Metriplectic Heavy Top: An Example of Geometrical Dissipation

By Michael Updike
Supervised by Dr. P. J. Morrison  Undergraduate Thesis
Department of Physics, The University of Texas at Austin

(December 2022)

## Abstract

Recently, Morrison and Updike [10] showed that many dissipative systems are
naturally described as possessing a Riemann curvature-like bracket, which
similar to the Poisson bracket, generates the dissipative equations of motion
once suitable generators are chosen. In this paper, we use geometry to
construct and explore the dynamics of these new brackets. Specifically, we
consider the dynamics of a heavy top with dissipation imposed by a Euclidian
contravariant curvature. We find that the equations of motion, despite their
rather formal motivation, naturally generalize the energy-conserving
dissipation considered by Matterasi and Morrison [4]. In particular, with
suitable initial conditions, we find that the geometrically motivated
equations of motion cause the top to relax to rotation about a principal axis.

## Definitions and Conventions

All repeated indices are assumed to be summed over unless otherwise specified.
We denote the phase space manifold of a system $\mathcal{Z}$. We use $z^{i}$
to denote the coordinates of $\mathcal{Z}$. The symbols $f,g,h,F,G,K,$ and $N$
are always used to represent smooth ($C^{\infty}(\mathcal{Z})$) functions. The
tangent bundle of $\mathcal{Z}$ is denoted $T\mathcal{Z}$, and the cotangent
bundle $T^{*}\mathcal{Z}$.

Here, $\alpha,\beta,$ and $\gamma$ always represent one-forms. The space of
one-forms is denoted $\Omega^{1}(\mathcal{Z})$. The symbol $\mathbf{d}$
represents the exterior derivative, which acts on a function $f(z)$ as

$\mathbf{d}f=\frac{\partial f}{\partial z^{i}}\mathbf{d}{z^{i}}.$

The exterior derivative of a coordinate function is strictly formal. In
coordinates, a general one-form may be written

$\alpha=f_{i}(z)\;\mathbf{d}z^{i}.$

A one-form $\alpha$ is said to be exact if $\alpha=\mathbf{d}f$ for some $f$.

Given a tangent vector $v$, we use $v[f]$ to mean $v$ acting on $f$. We use
$\left<\beta,v\right>$ to denote the canonical pairing between a vector and a
one-form.

We use the word bracket to mean a smooth map from some number of smooth
functions to a single smooth function. We demand that a bracket is both linear
and a derivation in each of its arguments. In particular, we use
$\\{\cdot,\cdot\\}\colon$ $C^{\infty}(\mathcal{Z})\times
C^{\infty}(\mathcal{Z})\to C^{\infty}(\mathcal{Z})$ to denote the Poisson
bracket. The Poisson bracket gives a bivector field

$J^{ij}=J(\mathbf{d}z^{i},\mathbf{d}z^{j})\colon=\\{z^{i},z^{j}\\}$

and also an anchor map $J^{\\#}:T^{*}\mathcal{Z}\to T\mathcal{Z}$ defined
implicitly by

$\left<\beta,J^{\\#}(\alpha)\right>\colon=J(\alpha,\beta).$

We use $D$ to represent a contravariant connection, also called a
contravariant derivative, which gives the derivative of a-one form with
respect to another one-form. A contravariant derivative is similar to but not
the same as a covariant derivative.

## Introduction

Given a degenerate Poisson bracket $\\{\cdot,\cdot\\}$, there exist
distinguished functions $S$ called Casimirs that Poisson commute with all
functions $f$

$\displaystyle\\{f,S\\}=0.$ (1)

Casimirs represent, potentially, the entropy of a system. Given that the
dynamics of any Hamiltonian system with the Poisson bracket
$\\{\cdot,\cdot\\}$ is constrained to surfaces of constant $S$, a dissipative
structure is required to create dynamics that respect the second law of
thermodynamics. First developed by Morrison and others, metriplectic dynamics
is a systematic way to add energy-conserving dissipation to an otherwise
Hamiltonian system (cf.[7]).

Suppose we have a Hamiltonian system with Hamiltonian $H$ and a Casimir $S$.
We consider a symmetric bracket $(\cdot,\cdot)$ such that, for all functions
$f$,

$\displaystyle(f,H)=0.$ (2)

The metriplectic equations of motion that both preserve energy and increase
$S$ are

$\displaystyle\dot{f}=\\{f,H\\}+(f,S).$ (3)

In this paper, we first use the framework of both Riemannian and Poisson
geometry to rephrase metriplectic dynamics as a geometrical theory. We then
use this new approach to metriplectic dynamics to recreate a bracket first
introduced by [6] and used in [4] for control of a rigid body. Afterward, we
naturally generalize the construction to the heavy top system, constructing a
family of dissipative theories which we collectively call the metriplectic
heavy top. Finally, we simulate the dynamics of a handful of these theories,
finding asymptotic relaxation of the heavy top.

## Riemann-Poisson Geometry

Suppose we have system described on some Poisson manifold
$(\mathcal{Z},\\{\cdot,\cdot\\})$ with coordinates $z^{i}$. The Poisson
bracket on functions is naturally extended to $T^{*}\mathcal{Z}$ by the Koszul
bracket $[\cdot,\cdot]$, which for exact forms reads [2][1]

$\displaystyle[\mathbf{d}z^{i},\mathbf{d}z^{j}]\colon=\mathbf{d}\\{z^{i},z^{j}\\}.$
(4)

The Kozul bracket gives us a natural Lie bracket on forms. Given a
(pseudo-)metric

$\displaystyle g^{ij}=g(\mathbf{d}z^{i},\mathbf{d}z^{j})$ (5)

we can define a contravariant Levi-Civita connection $D$ on $\mathcal{Z}$ via
the formula

$\displaystyle 2g(D_{\alpha}\beta,\gamma)=$ $\displaystyle
J^{\\#}(\alpha)[g(\beta,\gamma)]-J^{\\#}(\gamma)[g(\alpha,\beta)]$ (6)
$\displaystyle+J^{\\#}(\beta)[g(\gamma,\alpha)]+g([\alpha,\beta],\gamma)-g([\beta,\gamma],\alpha)+g([\gamma,\alpha],\beta)$

where $\alpha,\beta,\gamma\in\Omega^{1}(M)$ and
$J^{\\#}:T^{*}\mathcal{Z}\rightarrow T\mathcal{Z}$ is the anchor map defined
by

$\displaystyle\left<\beta,J^{\\#}(\alpha)\right>=J(\alpha,\beta).$ (7)

In coordinate form

$\displaystyle J^{\\#}(\mathbf{d}f)[g]=\frac{\partial f}{\partial
z^{i}}J^{ij}\frac{\partial g}{\partial z^{j}}.$ (8)

The Kozul bracket on general one-forms can be obtained using the formula

$\displaystyle[\alpha,f(z)\beta]=f(z)[\alpha,\beta]+J^{\\#}(\alpha)[f(z)]\beta.$
(9)

It should be noted $D$ is not a covariant derivative, which is defined with
tangent vectors. Even so, it can be shown $D$ satisfies similar linearity
properties

$\displaystyle D_{\alpha+\beta}\gamma$
$\displaystyle=D_{\alpha}\gamma+D_{\beta}\gamma$ (10) $\displaystyle
D_{f\alpha}\gamma$ $\displaystyle=fD_{\alpha}\gamma$ $\displaystyle
D_{\alpha}(\beta+\gamma)$ $\displaystyle=D_{\alpha}\beta+D_{\alpha}\gamma$
$\displaystyle D_{\alpha}(f\gamma)$
$\displaystyle=fD_{\alpha}\gamma+J^{\\#}(\alpha)[f]\gamma.$

Letting $\alpha=\mathbf{d}f$, $\beta=\mathbf{d}h$, and
$\gamma=\mathbf{d}z^{\delta}$, we can write the contravariant derivate as

$\displaystyle 2g^{\delta a}(D_{\mathbf{d}f}\mathbf{d}h)_{a}=\,$
$\displaystyle\frac{\partial f}{\partial z^{i}}J^{ij}\frac{\partial}{\partial
z^{j}}\left[g^{k\delta}\frac{\partial h}{\partial
z^{k}}\right]+J^{i\delta}\frac{\partial}{\partial
z^{i}}\left[g^{kl}\frac{\partial f}{\partial z^{k}}\frac{\partial h}{\partial
z^{l}}\right]$ (11) $\displaystyle+$ $\displaystyle\frac{\partial h}{\partial
z^{i}}J^{ij}\frac{\partial}{\partial z^{j}}\left[g^{k\delta}\frac{\partial
f}{\partial z^{k}}\right]+g^{k\delta}\frac{\partial}{\partial
z^{k}}\left[J^{ij}\frac{\partial f}{\partial z^{i}}\frac{\partial g}{\partial
z^{j}}\right]$ $\displaystyle-g^{kl}\frac{\partial}{\partial
z^{k}}\left[J^{i\delta}\frac{\partial h}{\partial z^{i}}\right]\frac{\partial
f}{\partial z^{l}}-g^{kl}\frac{\partial}{\partial
z^{k}}\left[J^{jk}\frac{\partial f}{\partial z^{j}}\right]\frac{\partial
h}{\partial z^{l}}.$

It is convenient to introduce

$\displaystyle\Gamma^{ij}_{l}\colon=(D_{\mathbf{d}z^{i}}\mathbf{d}z^{j})_{k}=\,$
$\displaystyle\frac{1}{2}g_{kl}\left[J^{ia}\frac{\partial g^{jk}}{\partial
z^{a}}-J^{ka}\frac{\partial g^{ij}}{\partial z^{a}}+J^{ja}\frac{\partial
g^{ik}}{\partial z^{a}}\right]$ (12)
$\displaystyle+\frac{1}{2}g_{kl}\left[g^{ka}\frac{\partial J^{ij}}{\partial
z^{a}}-g^{ai}\frac{\partial J^{jk}}{\partial z^{a}}-g^{aj}\frac{\partial
J^{ik}}{\partial z^{a}}\right].$

$D$ is the unique connection that is both torison-free

$\displaystyle D_{\alpha}\beta-D_{\beta}\alpha=[\alpha,\beta]$ (13)

and metric compatible

$\displaystyle
J^{\\#}(\alpha)[g(\beta,\gamma)]=g(D_{\alpha}\beta,\gamma)+g(\beta,D_{\alpha}\gamma).$
(14)

In coordinates, these conditions read

$\displaystyle\Gamma^{ij}_{k}-\Gamma^{ji}_{k}=\frac{\partial J^{ij}}{\partial
z^{k}}$ (15)

and

$\displaystyle J^{ia}\frac{\partial g^{jk}}{\partial
z^{j}}=g^{ka}\Gamma^{ij}_{a}+g^{ja}\Gamma_{a}^{ik}.$ (16)

We define the contravariant Riemann curvature tensor in a way formally
reminiscent of the usual curvature tensor

$\displaystyle R(\alpha,\beta)\gamma=D_{\alpha}D_{\beta}\gamma-
D_{\beta}D_{\alpha}\gamma-D_{[\alpha,\beta]}\gamma.$ (17)

This tensor gives us a natural 4-bracket
$(\cdot,\cdot;\cdot,\cdot):C^{\infty}(\mathcal{Z})\times
C^{\infty}(\mathcal{Z})\times C^{\infty}(\mathcal{Z})\times
C^{\infty}(\mathcal{Z})\rightarrow C^{\infty}(\mathcal{Z})$ defined by

$\displaystyle(F,K;G,N)\colon=g(R(\mathbf{d}F,\mathbf{d}K)\mathbf{d}G,\mathbf{d}N)$
(18)

or in coordinates, noting the raised index,

$\displaystyle(F,K;G,N)=R^{ijkl}\frac{\partial F}{\partial
z^{i}}\frac{\partial K}{\partial z^{j}}\frac{\partial G}{\partial
z^{k}}\frac{\partial N}{\partial z^{l}}.$ (19)

This bracket inherits the following symmetries

$\displaystyle(F,K;G,N)=-(K,F;G,N)=-(F,K;N,G)=(G,N;F,K)$ (20)

in addition to the symmetries obtained by the first and second Bianchi
identity (see e.g.[11]).

Given a Hamiltonian $H$, we can define a bracket analogous to the operator
describing geodesic deviation

$\displaystyle(F,G)\colon=(F,H;G,H).$ (21)

Notice, that by the symmetries of the 4-bracket, $(F,G)$ is both symmetric and
has $H$ in its kernel, precisely the necessary conditions for a metriplectic
bracket.

## Free Rigid Body Bracket

Before we use the formalism of Riemann-Poisson Geometry to add dissipation to
the heavy top, we first attempt to understand the geometry of the simpler free
rigid body. Using angular momentum as phase space coordinates ($z^{i}=L^{i}$,
$\mathcal{Z}=\mathbb{R}^{3}$), the Poisson bracket for the FRB system realizes
the $\mathfrak{so}(3)$ algebra

$\displaystyle\\{L^{i},L^{j}\\}=-\epsilon^{ijk}L^{k}.$ (22)

The Hamiltonian for this system is

$\displaystyle
H=\frac{1}{2}\left(\frac{L^{1}}{I_{1}}+\frac{L^{2}}{I_{2}}+\frac{L^{3}}{I_{3}}\right).$
(23)

To introduce dissipation to this system, we first consider a dissipative
metric for the system. The simplest possible choice is the Euclidian (or
equivalently the Cartan-killing) metric

$\displaystyle g^{ij}=\delta^{ij}$ (24)

For a general Lie-Poisson bracket corresponding to a semi-simple, compact
algebra with structure constants $C^{ij}_{k}$, the contravariant derivative
takes a very simple form when we use the Cartan-Killing form as the metric

$\displaystyle\Gamma^{ij}_{k}=\frac{1}{2}C^{ij}_{k}.$ (25)

The totally contravariant curvature tensor is also quite simple

$\displaystyle
R^{ijk}_{l}=\frac{1}{4}C^{jk}_{a}C^{ia}_{l}-\frac{1}{4}C^{ik}_{a}C^{jal}+\frac{1}{2}C^{ij}_{a}C^{ka}_{l}.$
(26)

For the $\mathfrak{so}(3)$ case at hand, the structure constants are given by
$C^{ij}_{k}=-\epsilon_{ijk}$ and our curvature tensor is given by

$\displaystyle
R^{ijkl}=\frac{1}{4}(\delta^{ik}\delta^{jl}-\delta^{il}\delta^{jk}).$ (27)

The metriplectic bracket is, up to a constant,

$\displaystyle(F,G)=(\omega^{2}\delta_{ij}-\omega_{i}\omega_{j})\frac{\partial
F}{\partial L^{i}}\frac{\partial G}{\partial L^{j}}.$ (28)

Amazingly, this is the bracket first constructed in [6] and later considered
in Matterasi and Morrison [4] as a way to model an energy-conserving torque
driving the body to rotate about a principle axis.

## Heavy Top Bracket

Inspired by the free-rigid body, we consider the heavy top system describing a
rigid body in a constant gravitational field, which is given by the body frame
vector $\mathbf{\Gamma}$. The Poisson bracket is [13][12]

$\displaystyle\\{F,G\\}=-\epsilon_{ijk}L^{i}\frac{\partial F}{\partial
L^{j}}\frac{\partial G}{\partial
L^{k}}-\epsilon_{ijk}\Gamma^{i}\left(\frac{\partial F}{\partial
L^{j}}\frac{\partial G}{\partial\Gamma^{k}}-\frac{\partial G}{\partial
L^{j}}\frac{\partial F}{\partial\Gamma^{k}}\right).$ (29)

In analogy with the free rigid body system from earlier, it is only natural
(and as we will see, desirable) to assume the heavy top also has a Euclidean
metric

$\displaystyle g^{ij}=\delta^{ij}.$ (30)

Organizing the phase space coordinates
$\mathbf{z}^{i}=(\mathbf{L},\mathbf{\Gamma})^{i}$, the Lie-Poisson structure
of the heavy top bracket allows us to write the bracket as

$\displaystyle\\{z^{i},z^{j}\\}=C^{ij}_{k}z^{k}$ (31)

where $C^{ij}_{k}$ are the structure constants of the semi-direct product
algebra $\mathfrak{so}(3)\ltimes\mathbb{R}^{3}$.

For any Lie-Poisson bracket with the Euclidean metric, the connection
coefficients are (a similar formula can be found in [5])

$\displaystyle\Gamma^{ij}_{k}=\frac{1}{2}(C^{ij}_{k}-C^{jk}_{i}+C^{ki}_{j})$
(32)

and the curvature tensor is

$\displaystyle R^{ijk}_{l}=$
$\displaystyle\frac{1}{4}(C^{jk}_{a}-C^{ka}_{j}+C^{aj}_{k})(C^{ia}_{l}-C^{al}_{i}+C^{li}_{a})$
(33)
$\displaystyle-\frac{1}{4}(C^{ik}_{a}-C^{ka}_{i}+C^{ai}_{k})(C^{ja}_{l}-C^{al}_{j}+C^{lj}_{a})-\frac{1}{2}C^{ij}_{a}(C^{ak}_{l}-C^{kl}_{a}+C^{la}_{k}).$

In general, a six-dimensional system like the heavy top will have $105$
independent components in the curvature tensor, and over a $1000$ dependant
terms. Fortunately, calculations like these are easily done computationally.
For the heavy top, we find that the curvature tensor is quite sparse with the
nonzero terms being

$\displaystyle\frac{1}{4}$
$\displaystyle=R^{1212}=R^{1313}=R^{2121}=R^{2323}=R^{3131}=R^{3232}$ (34)
$\displaystyle=-R^{2112}=-R^{3113}=-R^{1221}=-R^{3223}=-R^{1331}=-R^{2332}.$

This expression is even simpler than it appears. Every term involving a
$\mathbf{\Gamma}$ index vanishes. Furthermore, for $i,j,k,l\in\\{1,2,3\\}$

$\displaystyle
R^{ijkl}=\frac{1}{4}(\delta^{ik}\delta^{jl}-\delta^{il}\delta^{jk}).$ (35)

That is, restricted to $\mathbf{L}$ indices only, the curvature tensor is
exactly that of the free rigid body.

From here on, we assume the top is symmetric ($I_{1}=I_{2}$).The Hamiltonian
for the symmetric heavy top system is given by

$\displaystyle
H=\frac{1}{2}\left(\frac{L^{1}}{I_{1}}+\frac{L^{2}}{I_{1}}+\frac{L^{3}}{I_{3}}\right)+\xi\Gamma^{3}.$
(36)

Plugging this into the curvature tensor, the heavy top metriplectic bracket is
precisely the same bracket considered in Morrison and Matterasi to describe
the FRB

$\displaystyle(F,G)=(\omega^{2}\delta^{ij}-\omega^{i}\omega^{j})\frac{\partial
F}{\partial L^{i}}\frac{\partial G}{\partial L^{j}}.$ (37)

Unlike the free rigid body, the total angular momentum $L^{2}$ is no longer a
Casimir invariant. Rather, there are two new Casimirs up to composition with
an analytic function [13]

$\displaystyle\vec{\Gamma}\cdot\vec{L}$ (38)

and

$\displaystyle\vec{\Gamma}^{2}.$ (39)

The latter choice of Casimir trivializes the metriplectic dynamics, so the
only suitable choice of generating function for the dissipation has the form

$\displaystyle S=C(\vec{\Gamma}\cdot\vec{\omega})$ (40)

where $C$ is analytic.

The dissipative equations of motion are

$\displaystyle\dot{f}=\\{f,H\\}+(f,S)=$ (41)
$\displaystyle-\epsilon^{ijk}L^{i}\frac{\partial f}{\partial
L^{j}}\omega^{k}-\Gamma^{i}\left(\epsilon^{ij3}\frac{\partial f}{\partial
L^{j}}\xi-\epsilon^{ijk}\omega^{j}\frac{\partial
f}{\partial\Gamma^{k}}\right)+C^{\prime}(\vec{\Gamma}\cdot\vec{\omega})(\omega^{2}\delta^{ij}-\omega^{i}\omega^{j})\frac{\partial
f}{\partial L^{j}}\Gamma^{i}.$

We can expand this equation as a system of $6$ ODEs

$\displaystyle\dot{L}^{1}$
$\displaystyle=-L^{2}L^{3}\left(\frac{1}{I_{1}}-\frac{1}{I_{3}}\right)+\xi\Gamma^{2}+C^{\prime}(\mathbf{\Gamma}\cdot\mathbf{L})\left(\left[\frac{(L^{2})^{2}}{I_{1}^{2}}+\frac{(L^{3})^{2}}{I_{3}^{2}}\right]\Gamma^{1}-\frac{L^{1}L^{2}}{I^{2}_{1}}\Gamma^{2}-\frac{L^{1}L^{3}}{I_{1}I_{3}}\Gamma^{3}\right)$
(42) $\displaystyle\dot{L}^{2}$
$\displaystyle=L^{1}L^{3}\left(\frac{1}{I_{1}}-\frac{1}{I_{3}}\right)-\xi\Gamma^{1}+C^{\prime}(\mathbf{\Gamma}\cdot\mathbf{L})\left(-\frac{L^{1}L^{2}}{I^{2}_{1}}\Gamma^{1}+\left[\frac{(L^{1})^{2}}{I_{1}^{2}}+\frac{(L^{3})^{2}}{I_{3}^{2}}\right]\Gamma^{2}-\frac{L^{2}L^{3}}{I_{1}I_{3}}\Gamma^{3}\right)$
$\displaystyle\dot{L}^{3}$
$\displaystyle=C^{\prime}(\mathbf{\Gamma}\cdot\mathbf{L})\left(-\frac{L^{1}L^{3}}{I_{1}I_{3}}\Gamma^{1}-\frac{L^{2}L^{3}}{I_{1}I_{3}}\Gamma^{2}+\left[\frac{(L^{2})^{2}}{I_{1}^{2}}+\frac{(L^{1})^{2}}{I_{1}^{2}}\right]\Gamma^{3}\right)$
$\displaystyle\dot{\Gamma}^{1}$
$\displaystyle=\Gamma^{2}\frac{L^{3}}{I_{3}}-\Gamma^{3}\frac{L^{2}}{I_{1}}$
$\displaystyle\dot{\Gamma}^{2}$
$\displaystyle=\Gamma^{3}\frac{L^{1}}{I_{1}}-\Gamma^{1}\frac{L^{3}}{I_{3}}$
$\displaystyle\dot{\Gamma}^{3}$
$\displaystyle=\Gamma^{1}\frac{L^{2}}{I_{1}}-\Gamma^{2}\frac{L^{1}}{I_{1}}.$

Notice that the dissipation does not affect the equation of motion for
$\mathbf{\Gamma}$, as we should expect from any physically realistic system.
In fact, had we naively used a bracket projecting out the Hamiltonian, this
would no longer be true. It is worthwhile to note that $S$ is increased by the
dynamics since

$\displaystyle\dot{S}=(C^{\prime}(\vec{\Gamma}\cdot\vec{\omega}))^{2}\Gamma^{2}\omega^{2}\sin^{2}\theta\geq
0$ (43)

where $\theta$ measures the angle between $\vec{\omega}$ and $\vec{\Gamma}$.

## Dynamics of The Heavy Top

The first analytic function we try is $C(x)=\lambda x$. The equations of
motion are

$\displaystyle\dot{L}^{1}$
$\displaystyle=-L^{2}L^{3}\left(\frac{1}{I_{1}}-\frac{1}{I_{3}}\right)+\xi\Gamma^{2}+\lambda\left(\left[\frac{(L^{2})^{2}}{I_{1}^{2}}+\frac{(L^{3})^{2}}{I_{3}^{2}}\right]\Gamma^{1}-\frac{L^{1}L^{2}}{I^{2}_{1}}\Gamma^{2}-\frac{L^{1}L^{3}}{I_{1}I_{3}}\Gamma^{3}\right)$
(44) $\displaystyle\dot{L}^{2}$
$\displaystyle=L^{1}L^{3}\left(\frac{1}{I_{1}}-\frac{1}{I_{3}}\right)-\xi\Gamma^{1}+\lambda\left(-\frac{L^{1}L^{2}}{I^{2}_{1}}\Gamma^{1}+\left[\frac{(L^{1})^{2}}{I_{1}^{2}}+\frac{(L^{3})^{2}}{I_{3}^{2}}\right]\Gamma^{2}-\frac{L^{2}L^{3}}{I_{1}I_{3}}\Gamma^{3}\right)$
$\displaystyle\dot{L}^{3}$
$\displaystyle=\lambda\left(-\frac{L^{1}L^{3}}{I_{1}I_{3}}\Gamma^{1}-\frac{L^{2}L^{3}}{I_{1}I_{3}}\Gamma^{2}+\left[\frac{(L^{2})^{2}}{I_{1}^{2}}+\frac{(L^{1})^{2}}{I_{1}^{2}}\right]\Gamma^{3}\right)$
$\displaystyle\dot{\Gamma}^{1}$
$\displaystyle=\Gamma^{2}\frac{L^{3}}{I_{3}}-\Gamma^{3}\frac{L^{2}}{I_{1}}$
$\displaystyle\dot{\Gamma}^{2}$
$\displaystyle=\Gamma^{3}\frac{L^{1}}{I_{1}}-\Gamma^{1}\frac{L^{3}}{I_{3}}$
$\displaystyle\dot{\Gamma}^{3}$
$\displaystyle=\Gamma^{1}\frac{L^{2}}{I_{1}}-\Gamma^{2}\frac{L^{1}}{I_{1}}.$

It’s not hard to see that the metriplectic heavy top has an equilibrium when
$\Gamma^{1,2}=0$ and $L^{1,2}=0$. Given an initial configuration, we define
$L^{3}_{*}>0$ and $\Gamma_{*}^{3}>0$ to be the physically realizable
equilibrium consistent with $\Gamma^{1,2}=0$ and $L^{1,2}=0$. Linearizing the
equations of motion around this equilibrium, we get the equations

$\displaystyle\dot{\delta L}^{1}$
$\displaystyle=-L_{*}^{3}\left(\frac{1}{I_{1}}-\frac{1}{I_{3}}\right)\delta
L^{2}+\xi\delta\Gamma^{2}+\lambda\left(\frac{(L_{*}^{3})^{2}}{I_{3}^{2}}\delta\Gamma^{1}-\frac{L_{*}^{3}\Gamma_{*}^{3}}{I_{1}I_{3}}\delta
L^{1}\right)$ (45) $\displaystyle\dot{\delta L}^{2}$
$\displaystyle=L_{*}^{3}\left(\frac{1}{I_{1}}-\frac{1}{I_{3}}\right)\delta
L^{1}-\xi\delta\Gamma^{1}+\lambda\left(\frac{(L_{*}^{3})^{2}}{I_{3}^{2}}\delta\Gamma^{2}-\frac{L_{*}^{3}\Gamma_{*}^{3}}{I_{1}I_{3}}\delta
L^{2}\right)$ $\displaystyle\dot{\delta L}^{3}$ $\displaystyle=0$
$\displaystyle\dot{\delta\Gamma}^{1}$
$\displaystyle=\frac{L_{*}^{3}}{I_{3}}\delta\Gamma^{2}-\frac{\Gamma_{*}^{3}}{I_{1}}\delta
L^{2}$ $\displaystyle\dot{\delta\Gamma}^{2}$
$\displaystyle=\frac{\Gamma_{*}^{3}}{I_{1}}\delta
L^{1}-\frac{L_{*}^{3}}{I_{3}}\delta\Gamma^{1}$
$\displaystyle\dot{\delta\Gamma}^{3}$ $\displaystyle=0.$

We can express these equations compactly as

$\displaystyle\dot{\delta\mathbf{z}}=M\delta\mathbf{z}.$ (46)

The general solution for the spectrum of $M$ isn’t particularly enlightening
beyond the fact that the eigenvalues are always of the form
$(0,0,A,A^{*},B,B^{*})$. Depending on the initial conditions, either (a.)
$Re(A)$ and $Re(B)$ differ in sign or (b.) both $Re(A)<0$ and $Re(B)<0$. The
latter case is of particular interest since it implies the system is linearly
stable. If we impose that $I_{3}=2I_{1}$, the spectrum of $M$ is a lot easier
to understand with

$A=\frac{-\Gamma_{*}^{3}L_{*}^{3}\lambda-\sqrt{-4(I_{1})^{2}(L_{*}^{3})^{2}+(\Gamma_{*}^{3})^{2}(L_{*}^{3})^{2}\lambda^{2}+16\Gamma_{*}^{3}(I_{1})^{3}\xi}}{4(I_{1})^{2}}$
$B=\frac{-\Gamma_{*}^{3}L_{*}^{3}\lambda+\sqrt{-4(I_{1})^{2}(L_{*}^{3})^{2}+(\Gamma_{*}^{3})^{2}(L_{*}^{3})^{2}\lambda^{2}+16\Gamma_{*}^{3}(I_{1})^{3}\xi}}{4(I_{1})^{2}}.$

When the root is real, the system is linearly stable if and only if

$(L_{*}^{3})^{2}\geq 4\Gamma_{*}^{3}I_{1}\xi.$

When the root is complex, the system is always linearly stable, with the
dynamics oscillating towards the equilibrium $L^{1,2}=\Gamma^{1,2}=0$.

Computationally modeling the dynamics in Mathematica, it seems the relaxation
behavior of this system can be well understood in terms of its linearization,
provided the top does not fall. For example, in arbitrary units, we can let
$I_{1}=1=\frac{1}{2}I_{3}$, $\lambda=0.1$, and $\xi=1$. If we start with with
initial conditions such that $L_{*}^{3}\approx 5$, $\Gamma_{*}^{3}\approx 3$
then

$\displaystyle\text{Eigenvalues}(M)\approx(0,0,-.38+1.76i,-.38-1.76i,-.38+1.76i,-.38-1.76i).$
(47)

(a)

(b)

(c)

(d)

(e)

(f)

Figure 1: $\lambda=0.1;C(x)=\lambda x$

The linear dynamics predict that the metriplectic heavy top solutions will
decay towards equilibrium while oscillating. Modeling the full nonlinear
dynamics with initial conditions $\mathbf{\Gamma}(t=0)=(1,0,2.8)$ and
$\mathbf{L}(t=0)=(1,0,4.2)$, we see that the behavior of the system
qualitatively approximates the linear dynamics, even with large perturbations
from equilibrium (figure 1). This is an interesting point, since $S$ is not a
valid Lyapunov function. Likely, this stability comes from a sort of
constrained optimization on the surfaces of constant energy.

If we instead fix $\lambda=1$ then

$\displaystyle\text{Eigenvalues}(M)\approx(0,0,-7.03,-7.03,-.46,-.46),$ (48)

which is again stable, but this time we see no oscillatory behavior in the
linearized dynamics. Changing our initial conditions to
$\mathbf{\Gamma}(t=0)=(.3,0,3)^{T}$ and $\mathbf{L}(t=0)=(.5,0,5.2)$ so the
top doesn’t fall, we again observe the linear dynamics is qualitatively
similar to the relaxation behavior of the full nonlinear dynamics (figure 2).

(a)

(b)

(c)

(d)

(e)

(f)

Figure 2: $\lambda=1;C(x)=\lambda x$

Provided we change the phase space from $\mathbb{R}^{6}$ to
$\mathbb{R}^{6}-S$, where $S$ is the vanishing set of
$\mathbf{\Gamma}\cdot\mathbf{L}$, we can also consider the dynamics of the
heavy top with an otherwise singular generating function. In particular, we
explore when $S=\lambda\log(\mathbf{\Gamma}\cdot\mathbf{L})$. The equations of
motion are

$\displaystyle\dot{L}^{1}$
$\displaystyle=-L^{2}L^{3}\left(\frac{1}{I_{1}}-\frac{1}{I_{3}}\right)+\xi\Gamma^{2}+\frac{\lambda}{\mathbf{\Gamma}\cdot\mathbf{L}}\left(\left[\frac{(L^{2})^{2}}{I_{1}^{2}}+\frac{(L^{3})^{2}}{I_{3}^{2}}\right]\Gamma^{1}-\frac{L^{1}L^{2}}{I^{2}_{1}}\Gamma^{2}-\frac{L^{1}L^{3}}{I_{1}I_{3}}\Gamma^{3}\right)$
(49) $\displaystyle\dot{L}^{2}$
$\displaystyle=L^{1}L^{3}\left(\frac{1}{I_{1}}-\frac{1}{I_{3}}\right)-\xi\Gamma^{1}+\frac{\lambda}{\mathbf{\Gamma}\cdot\mathbf{L}}\left(-\frac{L^{1}L^{2}}{I^{2}_{1}}\Gamma^{1}+\left[\frac{(L^{1})^{2}}{I_{1}^{2}}+\frac{(L^{3})^{2}}{I_{3}^{2}}\right]\Gamma^{2}-\frac{L^{2}L^{3}}{I_{1}I_{3}}\Gamma^{3}\right)$
$\displaystyle\dot{L}^{3}$
$\displaystyle=\frac{\lambda}{\mathbf{\Gamma}\cdot\mathbf{L}}\left(-\frac{L^{1}L^{3}}{I_{1}I_{3}}\Gamma^{1}-\frac{L^{2}L^{3}}{I_{1}I_{3}}\Gamma^{2}+\left[\frac{(L^{2})^{2}}{I_{1}^{2}}+\frac{(L^{1})^{2}}{I_{1}^{2}}\right]\Gamma^{3}\right)$
$\displaystyle\dot{\Gamma}^{1}$
$\displaystyle=\Gamma^{2}\frac{L^{3}}{I_{3}}-\Gamma^{3}\frac{L^{2}}{I_{1}}$
$\displaystyle\dot{\Gamma}^{2}$
$\displaystyle=\Gamma^{3}\frac{L^{1}}{I_{1}}-\Gamma^{1}\frac{L^{3}}{I_{3}}$
$\displaystyle\dot{\Gamma}^{3}$
$\displaystyle=\Gamma^{1}\frac{L^{2}}{I_{1}}-\Gamma^{2}\frac{L^{1}}{I_{1}}.$

Again, this system has an equilibrium $\Gamma^{1,2}=0$ and $L^{1,2}=0$.
Letting be $L^{3}_{*}>0$ and $\Gamma_{*}^{3}>0$ be the realizable equilibrium
conditions, the linearized equations are

$\displaystyle\dot{\delta L}^{1}$
$\displaystyle=-L_{*}^{3}\left(\frac{1}{I_{1}}-\frac{1}{I_{3}}\right)\delta
L^{2}+\xi\delta\Gamma^{2}+\frac{\lambda}{L_{*}^{3}\Gamma_{*}^{3}}\left(\frac{(L_{*}^{3})^{2}}{I_{3}^{2}}\delta\Gamma^{1}-\frac{L_{*}^{3}\Gamma_{*}^{3}}{I_{1}I_{3}}\delta
L^{1}\right)$ (50) $\displaystyle\dot{\delta L}^{2}$
$\displaystyle=L_{*}^{3}\left(\frac{1}{I_{1}}-\frac{1}{I_{3}}\right)\delta
L^{1}-\xi\delta\Gamma^{1}+\frac{\lambda}{L_{*}^{3}\Gamma_{*}^{3}}\left(\frac{(L_{*}^{3})^{2}}{I_{3}^{2}}\delta\Gamma^{2}-\frac{L_{*}^{3}\Gamma_{*}^{3}}{I_{1}I_{3}}\delta
L^{2}\right)$ $\displaystyle\dot{\delta L}^{3}$ $\displaystyle=0$
$\displaystyle\dot{\delta\Gamma}^{1}$
$\displaystyle=\frac{L_{*}^{3}}{I_{3}}\delta\Gamma^{2}-\frac{\Gamma_{*}^{3}}{I_{1}}\delta
L^{2}$ $\displaystyle\dot{\delta\Gamma}^{2}$
$\displaystyle=\frac{\Gamma_{*}^{3}}{I_{1}}\delta
L^{1}-\frac{L_{*}^{3}}{I_{3}}\delta\Gamma^{1}$
$\displaystyle\dot{\delta\Gamma}^{3}$ $\displaystyle=0.$

Again, we may write

$\delta\dot{\mathbf{z}}=M\delta\mathbf{z}.$

Provided $2I_{1}=I_{3}$, the eigenvalues of $M$ are $(0,0,A,A^{*},B,B^{*})$
where

$A=\frac{-\Gamma_{*}^{3}\lambda-\sqrt{-4(\Gamma_{*}^{3})^{2}(I_{1})^{2}(L_{*}^{3})^{2}+(\Gamma_{*}^{3})^{2}\lambda^{2}+16\Gamma_{*}^{3}(I_{1})^{3}\xi}}{4(I_{1})^{2}\Gamma_{*}^{3}}$
$B=\frac{-\Gamma_{*}^{3}\lambda+\sqrt{-4(\Gamma_{*}^{3})^{2}(I_{1})^{2}(L_{*}^{3})^{2}+(\Gamma_{*}^{3})^{2}\lambda^{2}+16\Gamma_{*}^{3}(I_{1})^{3}\xi}}{4(I_{1})^{2}\Gamma_{*}^{3}}.$

(a)

(b)

(c)

(d)

(e)

(f)

Figure 3: $\lambda=1;C(x)=\lambda\log(x)$

The system is linearly stable when

$\Gamma_{*}^{3}(L_{*}^{3})^{2}>4I_{1}\xi$

For the sake of comparison to the $C(x)=\lambda x$ case, we use $I_{1}=1$,
$I_{1}=2$, $\xi=1$, and $\lambda=1$ with the initial conditions
$\Gamma=(.3,0,3)$ and $L=(.5,0,5.2)$. We see the nonlinear system stable
relaxes to $\Gamma^{3}_{*}\approx 3.0$ and $L^{3}_{*}\approx 5.2$ (figure 3).
This is reflected in the matrix for the linear equations of motion, which has
eigenvalues

$\text{Eigenvalues}(M)\approx(0,0,-0.25+2i,-0.25+2i,-0.25+2i,-0.25+2i).$

Also for the sake of comparison, we can also choose $\lambda=.1$,
$\mathbf{L}(t=0)=(1,0,4.2)$, and $\mathbf{\Gamma}(t=0)=(1,0,2.8)$. We see the
system relaxes to equilibrium, but at about a tenth the speed as before,
agreeing with the dynamics predicted by the eigenvalues of $M$ (figure 4)

$\text{Eigenvalue}(M)\approx(0,0,-0.025+1.67i,-0.025+1.67i,-0.025-1.67i,-0.025-1.67i).$

(a)

(b)

(c)

(d)

(e)

(f)

Figure 4: $\lambda=0.1;C(x)=\lambda\log(x)$

As a final example, we consider the case when $C(x)=\frac{\lambda}{2}x^{2}$.
The equations of motion are

$\displaystyle\dot{L}^{1}$
$\displaystyle=-L^{2}L^{3}\left(\frac{1}{I_{1}}-\frac{1}{I_{3}}\right)+\xi\Gamma^{2}+\lambda\mathbf{\Gamma}\cdot\mathbf{L}\left(\left[\frac{(L^{2})^{2}}{I_{1}^{2}}+\frac{(L^{3})^{2}}{I_{3}^{2}}\right]\Gamma^{1}-\frac{L^{1}L^{2}}{I^{2}_{1}}\Gamma^{2}-\frac{L^{1}L^{3}}{I_{1}I_{3}}\Gamma^{3}\right)$
(51) $\displaystyle\dot{L}^{2}$
$\displaystyle=L^{1}L^{3}\left(\frac{1}{I_{1}}-\frac{1}{I_{3}}\right)-\xi\Gamma^{1}+\lambda\mathbf{\Gamma}\cdot\mathbf{L}\left(-\frac{L^{1}L^{2}}{I^{2}_{1}}\Gamma^{1}+\left[\frac{(L^{1})^{2}}{I_{1}^{2}}+\frac{(L^{3})^{2}}{I_{3}^{2}}\right]\Gamma^{2}-\frac{L^{2}L^{3}}{I_{1}I_{3}}\Gamma^{3}\right)$
$\displaystyle\dot{L}^{3}$
$\displaystyle=\lambda\mathbf{\Gamma}\cdot\mathbf{L}\left(-\frac{L^{1}L^{3}}{I_{1}I_{3}}\Gamma^{1}-\frac{L^{2}L^{3}}{I_{1}I_{3}}\Gamma^{2}+\left[\frac{(L^{2})^{2}}{I_{1}^{2}}+\frac{(L^{1})^{2}}{I_{1}^{2}}\right]\Gamma^{3}\right)$
$\displaystyle\dot{\Gamma}^{1}$
$\displaystyle=\Gamma^{2}\frac{L^{3}}{I_{3}}-\Gamma^{3}\frac{L^{2}}{I_{1}}$
$\displaystyle\dot{\Gamma}^{2}$
$\displaystyle=\Gamma^{3}\frac{L^{1}}{I_{1}}-\Gamma^{1}\frac{L^{3}}{I_{3}}$
$\displaystyle\dot{\Gamma}^{3}$
$\displaystyle=\Gamma^{1}\frac{L^{2}}{I_{1}}-\Gamma^{2}\frac{L^{1}}{I_{1}},$

which we linearize

$\displaystyle\dot{\delta L}^{1}$
$\displaystyle=-L_{*}^{3}\left(\frac{1}{I_{1}}-\frac{1}{I_{3}}\right)\delta
L^{2}+\xi\delta\Gamma^{2}+\lambda
L_{*}^{3}\Gamma_{*}^{3}\left(\frac{(L_{*}^{3})^{2}}{I_{3}^{2}}\delta\Gamma^{1}-\frac{L_{*}^{3}\Gamma_{*}^{3}}{I_{1}I_{3}}\delta
L^{1}\right)$ (52) $\displaystyle\dot{\delta L}^{2}$
$\displaystyle=L_{*}^{3}\left(\frac{1}{I_{1}}-\frac{1}{I_{3}}\right)\delta
L^{1}-\xi\delta\Gamma^{1}+\lambda
L_{*}^{3}\Gamma_{*}^{3}\left(\frac{(L_{*}^{3})^{2}}{I_{3}^{2}}\delta\Gamma^{2}-\frac{L_{*}^{3}\Gamma_{*}^{3}}{I_{1}I_{3}}\delta
L^{2}\right)$ $\displaystyle\dot{\delta L}^{3}$ $\displaystyle=0$
$\displaystyle\dot{\delta\Gamma}^{1}$
$\displaystyle=\frac{L_{*}^{3}}{I_{3}}\delta\Gamma^{2}-\frac{\Gamma_{*}^{3}}{I_{1}}\delta
L^{2}$ $\displaystyle\dot{\delta\Gamma}^{2}$
$\displaystyle=\frac{\Gamma_{*}^{3}}{I_{1}}\delta
L^{1}-\frac{L_{*}^{3}}{I_{3}}\delta\Gamma^{1}$
$\displaystyle\dot{\delta\Gamma}^{3}$ $\displaystyle=0.$

Even with the usual caveat that $I_{3}=2I_{2}$, the eigenvalues of $M$ do not
have a particularly nice expression. We again let $\xi=1,I_{1}=1$, $I_{3}=2$,
and $\lambda=.1$. With the initial conditions $\mathbf{\Gamma}(t=0)=(1,0,2.8)$
and $\mathbf{L}(t=0)=(1,0,4.2)$ the system is linearly stable with

$\text{Eigenvalues}(M)\approx(-8.31,-8.31,-0.21,-0.21,0,0),$

as reflected by the nonlinear dynamics (figure 5).

It should always be noted, including in all the prior examples, that this
system is not stable for all initial conditions, as letting
$\mathbf{L}(t=0)=(1,0,2.8)^{T}$ and $\mathbf{\Gamma}(t=0)=(1,0,4.2)$
exemplifies (figure 6).

(a)

(b)

(c)

(d)

(e)

(f)

Figure 5: $\lambda=0.1$; $C(x)=\frac{\lambda}{2}x^{2}$;
$\mathbf{\Gamma}(t=0)=(1,0,2.8)$; $\mathbf{L}(t=0)=(1,0,4.2)$

(a)

(b)

(c)

(d)

(e)

(f)

Figure 6: $\lambda=0.1$; $C(x)=\frac{\lambda}{2}x^{2}$;
$\mathbf{L}(t=0)=(1,0,2.8)$; $\mathbf{\Gamma}(t=0)=(1,0,4.2)$

Keeping with tradition, we also try $\lambda=1$ with the initial conditions
$\Gamma=(.3,0,3)$ and $L=(.5,0,5.2)$. Like all the other systems, this initial
condition is linearly stable with the nonlinear dynamics following suit
(figure 7)

$\text{Eigevalues}(M)\approx(-11.85,-11.85,-0.317,-0.317,0,0).$

(a)

(b)

(c)

(d)

(e)

(f)

Figure 7: $\lambda=1$; $C(x)=\frac{\lambda}{2}x^{2}$

## Conclusion

In this paper, we showed how the formalism of Riemann-Poisson geometry can be
used to create dissipative dynamics. In particular, by considering the
simplest dissipative metric on the heavy top phase space, we constructed a
toy-model system that often relaxed to a stable equilibrium spinning about a
principle axis, extending the work by Matterasi and Morrison [4] to the case
where gravity cannot be neglected. We then computationally modeled the
dynamics of a few systems, which while possessing the same metriplectic
bracket, had their dynamics generated by different functions.

Beyond serving as a useful toy model and proof-of-concept for more
sophisticated work, the family of dynamical systems constructed in this paper
have an obvious control-theoretic utility. Without dissipating any energy,
beyond that required to monitor the orientation and rotation rate of a heavy
top, we showed how a torque can be applied in such a way as to align a
symmetric spinning body with its third moment of inertia. Perhaps more
importantly, this work highlights how the language of geometry can be fitted
to metriplectic dissipation, allowing for a wide class of geometric
constructions to be carried into the theory of dissipative systems.

This work can be extended in a number of ways. For one, it stands open to
exploration as to how different metrics can affect the allowed dynamics of
this and other systems. Another, perhaps harder question, is applying this
geometric formalism to field theories such as the Navier-Stokes equations. In
this paper, we also left many questions about the metriplectic heavy top
unanswered. Most notably, we did not fully address many interesting questions
related to the stability of the system, such as what regions of the phase
space are nonlinearly asymptotically stable.

## References

  * [1] B. Aliounea, M. Boucetta and A.. Lessiad  “On Riemann-Poisson Lie Groups”  In _Archivum Mathematicum_ 56.4, 2020, pp. 225–247  DOI: 10.5817/AM2020-4-225
  * [2] R. Fernandes  “Connection in Poisson Geometry: Holonomy and Invariants”  In _Differential Geometry_ 54.2, 2000, pp. 303–365  DOI: 10.4310/jdg/1214341648
  * [3] H. Goldstein, C. Poole and J. Safko  “Classical Mechanics”  Addison Wesley, 2001
  * [4] M. Materassi and P.. Morrison  “Metriplectic Torque for Rotation Control of A Rigid Body”  In _Cybernetics and Physics_ 7.2, 2018, pp. 78–86  DOI: 10.35470/2226-4116-2018-7-2-78-86
  * [5] J. Milnor  “Curvatures of Left Invariant Metrics on Lie Groups”  In _Advances in Mathematics_ 21.3, 1976, pp. 293–329  DOI: 10.1016/S0001-8708(76)80002-3
  * [6] P.. Morrison  “A Paradigm For Joined Hamiltonian and Dissipative Systems”  In _Physica D_ 18.1-3, 1986, pp. 410–419  DOI: 10.1016/0167-2789(86)90209-5
  * [7] P.. Morrison  “A Paradigm for Joined Hamiltonian and Dissipative Systems”  In _Physica D: Nonlinear Phenomena_ 18, 1986, pp. 410–419  DOI: 10.1016/0167-2789(86)90209-5
  * [8] P.. Morrison  “Hamiltonian Description of The Ideal Fluid”  In _Reviews of Modern Physics_ 70.2, 1998, pp. 467–519
  * [9] P.. Morrison  “Thoughts on brackets and dissipation: old and new”  In _Journal of Physics: Conference Series_ 169.012006, 2009  DOI: 10.1088/1742-6596/169/1/012006
  * [10] Philip J. Morrison and Michael H. Updike  “An inclusive curvature-like framework for describing dissipation: metriplectic 4-bracket dynamics”, 2023  arXiv:2306.06787 [math-ph]
  * [11] M. Nakahara  “Geometry, Topology, and Physics”, Graduate Student Series in Physics  Institute of Physics Publishing, 1981
  * [12] N. Sudarshan E.  “Classical Mechanics: A Modern Perspective”  Wspc, 2015
  * [13] J-L. Thiffeault and P.. Morrison  “Invariants and Labels for Lie-Poisson Systems”  In _Annals of the New York Academy of Sciences_ 867, 1998, pp. 109–119


* Cleynen and Lebarbier, (2014) Cleynen, A. and Lebarbier, E. (2014).  Segmentation of the Poisson and negative binomial rate models: a penalized estimator.  ESAIM Probab. Stat., 18, 750–769.
  * Cleynen and Lebarbier, (2017) Cleynen, A. and Lebarbier, E. (2017).  Model selection for the segmentation of multiparameter exponential family distributions.  Electron. J. Stat., 11, 800–842.
  * Fearnhead, (2006) Fearnhead, P. (2006).  Exact and efficient Bayesian inference for multiple changepoint.  Statist. Comput., 16, 203–213.
  * Fearnhead and Rigaill, (2019) Fearnhead, P. and Rigaill, G. (2019).  Changepoint detection in the presence of outliers.  J. Amer. Statist. Assoc., 114, 169–183.
  * Fearnhead and Rigaill, (2020) Fearnhead, P. and Rigaill, G. (2020).  Relating and comparing methods for detecting changes in mean.  Stat, e291.
  * Frick et al., (2013) Frick, K., Munk, A. and Sieling, H. (2013).  Multiscale change point inference.  J. Roy. Statist. Soc., Ser. B, 76, 495–580.
  * Fryźlewicz, (2014) Fryźlewicz, P. (2014).  Wild binary segmentation for multiple change-point detection.  Ann. Statist., 42, 2243–2281.
  * Gallagher et al., (2013) Gallagher, C., Lund, R. and Robbins, M. (2013).  Changepoint detection in climate time series with long-term trends.  J. Clim., 26, 4994–5006.
  * Green, (1995) Green, P. J. (1995).  Reversible jump Markov chain Monte Carlo computation and Bayesian model determination.  Biometrika, 82, 711–732.
  * Hotz et al., (2013) Hotz, T., Schütte, O. M., Sieling, H., Polupanow, T., Diederichsen, U., Steinem, C. and Munk, A. (2013).  Idealizing ion channel recordings by a jump segmentation multiresolution filter.  IEEE Trans. NanoBioscience, 12, 376–386.
  * Huang et al., (2005) Huang, T., Wu, B., Lizardi, P. and Zhao, H. (2005).  Detection of DNA copy number alterations using penalized least squares regression.  Bioinformatics, 21, 3811–3817.
  * Killick et al., (2012) Killick, R., Fearnhead, P. and Eckley, I. A. (2012).  Optimal detection of changepoints with a linear computational cost.  J. Amer. Statist. Assoc., 107, 1590–1598.
  * Kolaczyk and Nowak, (2005) Kolaczyk, E. D. and Nowak, R. D. (2005).  Multiscale generalised linear models for nonparametric function estimation.  Biometrika, 92, 119–133.
  * Le Cam, (1986) Le Cam, L. (1986).  Asymptotic Methods in Statistical Decision Theory.  Springer Series in Statistics. Springer, New York.
  * Le Cam and Yang, (1990) Le Cam, L. and Yang, G. L. (1990).  Asymptotics in Statistics : Some Basic Concepts.  Springer Series in Statistics. Springer, New York.
  * Li et al., (2016) Li, H., Munk, A. and Sieling, H. (2016)  FDR-control in multiscale change-point segmentation.  Electron. J. Stat., 10, 918–959.
  * Lloyd et al., (2015) Lloyd, C., Gunter, T., Osborne, M. A. and Roberts, S. J. (2015)  Variational inference for Gaussian process modulated Poisson processes.  In International Conference on Machine Learning, 37, 1814–1822.
  * Massart, (2007) Massart, P. (2007).  Concentration Inequalities and Model Selection, volume 1896 of Lecture Notes in Mathematics.  Springer, Berlin.  Lectures from the 33rd Summer School on Probability Theory held in Saint-Flour, July 6–23, 2003.
  * Muggeo and Adelfio, (2010) Muggeo, V. M. R. and Adelfio, G. (2010).  Efficient change point detection for genomic sequences of continuous measurements.  Bioinformatics, 27, 161–166.
  * Olshen et al., (2004) Olshen, A. B., Venkatraman, E., Lucito, R. and Wigler, M. (2004).  Circular binary segmentation for the analysis of array-based DNA copy number data.  Biostatistics, 5, 557–572.
  * Pierre-Jean et al., (2015) Pierre-Jean, M., Rigaill, G. and Neuvial, P. (2015).  Performance evaluation of DNA copy number segmentation methods.  Brief. Bioinformatics, 16, 600–615.
  * Redon et al., (2006) Redon, R., Ishikawa, S., Fitch, K., Feuk, L., Perry, G., Andrews, T., Fiegler, H., Shapero, M., Carson, A., Chen, W., Cho, E., Dallaire, S., Freeman, J., Gonzalez, J., Gratacòs, M., Huang, J., Kalaitzopoulos, D., Komura, D., Macdonald, J. and Hurles, M. (2006).  Global variation in copy number in the human genome.  Nature, 444, 444-454.
  * Reeves et al., (2007) Reeves, J., Chen, J., Wang, X. L., Lund, R. and Lu, Q. Q. (2007).  A review and comparison of changepoint detection techniques for climate data.  J. Appl. Meteorol. and Climatol., 46, 900–915.
  * Rigaill, (2015) Rigaill, G. (2015).  A pruned dynamic programming algorithm to recover the best segmentations with $1$ to $K_{max}$ change-points.  arXiv preprint, arXiv:1004.0887.
  * Scott and Knott, (1974) Scott, A. J. and Knott, M. (1974).  A cluster analysis method for grouping means in the analysis of variance.  Biometrics, 30, 507–512.
  * Spokoiny, (2009) Spokoiny, V. (2009).  Multiscale local change point detection with applications to value-at-risk.  Ann. Statist., 37, 1405–1436.
  * Tibshirani and Wang, (2007) Tibshirani, R. and Wang, P. (2007).  Spatial smoothing and hot spot detection for CGH data using the fused Lasso.  Biostatistics, 9, 18–29.
  * Truong et al., (2020) Truong, C., Oudre, L. and Vayatis, N. (2020).  Selective review of offline change point detection methods.  Signal Process., 167, 107299.
  * van der Vaart and Wellner, (1996) van der Vaart, A. W. and Wellner, J. A. (1996).  Weak Convergence and Empirical Processes. With Applications to Statistics.  Springer Series in Statistics. Springer-Verlag, New York.
  * Venkatraman and Olshen, (2007) Venkatraman, E. S. and Olshen, A. B. (2007).  A faster circular binary segmentation algorithm for the analysis of array CGH data.  Bioinformatics, 23, 657–663.
  * Verzelen et al., (2020) Verzelen, N., Fromont, M., Lerasle, M. and Reynaud-Bouret, P. (2020).  Optimal change-point detection and localization.  arXiv preprint, arXiv:2010.11470.
  * Wang et al., (2020) Wang, D., Yu, Y. and Rinaldo, A. (2020).  Univariate mean change point detection: penalization, CUSUM and optimality.  Electron. J. Stat., 14, 1917–1961.
  * Yang and Kuo, (2001) Yang, T. Y. and Kuo, L. (2001).  Bayesian binary segmentation procedure for a Poisson process with multiple changepoints.  J. Comput. Graph. Statist., 10, 772–785.
  * Zhang and Siegmund, (2007) Zhang, N. R. and Siegmund, D. O. (2007).  A modified Bayes information criterion with applications to the analysis of comparative genomic hybridization data.  Biometrics, 63, 22–32.


# A localization approach for autonomous underwater vehicles: A ROS-Gazebo
framework

Frederico C. Vaz1,2, David Portugal2, André Araújo2, Micael S. Couceiro1,2 and
Rui P. Rocha1 This work was supported by the Institute of Systems and Robotics
(ISR) – University of Coimbra funded by “Fundação para a Ciência e a
Tecnologia“ (FCT) under Grant No. UID/EEA/00048/2013.1F. Vaz, M. S. Couceiro
and R. P. Rocha are with Institute of Systems and Robotics, University of
Coimbra, 3030-290 Coimbra, Portugal
<EMAIL_ADDRESS>C. Vaz, D. Portugal, A.
Araújo and M. S. Couceiro are with Ingeniarius, Ltd., 3025-307 Coimbra,
Portugal<EMAIL_ADDRESS>
###### Abstract

Autonomous Underwater Vehicles (AUVs) have the ability to operate in harsh
underwater environments without endangering human lives in the process.
Nevertheless, just like their ground and aerial counterparts, AUVs need to be
able to estimate their own position. Yet, unlike ground and aerial robots,
estimating the pose of AUVs is very challenging, with only a few high-cost
technological solutions available in the market. In this paper, we present the
development of a realistic underwater acoustic model, implemented within the
Robot Operating System (ROS) and the Gazebo simulator framework, for
localization of AUVs using a set of water surface robots, time of flight of
underwater propagated acoustic waves, and a multilateration genetic algorithm
approach.

##  I INTRODUCTION

Autonomous Underwater Vehicles (AUVs) have been recommended for scientific
exploration and sub-aquatic inspection tasks. These robots need to face harsh
environmental conditions of oceans and rivers where visibility is reduced.
Many of the mobile robotics technologies are not suitable or do not work
correctly under these conditions. Moreover, the high cost of off-shore
missions makes it difficult to test hardware and software algorithms on real
scenarios. Robotics simulators are a practical and straightforward solution to
help developing and testing new ideas and prototypes on a virtual scenario for
localization, navigation and control algorithms. A robotic simulator for
underwater scenarios should be able to simulate robot motion, sensors and
environmental interactions of the AUVs. In addition, the assurance that the
developments on a robotic simulator can be realistically transfered to the
real robot scenario is an important prerequisite for robotics simulation
tools.

The focus in the present work is on acoustic localization for underwater ROS-
enabled robots. There are some high quality AUVs simulators available [1], but
they still do not cover adequately the AUVs localization using acoustic
beacons. The localization method proposed depends on the underwater acoustic
channel, therefore an important contribution of our work is on modeling
acoustic propagation in the water column. We make use of the Unmanned
Underwater Vehicle Simulator [2], a set of Gazebo plugins [3] and ROS nodes
under development, that are used for supporting and testing the simulation of
unnamed underwater vehicles (UUV), i.e AUVs. The UUV Simulator is based on
Gazebo, which provides a recognized physical fidelity, a programmatic
interface, sensor modeling and adequate documentation.

There are clear advantages of using the sound channel when compared to the
radio channel. Mainly, the radio channel in an underwater environment suffers
from large attenuation, resulting in very limited transmission range. The
underwater environment is a very complex medium for sound propagation,
specially in oceans, because of the heterogeneity and random fluctuations such
as the effect of surface seas and the ocean bottom variances. When the
acoustic signal propagates in the underwater environment there is a decrease
in the signal amplitude due to geometrical spreading, the distance traveled by
the signal, and the absorption due to the chemical properties of the sea or
river water. The range of acoustic waves are even more limited in high
frequencies due to the absorption factor. Additionally, the speed of acoustic
waves varies spatially in the sea water, due to variations in temperature and
pressure.

##  II Modeling the underwater acoustic channel

Figure 1: Underwater Acoustic Channel Software Architecture.

We modeled the acoustic channel by relating the local values of propagation
wave direction (Snell-Descartes law) and velocity. The modeling technique used
in our work is known as geometrical acoustics. Our method assumes that the
water column is divided into several layers that represent different depths
with distinct temperature, salinity and pH.

Figure 2: Localization System Architecture. Figure 3: The proposed system in a
simulation scenario: The AUV (1) and the 4 water surface robots (2).

In the first place, the method gets the source (beacon sensor) and receiver
(hydrophone sensor) positions to compute the ground truth distance. Then, the
method uses the environment conditions to compute the sound speed and the
absorption coefficient for each layer between the receiver and source depth.
Thus, the time-of-flight (TOF) between the source and the receiver takes into
account the distance and sound speed in each layer. The signal-to-noise-ratio
is the result from the environment noise and signal propagation that depends
on the absorption coefficient in each layer. The model architecture can be
seen in Fig. 1.

##  III Description of the system

The proposed system for AUV localization comprises the use of four water
surface robots as reference points, as shown in Fig. 3, and receptors of an
acoustic signal (through a hydrophone sensor) emitted by a beacon installed on
the AUV. With the time of flight (TOF) of the acoustic signal, we are able to
estimate the distances between each water surface robot to the AUV. Based on
the distance measurements to the AUV and on the water surface robots positions
given by the global position system (GPS), we developed a multilateration
genetic algorithm method to determine the 3-D position of the AUV. In this
way, the localization system will be tracking the position of the underwater
robot in earth-centered, earth-fixed (ECEF) to a East North Up (ENU) reference
frame. We used data from an inertial sensor (IMU) and a pressure sensor
mounted on the AUV to provide global orientation and reduce the localization
error of the AUV through sensor data fusion, using an Extended Kalman filter
[4]. The system architecture and the position given by the system can be seen
in Fig. 2 and in Fig. 4, respectively.

##  IV Conclusion and Future work

In this study, we develop an underwater localization system based on the UUV
Simulator [2]. Moreover, we have created a new plugin for Gazebo: the
Hydrophone-Beacon plugin, which provides a hydrophone sensor and a beacon as
actuator. The implementation of the underwater acoustic channel model allows
the UUV Simulator to provide an effective modeling and simulation tool for the
underwater robotics acoustic localization. This new Localization ROS package
enables testing underwater robot systems in a simulation environment before
going to real scenarios. Field experiments are foreseen in the close future,
in order to validate the system and the models that we use in the simulator
plugin.

Figure 4: AUV trajectory given by the system.

## References

  * [1] D. Cook, A. Vardy, and R. Lewis, “A survey of auv and robot simulators for multi-vehicle operations,” in _Autonomous Underwater Vehicles (AUV), 2014 IEEE/OES_. IEEE, 2014, pp. 1–8.
  * [2] M. M. M. Manhães, S. A. Scherer, M. Voss, L. R. Douat, and T. Rauschenbach, “UUV simulator: A gazebo-based package for underwater intervention and multi-robot simulation,” in _OCEANS 2016 MTS/IEEE Monterey_. IEEE, sep 2016\. [Online]. Available: https://doi.org/10.1109%2Foceans.2016.7761080
  * [3] N. Koenig and A. Howard, “Design and use paradigms for gazebo, an open-source multi-robot simulator,” in _Intelligent Robots and Systems, 2004.(IROS 2004). Proceedings. 2004 IEEE/RSJ International Conference on_ , vol. 3. IEEE, 2004, pp. 2149–2154.
  * [4] T. Moore and D. Stouch, “A generalized extended kalman filter implementation for the robot operating system,” in _Proceedings of the 13th International Conference on Intelligent Autonomous Systems (IAS-13)_. Springer, July 2014.