text:
\section{Contact and simulation}

\input{sections/06_particle_experiments}

\input{sections/07_system_experiments.tex}

\input{sections/08_discussion}

{
\small
\bibliographystyle{IEEEtran}

\section{Introduction}

\IEEEPARstart{C}{lassical} mechanical systems \cite{ClassicalMechanicsTaylor05} are the workhorse of controller design in geometric control\footnote{Geometric control encompasses operational space control \cite{khatib1987unified} and potential field methods \cite{KhatibPotentialFields1985}; we refer to them all here as geometric control.} \cite{bullo2004geometric}.
Their popularity derives from both powerful conservation properties \cite{LeonardSusskindClassicalMechanics} lending to simply Lyapunov stability and convergence analysis \cite{khalil2002nonlinear} and their link to Riemannian geometry (see geometric mechanics \cite{bullo2004geometric,ratliff2021geometry}), enabling the shaping of, for instance, straight-line end-effector behavior. However, classical mechanical systems have fundamentally limited expressivity.

When designing behavior in parts---such as combining target reaching with obstacle avoidance, joint limit avoidance, posturing, and speed regulation---these classical mechanical systems combine as metric-weighted averages of component policies (see Section~\ref{sec:ShortcomingsOfClassicalMechanics}). Three limitations are clear from this perspective:
(1) priority metrics depend only on position and not full system state; (2) the geometry is fundamentally linked to the metric making it impossible to design both independently; and (3) controllers, as a result, overly rely on potentials, dampers and unnatural task hierarchies.
Recent work in Riemannian Motion Policies (RMPs) \cite{cheng2018rmpflow},
demonstrated that removing the restrictions of classical mechanics can result in powerful design tools. However, RMPs sacrifice theoritcal insight for expressivity, forcing practitioners to gain experience before becoming proficient.

{\em Geometric fabrics}, a new physics of behavior, generalize classical mechanics while retaining the powerful mathematics of classical systems and promoting flexibility to circumvent the above issues. Geometric fabrics are what we call {\em bent} Finsler geometries \cite{ratliff2021geometry}---the Finsler system generalizes the Riemannian geometry of classical systems giving us velocity-dependent metrics, and the bending terms enable the independent shaping of associated policies.
Constraint forces are simple examples of bending terms commonly deployed to follow trajectories or trace surfaces \cite{Peters_AR_2008}, but our bending terms extend far beyond constraints. We derive fabrics that bend Finsler geometries to match any given nonlinear geometry using a technique called {\em geometry energization}.

One may view fabrics as a class of systems that instantaneously look like classical systems (with mass matrices, corresponding inertias, and forces), but which change from moment to moment as a function of state, similar to variable impedance systems; even concepts such as contact models generalize. In most cases, techniques for classical systems apply to fabrics, but now with the expressivity of RMPs.

Concretely, this work: (1) generalizes classical mechanics to capture the flexibility required for behavioral design; (2) develops an expressive class of RMPs that are stable, path consistent, and agnostic to parameterization (covariant); (3) develops a fabric component algebra analogous to the RMP algebra; and (4) presents experiments validating the theory and showing we can exploit the theoretical properties to out-perform a well-tuned implementation of RMPflow.

\iffalse
We start by reviewing background in Section~\ref{sec:Background} and discuss the limitations of classical systems.
We then define geometric fabrics in Section~\ref{sec:GeometricFabrics}, show that constrained fabrics are themselves fabrics, and present a simple stability and convergence theorem.
Subsection~\ref{sec:EnergizingHd2Geometries} presents energization, showing that fabrics are in a sense complete and can take on the shape of any geometry. Finally, we discuss computation in generalized coordinates in Subsection~\ref{sec:GeneralizedCoordinates} and use those results to define an algebra of {\em fabric components} in Section~\ref{sec:DesigningFabrics} which can be used to derive computational algorithms.
Section~\ref{sec:ControlledParticleExperiments} then analyzes empirically the advantages of geometric fabrics over standard unbent Riemannian and Finsler systems, and Section~\ref{sec:SystemExperiments} studies a full robot system implementation of fabrics, comparing it to a well-engineered and hardened implementation of RMPs.

The appendix and our website\footnote{\url{https://sites.google.com/nvidia.com/geometric-fabrics-behavior}} contains extended discussion, proofs, and additional experimental demonstrations.
\fi

\section{Related Work}
\label{sec:RelatedWork}
Second-order systems are natural for control since robots, themselves, are second-order systems (see \cite{IjspeertDMPs2013} for similar arguments), so we restrict our attention here to those.

Geometric control \cite{bullo2004geometric,khatib1987unified,Peters_AR_2008}, which uses energy shaping of classical mechanical systems for behavioral design, is the most directly related area.
We discuss the limitations of classical systems in Section~\ref{sec:ShortcomingsOfClassicalMechanics} and propose geometric fabrics to eliminate them. Energy shaping methods can leverage this broader class of systems for more expressive design.
Dynamic Movement Primitives (DMPs) \cite{IjspeertDMPs2013} are another closely related area. DMPs operate by perturbing stable systems with a forcing term that either vanishes or is periodic---the resulting system inherits its stability from the underlying stable system. This idea is again highly-compatible with geometric fabrics since we can swap out the underlying stable system with any stable second-order system, including geometric fabrics (see \cite{Wyk2021GeometriFabricsPhysicsArXivAbbrev} Appendix~\ref{sec:GeneralizingDmps} for details).
Optimal control \cite{OptimalControlEstimationStengel94} is another relevant area and again highly compatible. Any optimal controller optimizes over a dynamics function. Traditionally, that dynamics function is derived from classical mechanics, but fabrics are a viable generalized alternative that would enable building behavior into the dynamics to aid the optimization.

An alternative approach to fabrics, developed independent of this work, is the Pullback Bundle Dynamical System (PBDS), derived as a class of stable Riemannian systems on the tangent bundle \cite{BylardBonalliEtAl2021}.
PBDSs separate priority metric and policy design but in turn require task policies to be srictly Riemannian, a subset of the HD2 geometries of fabrics. Contact dynamics is also less clear, although such extensions are likely possible.
Another related area is Control Lyapunov Functions (CLF) \cite{Li2019LyapunovRmps} which constitute a general technique for stabilizing systems. CLFs work by projecting unstable systems onto a Lyapunov stable class introducing potentially large projection errors. With fabrics we directly design inherently stable nonlinear systems so the final system both captures the desired behavior and is stable by construction.
Finally, earlier work introduced Geometric Dynamical Systems (GDS) \cite{cheng2018rmpflow} as a stable class of RMPs. We now know that GDSs are better expressed as Finsler systems which are unbent fabrics. \cite{Wyk2021GeometriFabricsPhysicsArXivAbbrev} Appendix~\ref{apx:Gds} details this relation and Subsection~\ref{sec:ControlledParticleExperiments} shows empirically their limitations.

Other control approaches, including $\lambda_0$-PMP models \cite{tommasino2017task}, are first-order differential IK models related to classical Riemannian systems but lack the second-order fidelity and velocity dependence we address here. Fractal impedance methods \cite{tiseo2021passive} define locally optimal controllers that asymptotically track simultaneous position and force profiles. However, these are specialized controllers unlike geometric fabrics which describe a very general class of second order dynamical systems without any explicit notion of trajectory tracking or specific objectives. Moreover, composable energy policies \cite{urain2021composable} optimize actions over a product of stochastic policies to resolve multiple policy conflicts. Geometric fabrics address this issue as well through geometric policies. Alternative works try to further resolve multi-policy conflicts with hierarchical prioritization and null-space projections \cite{sentis2006whole, khatib1987unified,dietrich2018HierarchicalOpSpaceControl}. However, these structures are overly strict and relative priorities should instead continuously adjust.

\iffalse
Finally, Finsler and spray geometry in mathematics \cite{bucataru2007FinslerLagrangeDynamics} are, of course, highly relevant. That literature, however, remains insurmountably challenging for most roboticists. We work from the re-derivations presented in \cite{ratliff2021geometry}; it's an area of future work to lower the barrier-to-entry and leverage related ideas from the literature. To the best of our knowledge geometric fabrics and results on constrained fabrics and energization are novel, as are the fabric component algebra and resulting tools for putting them into practice.
\fi
\section{Background}
\label{sec:Background}
\subsection{Finsler and general nonlinear geometries}
\label{sec:FinslerAndHd2Geometries}

We briefly review some of the main ideas behind Finsler and general nonlinear geometries of paths (Homogeneous of Degree 2 (HD2) geometries)---see \cite{ratliff2021geometry} for details. We call any smooth differential equation ${\ddot{\q}} + \mathbf{h}(\mathbf{q}, {\dot{\q}}) = \mathbf{0}$ an HD2 geometry when $\mathbf{h}$ is HD2 in velocities so that $\mathbf{h}(\mathbf{q}, \alpha{\dot{\q}}) = \alpha^2\mathbf{h}(\mathbf{q}, {\dot{\q}})$. We often write it in resolved form ${\ddot{\q}} = -\mathbf{h}(\mathbf{q}, {\dot{\q}}) = \pi(\mathbf{q}, {\dot{\q}})$ and refer to $\pi(\mathbf{q}, {\dot{\q}})$ as a {\em geometric policy} to give policy semantics. The HD2 property enforces that integral curves are \textit{path consistent}, i.e., they follow geometric speed-independent paths \cite{ratliff2021geometry}.
{\em Finsler geometries} are HD2 geometries defined by their minimum length geodesics using the calculus of variations. Their generalized arc-length measures are defined by a Finsler structure $\mathcal{L}_g(\mathbf{q},{\dot{\q}})$:

\begin{definition} \label{def:FinslerStructure}
A Finsler structure is a Lagrangian with the following three properties:
\begin{enumerate}
    \item
    $\mathcal{L}_g(\mathbf{q},{\dot{\q}})\geq 0$
    with equality if and only if ${\dot{\q}}=\mathbf{0}$.
    \item $\mathcal{L}_g$ is positive homogeneous (HD1) in ${\dot{\q}}$ so that $\mathcal{L}_g(\mathbf{q},\alpha{\dot{\q}}) = \alpha\mathcal{L}_g(\mathbf{q}, {\dot{\q}})$ for $\alpha \geq 0$.
    \item $\partial^2_{{\dot{\q}}\qd}\mathcal{L}_e\succeq 0$ when ${\dot{\q}}\neq\mathbf{0}$ where $\mathcal{L}_e = \frac{1}{2}\mathcal{L}_g^2$.
\end{enumerate}
\end{definition}
The first two properties enforce that $\mathcal{L}_g(\mathbf{q}, {\dot{\q}})$ is a speed-independent measure of length (see \cite{ratliff2021geometry}). The third property ensures the equations of motion of the associated {\em Finsler energy} $\mathcal{L}_e$, derived from the Euler-Lagrange equation as $\frac{d}{dt}\partial_{\dot{\q}}\mathcal{L}_e - \partial_\mathbf{q}\mathcal{L}_e = \mathbf{M}(\mathbf{q},{\dot{\q}}){\ddot{\q}} + \bm{\xi}(\mathbf{q}, {\dot{\q}}) = \mathbf{0}$, with $\mathbf{M} = \partial_{{\dot{\q}}\qd}^2\mathcal{L}_e$ (mass / metric) and $\bm{\xi} = \partial_{{\dot{\q}}\mathbf{q}}\mathcal{L}_e{\dot{\q}} - \partial_\mathbf{q}\mathcal{L}_e$ (geometric terms), are well-formed with positive definite mass. These (closed) equations of motion (absent potential, damping, or external forces) define an HD2 (Finsler) geometry \cite{ratliff2021geometry}.

A simple example of a Finsler geometry is Riemannian geometry. The Riemannian Finsler structure measures length $\mathcal{L}_g(\mathbf{q}, {\dot{\q}}) = \big({\dot{\q}}^T\mathbf{M}(\mathbf{q}){\dot{\q}}\big)^{\frac{1}{2}} = \|{\dot{\q}}\|_\mathbf{M}$ using a symmetric positive definite metric $\mathbf{M}(\mathbf{q})$ dependent on position only. The corresponding energy is the well-known classical kinetic energy $\mathcal{L}_e = \frac{1}{2}\mathcal{L}_g^2 = \frac{1}{2}\|{\dot{\q}}\|_\mathbf{M}^2 = \frac{1}{2}{\dot{\q}}^T\mathbf{M}(\mathbf{q}){\dot{\q}}$ for mass $\mathbf{M}$. This link between geometry and system energy is reflected more generally in the link between Finsler geometry and Finsler energy, with equations of motion tracing geodesics.

Every Lagrangian system has an associated conserved energy (its Hamiltonian) $\mathcal{H} = \partial_{\dot{\q}}\mathcal{L}^\tr{\dot{\q}} - \mathcal{L} = \mathbf{p}^\tr{\dot{\q}} - \mathcal{L}$ where $\mathbf{p} = \partial_{\dot{\q}}\mathcal{L}$ is the {\em generalized momentum}. When the Lagrangian is a Finsler energy $\mathcal{L}_e = \frac{1}{2}\mathcal{L}_g^2$ with equations of motion $\mathbf{M}(\mathbf{q}, {\dot{\q}}){\ddot{\q}} + \bm{\xi}(\mathbf{q},{\dot{\q}}) = \mathbf{0}$, the Hamiltonian can be expressed $\mathcal{H}_e = \mathcal{L}_e = \frac{1}{2}{\dot{\q}}^\tr\mathbf{M}{\dot{\q}}$, showing that $\mathcal{L}_e$ itself is the energy and can be expressed in the classical form \cite{ratliff2021geometry}.

The Finsler metric $\mathbf{M}(\mathbf{q}, {\dot{\q}})$ is Homogeneous of Degree 0 (HD0) which means $\mathbf{M}(\mathbf{q}, \alpha{\dot{\q}}) = \mathbf{M}(\mathbf{q}, {\dot{\q}})$. $\mathbf{M}$ is, therefore, dependent on  velocity direction $\widehat{{\dot{\q}}}$ but not on speed. That means, for any system $\mathbf{M}(\mathbf{q},{\dot{\q}}){\ddot{\q}} + \mathbf{f}(\mathbf{q},{\dot{\q}}) = \mathbf{0}$ (where $\mathbf{f}$ may capture more than just $\bm{\xi}$), the resolved equation ${\ddot{\q}} = -\mathbf{M}^{-1}\mathbf{f}$ defines a geometric policy if and only if $\mathbf{f}$ is HD2.

\subsection{Energy Shaping with Classical Mechanical Systems}
\label{sec:ShortcomingsOfClassicalMechanics}

Classical mechanical systems have traditionally been important tools for energy shaping in controller design \cite{bullo2004geometric}.
The main idea of energy shaping is to design a virtual classical mechanical system that captures a desired behavior, then control the physical system to behave like it thereby inheriting the stability of the virtual system. This idea is powerful which is why it has been so important and popular in control. But we show here that classical systems have limited expressivity as tools for behavioral design.

Let $\mathcal{Q}$ be a generalized coordinate domain constrained by $\Cbb(\mathbf{q})=\mathbf{0}$ with Jacobian $\mathbf{J}(\mathbf{q})$.
A constrained classical mechanical system evolves according to the equation \cite{ClassicalMechanicsTaylor05}:
\begin{align} \nonumber
    \mathbf{M}(\mathbf{q}){\ddot{\q}} + \bm{\xi}(\mathbf{q}, {\dot{\q}}) + \mathbf{J}(\mathbf{q})^\tr\bm{\lambda}(\mathbf{q},{\dot{\q}}) = -\partial_\mathbf{q}\psi(\mathbf{q}) - \mathbf{B}(\mathbf{q},{\dot{\q}}){\dot{\q}},
\end{align}
where $\mathbf{M}(\mathbf{q}){\ddot{\q}} + \bm{\xi}(\mathbf{q}, {\dot{\q}}) = \mathbf{0}$ form the equations of motion of an unforced system with kinetic energy $\mathcal{K}(\mathbf{q}, {\dot{\q}}) = \frac{1}{2}{\dot{\q}}^\tr\mathbf{M}(\mathbf{q}){\dot{\q}}$ (the term $\bm{\xi}(\mathbf{q},{\dot{\q}})$ captures fictitious forces which we call here {\em geometric forces} in reference to their role in Finsler geometry), the function $\psi(\mathbf{q})$ is a potential function, $\mathbf{B}(\mathbf{q},{\dot{\q}})$ is a positive definite damping matrix, and $\bm{\lambda}(\mathbf{q},{\dot{\q}})$ are the Lagrange multipliers of the constraint force. With $\mathbf{B}=\mathbf{0}$ the system conserves total energy $\mathcal{E}(\mathbf{q},{\dot{\q}}) = \mathcal{K}(\mathbf{q},{\dot{\q}}) + \psi(\mathbf{q})$; with nonzero $\mathbf{B}$ the energy dissipates, and the system converges to a local minimum of the potential on the constrained domain and is therefore inherently stable.

Unfortunately, a simple analysis shows that these systems have limited expressivity. Consider multiple classical systems $\mathbf{M}_i(\mathbf{q}){\ddot{\q}} + \mathbf{f}_i(\mathbf{q},{\dot{\q}}) = \bm{\tau}_i$, where now $\mathbf{f}_i$ captures all forces excluding constraints (fictitious, potential, and damping---constraint forces are often added to the combined system after summing, so we ignore them here) and we single out $\bm{\tau}_i$ as the force the other systems place on system $i$. For instance, these component systems may express reaching a target, avoiding obstacles, joint limit avoidance, and posturing of the arm. The internal forces $\bm{\tau}_i$ must sum to zero $\sum_i\bm{\tau}_i = \mathbf{0}$, so we have $\sum_i\big(\mathbf{M}_i{\ddot{\q}} + \mathbf{f}_i\big) = \sum_i\bm{\tau}_i = \mathbf{0}$. Assuming each $\mathbf{M}_i$ is full rank, we can write $\pi_i = -\mathbf{M}_i^{-1}\mathbf{f}_i$ and express the resulting combined acceleration
\iffalse
as\footnote{A similar result holds if $\mathbf{M}_i$ is not full rank. Intuitively, we can add a small $\epsilon$ to the diagonal then take the limit $\epsilon\rightarrow 0$ in Equation~\ref{eqn:ClassicalMetricWeightedAverage}.}
\else
as
\fi
\begin{align} \label{eqn:ClassicalMetricWeightedAverage}
    {\ddot{\q}} = \left(\sum_i\mathbf{M}_i\right)^{-1}\sum_i\mathbf{M}_i\pi_i,
\end{align}
which is simply a metric-weighted average of constituent policies $\pi_i$, where each policy takes the form:
\begin{align} \nonumber
    \pi_i(\mathbf{q}, {\dot{\q}}) = -\mathbf{M}_i^{-1}\Big(\bm{\xi}_i(\mathbf{q},{\dot{\q}}) + \partial_\mathbf{q}\psi_i(\mathbf{q}) + \mathbf{B}_i(\mathbf{q},{\dot{\q}}){\dot{\q}}\Big).
\end{align}
Each component policy is defined by geometric forces, its potential function, and its damping matrix; they averaged together using the mass matrices as spectral priority weights (with Eigenvalues expressing different priorities along each Eigenvector). Three limitations manifest from this analysis.

First, priority metrics are functions of position only and not full position-velocity state. Many tasks are best expressed with velocity aware priorities which classical metrics cannot represent. E.g. any barrier avoidance task (obstacle, self, joint limit avoidance) should care about the barrier while moving toward it, but should forget about it when moving away.

Second, geometric forces can be powerful---they define the Riemannian geometry of the system associated with the metric $\mathbf{M}(\mathbf{q})$ \cite{bullo2004geometric,ratliff2021geometry}---but in classical systems they are intrinsically linked to the metric which also plays the role of the policy's priority matrix. Only one of the two can be designed, the other becomes an artifact.
This coupling results in fighting between components (when the metric is compromised) or forces an over-reliance on potentials and dampers (when the geometry is compromised). An example of a Riemannian geometry working against the policy is the robot's own (unshaped) inertia which will fight a potential.
Likewise, straight-line end-effector motion is a useful shaped geometry, but that defines a largely irrelevant priority metric.

Third, when the geometric forces cannot be used, we must shape the behavior entirely using potential functions and dampers. Potential functions add to the total potential energy affecting where local minima fall resulting in task conflict.
Moreover, potentials are only position-dependent and can induce spring-like oscillations when damping is insufficient. Alternatively, one can increase damping, or even entirely shape the behavior using dampers, but that can slow the system when mutliple contributions combine. Geometric forces, on the other hand, are path consistent and will shape the behavior without shifting the system local minima.

\iffalse
Note that the second and third limitations are largely the motivation for hierarchical task prioritization \cite{dietrich2018HierarchicalOpSpaceControl}---when multiple potentials may conflict or priority metrics are compromised, designers preserve task fidelity by instituting a task hierarchy whereby lower-priority tasks are projected into the null-space of higher-priority tasks \cite{sentis2006whole, khatib1987unified}. However, these task hierarchies are often pragmatic solutions to these fundamental problems and somewhat artificial in practice. For instance, the relative priorities of reaching targets and avoiding obstacles should adjust continuously; it is unclear which is globally higher priority.
\fi

Geometric fabrics remove these limitations by (1) introducing velocity-dependent metrics using Finsler geometries, and (2) decoupling system geometry from priority metric using bending terms. Now, the system's behavior is primarily directed by the fabric, the potential is applied for a specific task, and damping for convergence and speed control.
\section{Geometric Fabrics: A New Physics of Behavior}

Geometric fabrics are a strict generalization of classical mechanics (with classical systems being a type of fabric) which retain the powerful mathematical properties exploited in control while removing the expressivity limitations outlined in Subsection~\ref{sec:ShortcomingsOfClassicalMechanics}. In particular, fabrics' priority metrics depend on full state and are decoupled from their associated geometric policy. Stability of these systems is easily analyzed using energy conservation (see Theorem~\ref{thm:GeometricFabricStability}). We derive fabrics as bent Finsler geometries using what we call {\em bending terms} which generalize constraint forces. These bending terms are powerful and can be used to bend a Finsler geometry to align with any desired HD2 geometry using a process we call geometry {\em energization}. We additionally derive formulas for calculation in generalized coordinates which we use in Section~\ref{sec:DesigningFabrics} to construct an algebra of fabric components to aid systems engineering.

\subsection{Bending Finsler geometries}
\label{sec:BendingFinslerGeometries}

An important concept in geometric fabrics is the {\em bending} of a Finsler geometry. Geometric bending is implemented using system modification terms that both perform no work (conserve energy) and are geometric (HD2).
An example of such a bending term is the classical constraint force that keeps a system moving along a constraint surface. These forces are zero-work by design (see the Principal of Virtual Work \cite{ClassicalMechanicsTaylor05}) and are known to preserve the geometric character of the system. Bending terms, as defined and characterized in the following Lemma, generalize classical constraints.

\begin{lemma} \label{lma:EnergyConservation}
Let $\mathcal{L}_e$ be a Finsler energy with unforced equations of motion $\mathbf{M}(\mathbf{q},{\dot{\q}}){\ddot{\q}} + \bm{\xi}(\mathbf{q},{\dot{\q}}) = \mathbf{0}$, and let $\psi(\mathbf{q})$ be a potential function and $\mathbf{B}(\mathbf{q}, {\dot{\q}})$ be a positive definite damping matrix. The modified system $\mathbf{M}{\ddot{\q}} + \bm{\xi} + \mathbf{f}_f = -\partial\psi - \mathbf{B}{\dot{\q}}$ with modifying term $\mathbf{f}_f(\mathbf{q},{\dot{\q}})$ has total energy $\mathcal{H} = \mathcal{L}_e + \psi$ that decreases at a rate
\begin{align}
    \dot{\mathcal{H}} = -{\dot{\q}}^\tr\mathbf{B}{\dot{\q}}
\end{align}
if and only if ${\dot{\q}}^\tr\mathbf{f}_f = \mathbf{0}$ everywhere. We call such a term $\mathbf{f}_f$ a {\em zero-work modification}. When $\mathbf{f}_f$ is additionally HD2, the unforced system (with $\psi=0$ and $\mathbf{B}=\mathbf{0}$) is geometric and we say $\mathbf{f}_f$ is a {\em bending term} that bends the Finsler geometry. Moreover, if $\psi$ is lower bounded, $\mathcal{H}$ is lower bounded.
\end{lemma}
\begin{proof}
See \cite{Wyk2021GeometriFabricsPhysicsArXivAbbrev} Appendix~\ref{apx:LemmaEnergyConservation}.
\end{proof}
Bending terms are geometric forces that are everywhere orthogonal to the direction of motion and thereby conserve energy by performing no work.

\subsection{Geometric fabrics}
\label{sec:GeometricFabrics}

Geometric fabrics generalize classical systems in two steps: (1) they generalize the classical Riemannian geometries to the broader class of Finsler geometries; and (2) they bend the resulting Finsler geometry using bending terms which include, but are not limited to, constraint forces.

\begin{definition}[Geometric fabric]
Let $\mathcal{L}_e = \frac{1}{2}\mathcal{L}_g^2$ be a Finsler energy derived from Finsler structure $\mathcal{L}_g$
with equations of motion $\mathbf{M}(\mathbf{q}, {\dot{\q}}){\ddot{\q}} + \bm{\xi}(\mathbf{q}, {\dot{\q}}) = \mathbf{0}$.
A geometric fabric is a system that evolves according to
\begin{align} \label{eqn:GeometricFabricEquationsOfMotion}
    &\mathbf{M}(\mathbf{q}, {\dot{\q}}){\ddot{\q}} + \bm{\xi}(\mathbf{q}, {\dot{\q}}) + \mathbf{f}_f(\mathbf{q}, {\dot{\q}}) = \mathbf{0},
\end{align}
where $\mathbf{f}_f(\mathbf{q}, {\dot{\q}})$ is any bending term. Geometric fabrics are denoted $\mathcal{F} = (\mathcal{L}_e, \mathbf{f}_f)$.
\end{definition}
The fabric defines a nominal behavior independent of a specific task. Potential functions then push away from that nominal behavior to drive the system toward task goals, and damping dissipation ensures convergence. The full equations of motion of the system then
\iffalse
become:\footnote{A broader theory of {\em optimization} fabrics \cite{optimizationFabricsForBehavioralDesignArXiv2020} drops the geometric requirement but retains the optimization properties outlined in Theorem~\ref{thm:GeometricFabricStability}. Geometric fabrics are a type of fabric that most naturally generalizes classical mechanics. We refer to geometric fabrics simply as fabrics.}
\else
become:
\fi
\begin{align}\label{eqn:ForcedFabricEquationsOfMotion}
    &\mathbf{M}(\mathbf{q}, {\dot{\q}}){\ddot{\q}} + \bm{\xi}(\mathbf{q}, {\dot{\q}}) + \mathbf{f}_f(\mathbf{q}, {\dot{\q}}) + \mathbf{J}(\mathbf{q})^\tr\bm{\lambda}(\mathbf{q},{\dot{\q}}) \\\nonumber
    &\ \ \ \ = -\partial_\mathbf{q}\psi(\mathbf{q}) - \mathbf{B}(\mathbf{q}, {\dot{\q}}){\dot{\q}},
\end{align}
where $\psi(\mathbf{q})$ is a potential function, $\mathbf{B}(\mathbf{q}, {\dot{\q}})$ is a positive definite damping matrix, and $\bm{\lambda}$ are the Lagrange multipliers of constraint forces with constraint Jacobian $\mathbf{J}$. We see many similarities to the classical equations given in Subsection~\ref{sec:ShortcomingsOfClassicalMechanics}. The potential function injects energy into the system and drives it toward its local minimum, while the damper bleeds energy off for stability. Lemma~\ref{lma:ConstrainedGeometricFabrics} characterizes constrained geometric fabrics and shows that a constrained geometric fabric is itself a geometric fabric. Theorem~\ref{thm:GeometricFabricStability} proves stability of these systems.

\begin{lemma} \label{lma:ConstrainedGeometricFabrics}
Let $(\mathcal{L}_e, \mathbf{f}_f)$ be a geometric fabric with Finsler energy $\mathcal{L}_e$, bending term $\mathbf{f}_f$, and associated equations of motion $\mathbf{M}{\ddot{\q}} + \bm{\xi} + \mathbf{f}_f = -\partial\psi - \mathbf{B}{\dot{\q}}$ for potential $\psi$ and damping matrix $\mathbf{B}$. Under equality constraint $\Cbb(\mathbf{q}) = \mathbf{0}$, assuming $(\mathbf{q},{\dot{\q}})$ already satisfy the constraint, the resulting constrained equations of motion are
\begin{align}\label{eqn:ConstrainedGeometricFabrics}
    \mathbf{M}{\ddot{\q}} + \mathbf{P}_{\!/\mkern-5mu/\!}\big(\bm{\xi}+\mathbf{f}_f\big) + \mathbf{P}_\perp \mathbf{M}{\dot{\J}}{\dot{\q}}
    = -\mathbf{P}_{\!/\mkern-5mu/\!}\big(\partial\psi + \mathbf{B}{\dot{\q}}\big)
\end{align}
where, $\mathbf{J} = \partial_\mathbf{q}\Cbb$ denotes the constraint Jacobian, $\mathbf{P}_\perp = \mathbf{J}^\tr(\mathbf{J}\mathbf{M}^{-1}\mathbf{J}^\tr)^{-1}\mathbf{J}\mathbf{M}^{-1}$ projects orthogonally to the constraints, and $\mathbf{P}_{\!/\mkern-5mu/\!} = \mathbf{I} - \mathbf{P}_\perp$ projects parallel to the constraints. Importantly, the left hand side defines a fabric $(\mathcal{L}_e, \wt{\mathbf{f}}_f)$ where $\wt{\mathbf{f}}_f = \mathbf{P}_\perp(\mathbf{M}{\dot{\J}}{\dot{\q}} - \bm{\xi}) + \mathbf{P}_{\!/\mkern-5mu/\!}\mathbf{f}_f$ is a {\em constrained} fabric.
\end{lemma}
\begin{proof}
See \cite{Wyk2021GeometriFabricsPhysicsArXivAbbrev} Appendix~\ref{apx:LemmaConstrainedGeometricFabrics}.
\end{proof}

The constrained system in Equation~\ref{eqn:ConstrainedGeometricFabrics} is intuitive. $\mathbf{P}_{\!/\mkern-5mu/\!}\big(\bm{\xi}+\mathbf{f}_f\big)$ projects the original fabric to operate parallel to the constraint and $-\mathbf{P}_{\!/\mkern-5mu/\!}(\partial\psi+\mathbf{B}{\dot{\q}})$ does the same with the forcing terms. $\mathbf{P}_\perp\mathbf{M}{\dot{\J}}{\dot{\q}}$ independently creates the required curvature force to keep the system moving along the constraint, modifying the original fabric's bending term.

\subsection{Stability analysis}
\label{sec:StabilityAnalysis}

The following theorem presents a general result on the stability and convergence of forced geometric fabrics.
\begin{theorem}[Geometric fabric stability] \label{thm:GeometricFabricStability}
Suppose $\mathcal{F} = (\mathcal{L}_e, \mathbf{f}_f)$ is geometric fabric operating under equality constraints $\mathbf{c}(\mathbf{q})=\mathbf{0}$ and inequality constraints $\mathbf{g}(\mathbf{q})\geq\mathbf{0}$ (both of which may be empty). Then for any $\psi(\mathbf{q})$ lower-bounded on the constrained domain and any bounded positive-definite damping matrix $\mathbf{B}(\mathbf{q}, {\dot{\q}})$ the system evolving according to Equation~\ref{eqn:ForcedFabricEquationsOfMotion} (with $\Cbb$ consisting of $\mathbf{c}$ and the active subset of $\mathbf{g}$)
converges to a local minimum of the optimization problem
\begin{align}
    &\min_{\mathbf{q}\in\mathcal{Q}} \psi(\mathbf{q})
    \ \ \mathrm{s.t.}\ \ \mathbf{c}(\mathbf{q}) = \mathbf{0}\ \ \mathrm{and}\ \ \mathbf{g}(\mathbf{q})\geq\mathbf{0}
\end{align}
under any contact model contributing as equality constraints with dissipation when in contact and conserving or dissipating energy on impact.
\end{theorem}
\begin{proof}
See \cite{Wyk2021GeometriFabricsPhysicsArXivAbbrev} Appendix~\ref{apx:TheoremGeometricFabricStability}
\end{proof}

\subsection{Completeness: Energizing HD2 geometries}
\label{sec:EnergizingHd2Geometries}

The next theorem is central to the theory of fabrics. It shows fabrics are complete in the sense that they can be bent to align with any given geometry. Completeness shows we have independent control over priority metric design (Finsler metric) and geometric policy design (desired geometry).
\begin{theorem}[Energized fabrics] \label{thm:SystemEnergization}
Let ${\ddot{\q}} = \pi(\mathbf{q}, {\dot{\q}})$ be an HD2 geometry, and suppose $\mathcal{L}_e$ is a Finsler energy with equations of motion $\mathbf{M}{\ddot{\q}} + \bm{\xi} = \mathbf{0}$ and energy $\mathcal{H}_e=\mathcal{L}_e$. Then ${\ddot{\q}} = \pi(\mathbf{q}, {\dot{\q}}) + \alpha_{\mathcal{H}_e}{\dot{\q}}$ is energy conserving when
\begin{align}\label{eqn:EnergizationTransformAlpha}
    \alpha_{\mathcal{H}_e} = -({\dot{\q}}^\tr\mathbf{M}{\dot{\q}})^{-1}{\dot{\q}}^\tr\big[\mathbf{M}\pi + \bm{\xi}\big],
\end{align}
and differs from the original system by only an acceleration along the direction of motion.
The new system is a geometric fabric, known as an {\em energized fabric}, expressed as:
\begin{align} \label{eqn:ZeroWorkEnergizationForm}
  \mathbf{M}{\ddot{\q}} + \bm{\xi} - \mathbf{P}_e\big[\mathbf{M}\pi + \bm{\xi}\big] = \mathbf{0},
\end{align}
where $\mathbf{P}_e = \mathbf{M}^{\,\frac{1}{2}}\Big[\mathbf{I} - \hat{\mathbf{v}} \hat{\mathbf{v}}^\tr\Big]\mathbf{M}^{-\frac{1}{2}}$ is a metric-weighted projection matrix
with $\mathbf{v} = \mathbf{M}^{\,\frac{1}{2}}{\dot{\x}}$ and $\hat{\mathbf{v}} = \mathbf{v}/\|\mathbf{v}\|$.
\end{theorem}
\begin{proof}
See \cite{Wyk2021GeometriFabricsPhysicsArXivAbbrev} Appendix~\ref{apx:TheoremOnEnergizedFabrics}.
\end{proof}
We denote the energization operation using an operator $\mathrm{energize}_{\mathcal{L}_e}[\cdot]$ that outputs a geometric fabric given an HD2 geometry $\pi$ so that
\begin{align}
    \mathrm{energize}_{\mathcal{L}_e}[\pi] = \Big(\mathcal{L}_e,\: -\mathbf{P}_e\big[\mathbf{M}\pi+\bm{\xi}\big]\Big)
\end{align}
is the energized geometry. Note that the energization operator can be applied to a non-geometric policy as well. In that case, we write the output in the same form using the same energization coefficient calculation in Equation~\ref{eqn:EnergizationTransformAlpha}, but with the knowledge that the resulting transformed system will conserve energy but not be a geometric fabric.

An informative example is when $\mathcal{L}_e = \frac{1}{2}\|{\dot{\q}}\|^2$ measures the squared velocity norm. Here $\mathbf{M} = \mathbf{I}$ and $\bm{\xi} = \mathbf{0}$, so under Equation~\ref{eqn:EnergizationTransformAlpha} the system reduces to
\begin{align}
    {\ddot{\q}} = \left[\mathbf{I} - \widehat{{\dot{\q}}}\,\widehat{{\dot{\q}}}^\tr\right]\pi(\mathbf{q}, {\dot{\q}}),
\end{align}
where $\widehat{{\dot{\q}}} = {\dot{\q}} / \|{\dot{\q}}\|$ is the normalized velocity.
In other words, the energized system projects accelerations of the geometric policy $\pi$ orthogonal to the direction of motion so they curve the system without changing its speed.

\subsection{Constraints in generalized coordinates}
\label{sec:GeneralizedCoordinates}

We showed above that we can construct fabrics by energizing geometries (see Theorem~\ref{thm:SystemEnergization}) and constraining them (see Theorem~\ref{lma:ConstrainedGeometricFabrics}).
There we calculated explicit constraint forces in the ambient space for analysis, but in many cases, it's natural to express constrained systems in generalized coordinates such as the joint angles of the robot. This section presents a fundamental result showing how to calculate an energized policy in generalized coordinates.

First we formally derive how to constrain systems using generalized coordinates. Let $\mathbf{M}(\mathbf{x},{\dot{\x}}){\ddot{\x}} + \mathbf{f}(\mathbf{x},{\dot{\x}}) = \mathbf{0}$ be any system of this form defined on a space $\mathcal{X}$ ($\mathbf{M}$ need not derive from a Finsler energy and $\mathbf{f}$ captures all system forces). Suppose $\phi:\mathcal{Q}\rightarrow\mathcal{X}$ defines a surface in $\mathcal{X}$ with generalized coordinates $\mathcal{Q}$ and Jacobian $\mathbf{J}=\partial_\mathbf{q}\phi$ (note that the rows of $\mathbf{J}$ here are orthogonal to those of the implicit surface Jacobian due to the differing representation). We have $\mathbf{x} = \phi(\mathbf{q})$, ${\dot{\x}} = \mathbf{J}{\dot{\q}}$, and ${\ddot{\x}} = \mathbf{J}{\ddot{\q}} + {\dot{\J}}{\dot{\q}}$; plugging that last expression for ${\ddot{\x}}$ into the system equation and multiplying through by $\mathbf{J}$ gives an expression for the constrained system in generalized coordinates:
\begin{align} \label{eqn:PullbackDerived}
    \Big(\mathbf{J}^\tr\mathbf{M}\mathbf{J}\Big) {\ddot{\q}} + \mathbf{J}^\tr\big(\mathbf{f} + \mathbf{M}{\dot{\J}}{\dot{\q}}\big) = \mathbf{0}.
\end{align}
This equation has the same form as the original system in $\mathcal{X}$, but with $\mathbf{M}\rightarrow \mathbf{J}^\tr\mathbf{M}\mathbf{J}$ and $\mathbf{f}\rightarrow \mathbf{J}^\tr\big(\mathbf{f} + \mathbf{M}{\dot{\J}}{\dot{\q}}\big)$. This transformation is known as the {\em pullback} of the system into the generalized coordinates, and can be expressed as an operation on system parameters
\begin{align} \label{eqn:SpecPullback}
    \mathrm{pull}_\phi\big[(\mathbf{M},\mathbf{f})\big] = \Big(\mathbf{J}^\tr\mathbf{M}\mathbf{J},\,\mathbf{J}^\tr\big(\mathbf{f} + \mathbf{M}{\dot{\J}}{\dot{\q}}\big)\Big).
\end{align}
For Lagrangian systems, this pullback operation is consistent with the generalized coordinate calculation of the Euler-Lagrange equation. Specifically, if $\mathcal{L}(\mathbf{x},{\dot{\x}})$ is any Lagrangian in $\mathcal{X}$, we can define
\begin{align}
    \mathrm{pull}_\phi\big[\mathcal{L}\big] = \mathcal{L}\Big(\phi(\mathbf{q}), \frac{d}{dt}\phi(\mathbf{q})\Big).
\end{align}
Then if $\mathrm{EL}[\mathcal{L}] = \big(\mathbf{M},\bm{\xi}\big)$ denotes the application of the Euler-Lagrange equation to construct equations of motion $\mathbf{M}{\ddot{\x}} + \bm{\xi} = \mathbf{0}$, we have
\begin{align} \label{eqn:PullbackCommutesWithEL}
    \mathrm{pull}_\phi\big[\mathrm{EL}[\mathcal{L}]\big]
    = \mathrm{EL}\big[\mathrm{pull}_\phi[\mathcal{L}]\big]
    = \mathrm{pull}_\phi\big[(\mathbf{M},\bm{\xi})\big].
\end{align}

The following theorem reveals that a constrained energized system in generalized coordinates can be constructed by first constraining the metric-weighted geometry then energizing directly in the generalized coordinates.

\begin{theorem}\label{thm:FabricsInGeneralizedCoords}
Let $\mathcal{L}_e^\mathcal{X}(\mathbf{x},{\dot{\x}})$ be a Finsler energy on $\mathcal{X}$ with equations of motion $\mathbf{M}{\ddot{\x}} + \bm{\xi} = \mathbf{0}$ and let $\pi(\mathbf{x},{\dot{\x}})$ be an HD2 geometry on $\mathcal{X}$. If $\phi:\mathcal{Q}\rightarrow\mathcal{X}$ is a full rank differentiable map defining a constraint in $\mathcal{X}$ with generalized coordinates $\mathcal{Q}$, then
\begin{align}\nonumber
    \mathrm{pull}_\phi\big[\mathrm{energize}_{\mathcal{L}_e^\mathcal{X}}[\pi]\big]
    = \mathrm{energize}_{\mathcal{L}_e^\mathcal{Q}}
        \Big[\mathrm{pull}_\phi\big[\big(\mathbf{M}, -\mathbf{M}\pi\big)\big]\Big],
\end{align}
where $\mathcal{L}_e^\mathcal{Q}=\mathrm{pull}_\phi[\mathcal{L}_e^\mathcal{X}]$ and $\big(\mathbf{M}, -\mathbf{M}\pi\big)$ denotes the parameters of the metric weighted geometry $\mathbf{M}\big({\ddot{\x}} - \pi(\mathbf{x},{\dot{\x}})\big) = \mathbf{0}$.
\end{theorem}
\begin{proof}
See \cite{Wyk2021GeometriFabricsPhysicsArXivAbbrev} Appendix~\ref{apx:TheoremFabricsInGeneralizedCoords}.
\end{proof}

\section{Designing with fabrics}
\label{sec:DesigningFabrics}

Theorem~\ref{thm:FabricsInGeneralizedCoords} is fundamental to fabric design because it enables us to design using transform trees. A transform tree is a tree of differentiable maps, where each node is a task space and each edge is a differentiable map mapping a parent space (domain) to a child space (co-domain). The tree-structure expresses sparsity patterns important for computational efficiency, but conceptually we can think of the tree as a single differentiable map linking the root space to the stacked joint output space of all nodes. From that view it's clear that we can exploit Theorem~\ref{thm:FabricsInGeneralizedCoords} to say that the resulting constrained system constraining the joint ambient space fabric is equivalent to the system defined by pulling the ambient space fabric to the root.
In light that, it's often useful to design geometric fabrics in parts, viewing each part as an HD2 geometric policy $\pi$ with an associated priority metric defined by a Finsler energy $\mathcal{L}_e$, each living on a node space of a transform tree. Then computationally we can pull the pieces back to the root exploiting sparsity in the same way as in earlier work \cite{cheng2018rmpflow}.
When designing we can reason about what the geometry of each part should be and how we should prioritize it,
understanding that they will combine a metric-weighted average of geometries.

Let $\mathcal{G} = \{\mathcal{L}, \mathbf{f}\}$ compactly denote a system of the form $\mathbf{M}{\ddot{\x}} + \mathbf{f} = \mathbf{0}$ paired with a Finsler energy from which $\mathbf{M}$ derives. We call $\mathcal{G}$ a {\em fabric component}. Note that this system by itself is not a geometric fabric, which must be a bent Finsler system. Here $\mathbf{f}$ can be arbitrary. For instance, it often denotes $\mathbf{f} = -\mathbf{M}\pi$ derived from geometric policy $\pi$.

The scaled sum of two fabric components $\mathcal{G}_1 = \{\mathcal{L}_1,\mathbf{f}_1\}$ and $\mathcal{G}_2 = \{\mathcal{L}_2,\mathbf{f}_2\}$ is given by
\begin{align}
    \alpha \mathcal{G}_1 + \beta \mathcal{G}_2
    = \big\{\alpha \mathcal{L}_1 + \beta \mathcal{L}_2, \alpha \mathbf{f}_1 + \beta \mathbf{f}_2\big\},
\end{align}
where $\alpha,\beta\in\Rbb_+$. Likewise, the pullback of a fabric component across a map $\mathbf{x} = \phi(\mathbf{q})$ with Jacobian $\mathbf{J}$ is
\begin{align}
    \mathrm{pull}_\phi\big[\{\mathcal{L}, \mathbf{f}\}\big]
    = \Big\{\mathrm{pull}_\phi[\mathcal{L}], \mathbf{J}^T\big(\mathbf{f} + \mathbf{M}{\dot{\J}}{\dot{\q}}\big)\Big\}.
\end{align}
\iffalse
See \cite{Wyk2021GeometriFabricsPhysicsArXivAbbrev} Appendix~\ref{apx:DerivationsOfFabricComponentAlgebra} for formal derivations of these formulas.
\fi
These expressions are easily derived using the basic definitions of Finsler energy pullback and fabric pullback from Section~\ref{sec:GeneralizedCoordinates}. It is easy to show that summation is linear (which implies associativity and commutativity), and pullback across the composition of functions is equal to the pullback across each $\mathrm{pull}_{\phi_2\circ\phi_1}[\mathcal{G}] = \mathrm{pull}_{\phi_1}[\mathrm{pull}_{\phi_2}[\mathcal{G}]]$.

Given a transform tree, we populate the nodes with these fabric components and exploit this algebra to
recursively pull each component from child to parent and sum the results at the parent (with optional scalar weights). What results is a final fabric component at the root $\wt{\mathcal{G}} = \{\wt{\mathcal{L}}, \wt{\mathbf{f}}\}$ and associated Finsler energy $\wt{\mathcal{L}}$. $\wt{\mathcal{G}}$ represents system $\wt{\mathbf{M}}{\ddot{\q}} + \wt{\mathbf{f}} = \mathbf{0}$ in generalized coordinates defining an HD2 geometric policy $\wt{\pi} = -\wt{\mathbf{M}}^{-1}\wt{\mathbf{f}}$ since it is formed as a metric-weighted average of pullback geometries from each component (and the metrics are HD0). The final fabric is simply the energization of this geometry with the associated Finsler energy:
\begin{align}
    \mathcal{F} = \mathrm{energize}_{\wt{\mathcal{L}}}\big[\wt{\pi}\big].
\end{align}
In practice, we also leverage the geometric path consistency of the final geometric policy to regulate speed without affecting the path traced during execution (see \cite{ratliff2021geometry}).

\iffalse
\subsection{Speed control through damping regulation}
\label{sec:SpeedControl}

Because we have the flexibility now to encode the nominal behavior into the fabric itself we can leverage accelerations along the direction of motion (known to leave the geometric paths unchanged) to design dampers that regulate speed without affecting the behavior.
We use dampers of the form
\begin{align}
    \mathbf{f}_\mathrm{damp} = -\beta(\mathbf{q},{\dot{\q}})\mathbf{M}(\mathbf{q},{\dot{\q}}){\dot{\q}},
\end{align}
with $\beta\in\Rbb_+$ so that the resulting acceleration $\mathbf{M}^{-1}\mathbf{f}_\mathrm{damp} = -\beta{\dot{\q}}$ acts solely along the direction of motion.
The potential then has a dual role to (1) speed the system (the acceleration component parallel to the direction of motion) and (2) force the system off a nominal path toward the goal (the acceleration component orthogonal to the direction of motion).
See Appendix~\ref{apx:FormulasSpeedControl} for a discussion on how to calculate $\beta$ to effectively control a target measure of speed.

\subsection{Tools for fabric design}
\label{sec:ToolsForFabricDesign}

Finsler energies and HD2 geometries share a common requirement of being Homogeneous of Degree 2 (HD2). \cite{Wyk2021GeometriFabricsPhysicsArXivAbbrev} Appendix~\ref{apx:DesigningHD2Terms} presents a number of methods for designing such functions which we use throughout our experiments. At a high-level these tend to revolve around designing HD0 terms (which can rely on velocity directionality but not the norm (speed)), and scaling them by some measure of Finsler energy (such as the squared norm of the velocity) to bring them from HD0 to HD2. We use these techniques to design the fabrics used in our experiments as described in Sections \ref{sec:ControlledParticleExperiments} and  \ref{subsec:FabricsDesign}.
\fi

\section{Particle Experiments} \label{sec:ControlledParticleExperiments}
Geometric fabrics have two parts to their generalization: (1) moving from Riemannian geometries to Finsler geometries; (2) bending those Finsler geometries. These experiments explore how both of those pieces affect the system by comparing four types of systems: (i) strictly Riemannian systems (classical mechanical systems); (ii) standard Finsler systems (enabling velocity-dependent metrics, but without bending terms); (iii,iv) bent variants of each of the first two. In the discussion below, we often refer to both Riemannian and Finsler systems as {\em unbent fabrics} and the bent variants as geometric fabrics; there are four combinations in total: unbent Riemannian, unbent Finsler, bent Riemannian, bent Finsler (this last one being a full geometric fabric).

We start with a collection of particles at rest to the right of a circular obstacle and we observe the behavior of different variants controlled to two different speeds attracting the particles to a point to the left of the obstacle (small square) under the influence of obstacle avoidance.

\subsection{Point Attraction}
\label{subsec:point_attraction}
This term uses task map $\mathbf{x} = \phi(\mathbf{q}) = \mathbf{q} - \mathbf{q}_d$ where $\mathbf{q}$, $\mathbf{q}_d \in \mathbb{R}^2$ are the current and desired particle position in Euclidean space. The metric is designed as $\mathbf{G}_\psi(\mathbf{x}) = (\widebar{m} - \underline{m}) e^{-(\alpha_m \|\mathbf{x}\|) ^ 2} \mathbf{I} + \underline{m} \mathbf{I}$, where $\widebar{m}$, $\underline{m} \in \mathbb{R}^+$ are the upper and lower isotropic masses, respectively, and $\alpha_m \in \mathbb{R}^+$ controls the width of the radial basis function.
We design the potential gradient directly as $\partial_\mathbf{q} \psi(\mathbf{x}) = \mathbf{M}_\psi(\mathbf{x}) \partial_\mathbf{q} \psi_1(\mathbf{x})$, with $\psi_1(\mathbf{x}) = k \left( \|\mathbf{x}\| + \frac{1}{\alpha_\psi}\log(1 + e^{-2\alpha_\psi \|\mathbf{x}\|} \right)$, where $k \in \mathbb{R}^+$ controls the overall gradient strength, $\alpha_\psi \in \mathbb{R}^+$ controls the transition rate of $\partial_\mathbf{q} \psi_1(\mathbf{x})$ from a positive constant to 0.
For the unbent fabric, the system derives from $\mathcal{L} = {\dot{\x}}^T \mathbf{G}_\psi(\mathbf{x}) {\dot{\x}} - \psi(\mathbf{x})$. For the geometric fabric, we use fabric component defined by geometry ${\ddot{\x}} = -\partial_\mathbf{x} \psi_1(\mathbf{x})$ and energy $\mathcal{L}_e = {\dot{\x}}^T \mathbf{G}_\psi(\mathbf{x}) {\dot{\x}}$.

\subsection{Circular Object Repulsion}
\label{subsec:circular_object_repulsion}
This term uses task map $x = \phi(\mathbf{q}) = \frac{\|\mathbf{q} - \mathbf{q}_o\|}{r} - 1$, where $\mathbf{q}_o$ is the origin of the circle and $r$ is its radius. Two different metrics are designed to prioritize object repulsion, $G_b (x) = \frac{k_b}{x^2}$ and $G_b (x, \dot{x}) = s(\dot{x}) \frac{k_b}{x^2}$. Moreover, $k_b \in \mathbb{R}^+$ is a barrier gain and $s(\dot{x}) = 1$ if $\dot{x} < 0$ and $s(\dot{x}) = 0$, otherwise.
We design the potential gradient as $\partial_\mathbf{q} \psi_b(x) = M_b (x) \partial_\mathbf{q} \psi_{1,b}(x)$ with $\psi_{1,b}(x) = \frac{\alpha_b}{2 x^8}$, and $\alpha_b \in \mathbb{R}^+$
For the unbent fabric, the system derives from either $\mathcal{L}_1 = G_b(x) \dot{x}^2 - \psi_b(x)$ (Riemannian) or  $\mathcal{L}_2 = G_b(x, \dot{x}) \dot{x}^2 -
\psi_b(x)$ (Finsler). For the geometric fabric, we use fabric component defined by $\ddot{x} = - s(\dot{x}) \dot{x}^2 \partial_x \psi_{1,b}(x)$ and energy either $\mathcal{L}_1$ or $\mathcal{L}_2$.

\subsection{Discussion}
\label{subsec:particle_discussion}

Fig. \ref{fig:gf_vs_lds_particles} visualizes the results of all for variants each under two different speeds $v_d=2,4$ to analyze degree of consistency under speed changes. For each variant, all particles avoid the object and reach the target position except the one shot directly at the object's local minimum. However, we observe several differences in behavior.

First, the unbent fabrics (left column) have more pronounced changes in paths across the different speed levels. This is amplified for Finsler systems (bottom two rows) using $G_b(x, \dot{x})$ since the velocity gate does not modulate the obstacle avoidance policy. Instead, the mass of the obstacle avoidance policy vanishes, while components of its force remain. This effect is amplified when traveling at a higher velocity, producing the ``launching'' artifacts. In contrast, geometric fabrics (right column) produce more consistent paths across speed levels without any launching artifacts.

Second, Finsler metrics (bottom two rows) facilitate straight-line motion to the desired location once past the obstacle by allowing the system to forget about it as a function of directionality. Riemannian metrics are unable to represent that directional dependence (top two rows) and can be seen clinging to the obstacle even once past it.

Finally, the repulsive forces from the unbent fabrics (right column) consistently push the center particle farther out from the object than with geometric fabrics (right column). Importantly, geometric fabrics allow for heightened obstacle avoidance behavior without shifting the system minima.

\begin{figure}[!t]
  \centering
  \iffalse
  \includegraphics[width=.8\linewidth]{figures/lds_vs_gf.png}
  \else
  \includegraphics[width=.6\linewidth]{figures/lds_vs_gf.png}
  \fi
  \vspace{-4mm}
  \caption{Particle behavior for an unbent fabric and a geometric fabrics with different metric designs. See Section \ref{subsec:particle_discussion} for discussion. }
  \vspace{-6mm}
  \label{fig:gf_vs_lds_particles}
\end{figure}
\section{Franka Panda Experiments}
\label{sec:SystemExperiments}

\iftrue
\begin{figure*}[!t]
  \centering
  \includegraphics[height=.15\linewidth]{figures/wall_v2.png}
  \hfill
  \includegraphics[height=.15\linewidth]{figures/ring_collage_v2.png}
  \hfill
  \includegraphics[height=.15\linewidth]{figures/moving_balls.png}
  \hfill
  \vspace{-4mm}
  \caption{Left: Reaching a target blocked by a wall (see Section~\ref{subsec:WallBarrierApproach}). Middle: Ring reaching experiment comparing fabrics (top) and RMPs (bottom) (see Section~\ref{subsec:RingNav}). Right: dynamic obstacle avoidance experimental setup (see Section~\ref{subsec:DynamicReaching}).}
  \vspace{-4mm}
  \label{fig:experiments}
\end{figure*}
\fi

Our second set of experiments is a study of geometric fabrics on a full 7 degree-of-freedom robot arm. The geometric fabric is designed once, and tested across a variety of problems. We compare its performance to our current, best performing Riemmanian Motion Policy (RMP) that has powered our manipulation processes for a variety of settings. Both policies are evaluated at 100 Hz and integrated with first-order Euler routines. Note that since this is purely a motion generation problem, the sim-to-real gap is minimal since smooth trajectories can easily be executed on physical robots using standard trajectory following controllers. We, therefore, focus only a simulated analysis here.

\subsection{Geometric Fabrics Design}
\label{subsec:FabricsDesign}
The geometric fabric policy ultimately consists of only three differently designed components. Each component consists of a Finsler energy $\mathcal{L}_e(\mathbf{x}, {\dot{\x}})$ and policy $\pi(\mathbf{x}, {\dot{\x}})$.

\subsubsection{Joint Space Attraction}
The Finsler energy is designed as $\mathcal{L}_e = \frac{m}{2} \|{\dot{\x}}\|^2$, where $m \in \mathbb{R}^+$ is the amount of isotropic mass. The policy is designed as $\pi = -\|{\dot{\x}}\|^2 \partial_\mathbf{x} \psi(\mathbf{x})$ with $\psi(\mathbf{x}) = k \log \left(\cosh (\alpha (\mathbf{x}_d - \mathbf{x}))\right) $ where $k$, $\alpha \in \mathbb{R}^+$ are gains that control scaling and sharpness, respectively, and $\mathbf{x}$ and $\mathbf{x}_d$ are the current and desired joint space positions. This component acts to posture the arm largely within the null space of higher-priority tasks such as target reaching.

\subsubsection{Distance Space Repulsion}
The Finsler energy is designed as $\mathcal{L}_e = \frac{k}{2 x} s(\dot{x}) \dot{x}^2$, where $k \in \mathbb{R}^+$ scales the barrier function and $s(\dot{x}) = 1$ if $\dot{x} < 0$ and $s(\dot{x}) = 0$, otherwise.
We define a potential
$\psi(x) = \frac{k_b}{x} + \frac{k_r}{\alpha} \log \left(1 + e^{-\alpha (x - x_o)} \right)$, where $k_b$, $k_r \in \mathbb{R}^+$ control the scaling of the barrier and soft rectified linear unit (ReLU) functions, respectively and $x \in \mathbb{R}^+$ is a signed distance space (e.g., distance to object). Moreover, $\alpha \in \mathbb{R}^+$ controls the sharpness in the ReLU and $x_o \in \mathbb{R}^+$ shifts the activation of the ReLU. We deploy this potential alongside a geometric variant $\pi = -\dot{x}^2 \partial_x \psi(x)$, for joint-limit avoidance, self-collision avoidance, and object-collision avoidance. The geometric policies are weighed significantly higher than potentials, which play the role of soft-penalty on the obstacle constraint here. The potential policy is included to shift the system minima away from numerical singularities, e.g., as $x \to 0$, then $\dot{x} \to 0$ and $\partial_x \psi(x) \to \infty$, which generates numerical issues for calculating the geometric policy.

\subsubsection{End-effector Space Attraction}
\label{subsubsec:attraction}
Using $\overline{m}$, $\underline{m} \in \mathbb{R}^+$ to denote maximum and minimum isotropic masses, respectively,
the Finsler energy is designed as $\mathcal{L}_e = \|{\dot{\x}}\|^2\left( \frac{1}{2} (\overline{m} - \underline{m}) (\tanh(-\alpha \|\mathbf{x}_d - \mathbf{x}\|) + 1) + \underline{m} \right)$, where $\alpha \in \mathbb{R}^+$ controls the sharpness of the hyperbolic tangent function, and $\mathbf{x},\mathbf{x}_d \in \mathbb{R}^3$ are the position and target position of the end-effector, respectively. The
potential is defined as
$\psi(\mathbf{x}) = k \log \left(\cosh (\alpha (\mathbf{x}_d - \mathbf{x}))\right) $, where $k,\alpha \in \mathbb{R}^+$ are gains that control scaling and sharpness, respectively. We deploy both the potential and a geometric variant $\pi = -\|{\dot{\x}}\|^2 \partial_\mathbf{x} \psi(\mathbf{x})$. Each are individually tuned, with the geometric variant being generally more active and the potential there to ensure convergence in alignment with our theory. Applying this design to two axes at the end-effector can also realize orientation control as in \cite{ratliff2018rmps}.

\subsection{Wall Barrier Approach} \label{subsec:WallBarrierApproach}

\iffalse
\begin{figure}[!t]
  \centering
  \includegraphics[width=.9\linewidth]{figures/wall_v2.png}
  \vspace{-4mm}
  \caption{Robot reaching for a target (orange sphere) while blocked by an obstacle (blue) under a geometric fabric policy (left) and an RMP (right).}
  \label{fig:wall_approach}
\end{figure}
\fi

To assess how well the policies can reach 3D end-effector position targets in the presence of boundary constraints, we incrementally advance the end-effector target every five seconds towards a wall boundary and monitor convergence
\iftrue
(see Fig. \ref{fig:experiments} (left)).
\else
(see Fig. \ref{fig:wall_approach}).
\fi
Ideally, a policy should accurately converge to the target while respecting the boundary constraint, but in practice boundary potentials can prevent that. In this experiment, we see geometric fabrics converge more accurately in general to these targets than the RMP policy. Both policies exhibit sub-millimeter convergence error when the target is 10 cm or more from the wall, but convergence incrementally degrades as the target approaches the wall due to conflicting potentials. RMPs already start at degrading at $\sim$9 cm while geometric fabrics begin degrading only at $\sim$6 cm since much of the boundary policy is encoded in a geometric policy. Moreover, geometric fabrics on average approach the wall about 50\% closer than the RMP.
\iftrue
Figure~\ref{fig:experiments} (left)
\else
Fig.~\ref{fig:wall_approach}
\fi
depicts this difference in convergence near the wall.

\subsection{Ring Constricted Navigation} \label{subsec:RingNav}

\iffalse
\begin{figure}[!t]
  \centering
  \includegraphics[width=.9\linewidth]{figures/ring_collage_v2.png}
  \vspace{-4mm}
  \caption{Robot executing reach to target (orange spheres) through a series of obstacle rings (blue) of increasingly smaller diameter with a geometric fabric (top row) and RMP (bottom row). Green spheres approximate robot geometry for collision avoidance. }
  \label{fig:ring_navigation}
\end{figure}
\fi

To assess convergence to end-effector targets in the presence of increasingly constricted passages, we create a sequence of rings with incrementally smaller radii that the end-effector must pass through to reach its target. The ring radius starts at 17 cm and decreases to 13 cm. The rings are constructed from small spheres (each with radius 2.5 cm) as shown in
\iftrue
Figure \ref{fig:experiments} (middle).
\else
Fig.~\ref{fig:ring_navigation}.
\fi
Distance spaces are pair-wise between robot spheres and obstacle spheres. As shown, the navigation from a geometric fabric allows the robot to successfully pass through all rings and reach its target. In contrast, the RMP cannot pass through the last two rings of the smallest diameter. Again, this difference arises with primarily encoding avoidance into the fabric's geometry alongside a diminished repulsive potential function.

\subsection{Dynamic Obstacle Reaching} \label{subsec:DynamicReaching}

\iffalse
\begin{figure}[!t]
  \centering
  \includegraphics[width=.8\linewidth]{figures/moving_balls.png}
  \vspace{-4mm}
  \caption{Robot reaching end-effector targets while avoiding dynamic objects.}
  \label{fig:dynamic_navigation}
\end{figure}
\fi

By nature, both the geometric fabric policy and RMP are reactive. To assess the quality of reactivity, we create a reaching task with both static and dynamic obstacles. The goal is for the policies to reach 19 randomly generated end-effector targets (with one target issued every 5 seconds) while avoiding collision with both static and dynamic obstacles
\iftrue
(see Fig. \ref{fig:experiments} (right)).
\else
(see Fig. \ref{fig:dynamic_navigation}).
\fi
The policies are only aware of the object positions while objects achieve speeds of up to $0.39$ $m/s$. With the geometric fabric policy, the robot reached 16 targets whereas the robot with the RMP reached 11 targets. Moreover, the geometric fabric policy yielded a $0.4 \%$ collision rate in contrast to the RMP's $8.18 \%$,  where collision rate is the percentage of time steps in collision. Since the geometric fabric obstacle avoidance component is predominantly geometric in nature, it can simultaneously be highly influential while not preventing optimization of the end-effector potential, resulting in higher target acquisition rates. The forcing component of the geometric fabric policy aggressively engages only when the robot is close to collision. Altogether, this allows strong obstacle avoidance behavior both close to and far from objects in a way that minimally impedes target acquisition. In contrast, the RMP had to strike a balance between forcing repulsion without impairing the ability to acquire end-effector targets too much. This compromise resulted in diminished performance.

\section{Conclusion}
Geometric fabrics naturally generalize classical mechanical systems to form a new physics of behavior. These systems provide important design expressivity while maintaining the theoretical properties that make classical systems powerful in control. Many existing techniques in the literature can naturally leverage fabrics (e.g., reinforcement learning as in \cite{LiRmp2RobotLearning2021}, imitation learning, planning), and it is a point of future work to explore many of these areas. Our theoretical and experimental analyses already demonstrate their utility but we believe the bulk of their application has yet to come.
\iftrue
For extended discussions, proofs, and additional exhibitions of geometry fabrics (including full pose control and equality constraint handling), see the extended version \cite{Wyk2021GeometriFabricsPhysicsArXivAbbrev}, video, and website (\url{https://sites.google.com/nvidia.com/geometric-fabrics-behavior}).
\else
\fi
\section{Preliminaries}
\input{sections/appendix_manifolds}
\input{sections/appendix_calculus}

\section{Additional comments on related work}

\subsection{Generalization of Dynamic Movement Primitive (DMP) principles}
\label{sec:GeneralizingDmps}

In the Related Work section (see Section~\ref{sec:RelatedWork}), we discussed how Dynamic Movement Primitives (DMPs) are perturbations from a stable second-order differential equation and, in principle, can easily be generalized to leverage stable differential equations beyond simple linear equations. We describe that connection in more detail here. Studying these generalized perturbed systems is an area of future work, but we note here that geometric fabrics are not inherently competitive with the principles of DMPs, but rather complimentary.

A DMP is a nonlinear force perturbation away from an intrinsically stable linear second-order differential equation \cite{IjspeertDMPs2013}. When the force perturbation vanishes over time the system dynamics is asymptotically governed by the linear differential equation and therefore becomes a point attractor; likewise, when the nonlinear force is periodic the system converges to a stable limit cycle. In all cases, the key innovation is perturbing a stable differential equation using either a vanishing or cyclic perturbation so that stability of the combined system is governed by the stability of the underlying differential equation.

These ideas can easily be generalized to geometric mechanics noting that inherently stable classical mechanical systems can be used as a drop-in replacement for the linear second-order differential equation to inherit stability. Since geometric fabrics are a generalization of those classical mechanical systems, the same applies for fabrics. Geometric fabrics can act as a stable underlying nonlinear second-order differential equation that we can perturb away from using a trained perturbation component generalizing the techniques studies in DMPs while gaining more expressive underlying system behaviors from this broader class of geometric fabric equations.

\subsection{Geometric Dynamical Systems vs Finsler Systems}
\label{apx:Gds}

Geometric Dynamical Systems (GDSs) were introduced in \cite{cheng2018rmpflow} as a stable class of RMPs for stable system design. But we now understand GDSs are more elegantly and concretely expressed as Finsler systems characterized in Subsection~\ref{sec:FinslerAndHd2Geometries}. Moreover, such systems are inherently covariant since they derive from the Euler-Lagrange equations. This section presents that connection in detail.

We start by reviewing Geometric Dynamical Systems (GDS) and express them in a (tensor notion) form that makes it easier to see the relationship to Lagrangian and Finsler systems. Let $\mathbf{G}(\mathbf{q},{\dot{\q}})$ be a velocity dependent metric tensor and define an energy function as $\mathcal{K}(\mathbf{q},{\dot{\q}}) = \frac{1}{2}{\dot{\q}}^\tr\mathbf{G}(\mathbf{q},{\dot{\q}}){\dot{\q}}$. Following the notation of \cite{cheng2018rmpflow}, we denote the columns of $\mathbf{G}$ by $\mathbf{g}_i$ so that $\mathbf{G} = [\mathbf{g}_1,\cdots,\mathbf{g}_d]$. An unforced GDS is defined as
\begin{align*}
    \mathbf{M}(\mathbf{q},{\dot{\q}}) {\ddot{\q}} + \bm{\xi}(\mathbf{q},{\dot{\q}}) = \mathbf{0}.
\end{align*}
where $\mathbf{M} = \mathbf{G} + \bm{\Xi}$ has additional mass defined by $\bm{\Xi} = \frac{1}{2}\sum_i\dot{q}^i\partial_{\dot{\q}}\mathbf{g}_i$ (here $\dot{q}^i$ denotes the $i$th component of ${\dot{\q}}$) and $\bm{\xi} = \mathring{\mathbf{G}}{\dot{\q}} - \frac{1}{2}\partial_\mathbf{q}\big({\dot{\q}}^\tr\mathbf{G}{\dot{\q}}\big)$ using $\mathring{\mathbf{G}} = [\mathbf{h}_1,\cdots,\mathbf{h}_d]$ with $\mathbf{h}_i = \big(\partial_\mathbf{q}\mathbf{g}_i\big){\dot{\q}}$. When $\mathbf{G}$ is independent of ${\dot{\q}}$, these equations of motion reduce to the classical equations since $\mathbf{M}\rightarrow\mathbf{G}$ and $\mathring{\mathbf{G}} \rightarrow \dot{\mathbf{G}}$.

GDS systems conserve total energy $\mathcal{E} = \mathcal{K} + \psi(\mathbf{q})$ when forced by a potential function $\psi(\mathbf{q})$ and are shown to optimize $\psi$ when energy is bled off using a damping term $-\mathbf{B}(\mathbf{q},{\dot{\q}}){\dot{\q}}$ for positive definite $\mathbf{B}$ similar to the setting outlined in Theorem~\ref{thm:GeometricFabricStability}.

Using tensor notation with the Einstein summation convention, this system can be expressed $M_{kl}\ddot{q}^l + \xi_{kij}\dot{q}^i\dot{q}^j = 0$ where $M_{kl} = G_{kl} + \Xi_{kl}$ with
\begin{align}
    \label{eqn:GdsEinstein1}
    \Xi_{kl} &= \frac{1}{2}\frac{\partial G_{ik}}{\partial \dot{q}^l}\dot{q}^i \\
    \label{eqn:GdsEinstein2}
    \xi_{kij} &= \frac{1}{2} \left(
                \frac{\partial G_{kj}}{\partial q^i}
                + \frac{\partial G_{ik}}{\partial q^j}
                - \frac{\partial G_{ij}}{\partial q^k}
            \right),
\end{align}
where $\xi_{kij}$ are the standard Christoffel symbols of the first kind \cite{LeeRiemannianManifolds97}. The only augmentation to the classical system is the term $\Xi$ which vanishes when $\mathbf{G}$ is independent of velocity.

A GDS by itself is not covariant to change of coordinates \cite{BylardBonalliEtAl2021}, but can be made covariant with the use of transform trees (see {\em structured GDS} in \cite{cheng2018rmpflow} and the discussion above in Appendix~\ref{apx:NoteOnCovariance}).
However, the next theorem shows that Lagrangian systems, and specifically Finsler systems derived from energies of the form $\mathcal{L}_e(\mathbf{q},{\dot{\q}}) = \frac{1}{2}{\dot{\q}}^\tr\mathbf{G}(\mathbf{q},{\dot{\q}}){\dot{\q}}$, make a qualitatively similar, but inherently covariant, alternative to GDS systems.

\begin{theorem} \label{thm:LagrangianEquivalent}
    Let $\mathcal{L} = \frac{1}{2}G_{ij}\dot{q}^i\dot{q}^j$
    be a Finsler energy with $G_{ij}(q,\dot{q})$ a function of both position
    and velocity. Then the Finsler system's equations of motion (defined by the Euler-Lagrange equation) can be expressed
    \begin{align}
        M_{kl}\ddot{q}^l + \Upsilon_{kij} \dot{q}^i\dot{q}^j
        = 0,
    \end{align}
    where $M_{kl} = G_{kl} + \Xi_{kl}$ with velocity curvature mass
    \begin{align}
        \Xi_{kl} =
            \frac{\partial G_{lj}}{\partial \dot{q}^k}\dot{q}^j
            + \frac{\partial G_{kj}}{\partial \dot{q}^l}\dot{q}^j
            + \frac{1}{2}\frac{\partial^2 G_{ij}}{\partial \dot{q}^l\partial \dot{q}^k}\dot{q}^i\dot{q}^j,
    \end{align}
    and velocity-augmented Christoffel symbols (of the first kind) given by
    \begin{align}
        \Upsilon_{kij} =
            \frac{1}{2} \left(
                \frac{\partial G_{kj}}{\partial q^i}
                + \frac{\partial G_{ik}}{\partial q^j}
                - \frac{\partial G_{ij}}{\partial q^k}
                + \frac{\partial^2 G_{ij}}{\partial q^l\partial \dot{q}^k}\dot{q}^l
            \right).
    \end{align}
\end{theorem}
\begin{proof}
    This proof amounts primarily to calculation and gathering together the
    terms appropriately to construct $M_{kl}$ and $\Upsilon_{kij}$.

    The Lagrangian is $\mathcal{L} = \frac{1}{2}G_{ij}\dot{q}^i\dot{q}^j$
    and the Euler-Lagrange equation defines
    $\frac{d}{dt}\frac{\partial \mathcal{L}}{\partial \dot{q}^k}
    - \frac{\partial\mathcal{L}}{\partial q^k} = 0$. That gives
    \begin{align*}
        \underbrace{\frac{d}{dt}\left[\frac{\partial}{\partial \dot{q}^k}\left(
            \frac{1}{2} G_{ij} \dot{q}^i\dot{q}^j
        \right)\right]}_{\mathrm{term 1}}
        - \frac{1}{2}\frac{\partial G_{ij}}{\partial q^k}\dot{q}^j\dot{q}^i
        = 0.
    \end{align*}
    Term 1 becomes
    \begin{align*}
    &\frac{d}{dt}\left[\frac{\partial}{\partial \dot{q}^k}\left(
        \frac{1}{2} G_{ij} \dot{q}^i\dot{q}^j
    \right)\right]
    =
    \frac{d}{dt}\left[
        \frac{1}{2}\frac{\partial G_{ij}}{\partial \dot{q}^k} \dot{q}^i\dot{q}^j
        + G_{kj}\dot{q}^j
    \right] \\
    &\ \ \ = \frac{1}{2}\left(
                \frac{\partial^2G_{ij}}{\partial q^l\partial \dot{q}^k}\dot{q}^l
                + \frac{\partial^2G_{ij}}{\partial \dot{q}^l\partial \dot{q}^k}\ddot{q}^l
            \right)\dot{q}^i\dot{q}^j
            + \frac{1}{2}\frac{\partial G_{ij}}{\partial \dot{q}^k} \ddot{q}^i\dot{q}^j \\
        &\ \ \ \ \ + \frac{1}{2}\frac{\partial G_{ij}}{\partial \dot{q}^k} \dot{q}^i\ddot{q}^j
        + \left(\frac{\partial G_{kj}}{\partial q^l}\dot{q}^l + \frac{\partial G_{kj}}{\partial \dot{q}^l}\ddot{q}^l\right)\dot{q}^j
        + G_{kj} \ddot{q}^j.
    \end{align*}
    Putting this together with the other terms and grouping them as needed,
    we get
    \begin{align*}
        &\frac{d}{dt}\frac{\partial \mathcal{L}}{\partial \dot{q}^k}
            - \frac{\partial\mathcal{L}}{\partial q^k}
        = \Bigg(
            G_{kl}
        + \frac{1}{2}\left[
                \frac{\partial G_{lj}}{\partial \dot{q}^k} \dot{q}^j
                + \frac{\partial G_{il}}{\partial \dot{q}^k} \dot{q}^i
            \right]
            + \frac{\partial G_{kj}}{\partial \dot{q}^l} \dot{q}^j \\
        &\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \,
        + \frac{1}{2}\frac{\partial^2 G_{ij}}{\partial \dot{q}^l\dot{q}^k}
                \dot{q}^i\dot{q}^j
            \Bigg) \ddot{q}^l \\
        &\ \ \ \ \ \ + \left(
            \frac{1}{2}\frac{\partial^2G_{ij}}{\partial q^l\partial \dot{q}^k}\dot{q}^l
            + \frac{\partial G_{kj}}{\partial q^i}
            - \frac{1}{2} \frac{\partial G_{ij}}{\partial q^k}
        \right) \dot{q}^i\dot{q}^j.
    \end{align*}
    Since $G_{ij}$ is symmetric, we have
    $\frac{\partial G_{lj}}{\partial \dot{q}^k} \dot{q}^j
    = \frac{\partial G_{il}}{\partial \dot{q}^k} \dot{q}^i$, and
    \begin{align*}
        \frac{\partial G_{kj}}{\partial q^i}\dot{q}^i\dot{q}^j
        &=
        \frac{1}{2}\frac{\partial G_{kj}}{\partial q^i}\dot{q}^i\dot{q}^j
        + \frac{1}{2}\frac{\partial G_{jk}}{\partial q^i}\dot{q}^i\dot{q}^j \\
        &= \frac{1}{2}\left(
            \frac{\partial G_{kj}}{\partial q^i}
            + \frac{\partial G_{ik}}{\partial q^j}
        \right)\dot{q}^i\dot{q}^j,
    \end{align*}
    so we can rewrite the above expression in its final form
    \begin{align*}
        &\frac{d}{dt}\frac{\partial \mathcal{L}}{\partial \dot{q}^k}
            - \frac{\partial\mathcal{L}}{\partial q^k} \\
        &\ \ \ = \left(
            G_{kl}
            + \frac{\partial G_{lj}}{\partial \dot{q}^k} \dot{q}^j
            + \frac{\partial G_{kj}}{\partial \dot{q}^l} \dot{q}^j
            + \frac{1}{2}\frac{\partial^2 G_{ij}}{\partial \dot{q}^l\dot{q}^k}
                \dot{q}^i\dot{q}^j
            \right) \ddot{q}^l \\
        &\ \ \ \ \ \ + \frac{1}{2}\left(
            \frac{\partial^2G_{ij}}{\partial q^l\partial \dot{q}^k}\dot{q}^l
            + \frac{\partial G_{kj}}{\partial q^i}
            + \frac{\partial G_{ik}}{\partial q^j}
            - \frac{\partial G_{ij}}{\partial q^k}
        \right) \dot{q}^i\dot{q}^j \\
        &\ \ = \Big(G_{kl} + \Xi_{kl}\Big)\ddot{q}^l
            + \Upsilon_{kij}\dot{q}^i\dot{q}^j,
    \end{align*}
    with $\Xi_{kl}$ and $\Upsilon_{kij}$ as given in the above theorem.
\end{proof}
Comparing the terms of Theorem~\ref{thm:LagrangianEquivalent} to Equations~\ref{eqn:GdsEinstein1} and \ref{eqn:GdsEinstein2} we see that Finsler systems of this form contains additional terms deriving from the velocity dependence in both the added mass and the augmented Christoffel symbols. Since these equations are derived using the Euler-Lagrange equation, they're known to be covariant to reparameterization; the added terms are critical for that covariance property to hold. Moreover, the class of all Finsler systems is a superset of these GDS-like Finsler systems, so Finsler systems broadly have many advantages over GDSs.

Both GDSs and Finsler systems, however, generalize metrics to be velocity dependent but remain subject to the other limitations outlined in Subsection~\ref{sec:ShortcomingsOfClassicalMechanics}. Namely, their metrics and associated policies remain coupled. In practice, in the earlier results of \cite{cheng2018rmpflow}, high gains on forcing potentials and damping terms were used in order to leverage the systems for metric design while also reject the resulting irrelevant geometric terms that resulted. Geometric fabrics added a new dimension of flexibility by introducing bending terms to circumvent these issues.

\section{Additional Lemmas}

This Lemma gives the time derivative of a Lagrangian's energy (the Hamiltonian), which is a common calculation needed for multiple proofs below.

\begin{lemma}\label{lma:EnergyTimeDerivative}
Let $\mathcal{L}(\mathbf{q},{\dot{\q}})$ be any Lagrangian with Hamiltonian (energy) $\mathcal{H}_\mathcal{L} = \partial_{\dot{\q}}\mathcal{L}^\tr {\dot{\q}} - \mathcal{L}$. The energy time derivative is
\begin{align} \label{eqn:EnergyTimeDerivative}
    \dot{\mathcal{H}}_\mathcal{L} = {\dot{\q}}^\tr\big(\mathbf{M}{\ddot{\q}} + \bm{\xi}\big),
\end{align}
where $\mathbf{M}$ and $\bm{\xi}$ come from the Lagrangian's equations of motion $\mathbf{M}(\mathbf{q},{\dot{\q}}){\ddot{\q}} + \bm{\xi}(\mathbf{q},{\dot{\q}}) = \mathbf{0}$.
\end{lemma}
\begin{proof}
The calculation is a straightforward time derivative of the Hamiltonian:
\begin{align*}
    \dot{\mathcal{H}}_{\mathcal{L}}
    &= \frac{d}{dt}\Big[\partial_{\dot{\q}}\mathcal{L}^\tr{\dot{\q}} - \mathcal{L}\Big] \\
    &= \big(\partial^2_{{\dot{\q}}\qd}\mathcal{L}\,{\ddot{\q}} + \partial_{{\dot{\q}}\mathbf{q}}\mathcal{L}\,{\dot{\q}}\big)^\tr {\dot{\q}} + \partial_{\dot{\q}}\mathcal{L}^\tr{\ddot{\q}} \\
    &\ \ \ \ \ \ \ \ - \big(\partial_{\dot{\q}}\mathcal{L}^\tr{\ddot{\q}} + \partial_\mathbf{q}\mathcal{L}^\tr{\dot{\q}} \big) \\
    &= {\dot{\q}}^\tr\Big(\partial^2_{{\dot{\q}}\qd}\mathcal{L}\,{\ddot{\q}} + \partial_{{\dot{\q}}\mathbf{q}}\mathcal{L}\,{\dot{\q}} - \partial_\mathbf{q}\mathcal{L}\Big) \\
    &= {\dot{\q}}^\tr\big(\mathbf{M}{\ddot{\q}} + \bm{\xi}\big).
\end{align*}
where $\mathbf{M} = \partial^2_{{\dot{\q}}\qd}\mathcal{L}$ and $\bm{\xi} = \partial_{{\dot{\q}}\mathbf{q}}\mathcal{L}\,{\dot{\q}} - \partial_\mathbf{q}\mathcal{L}$.
\end{proof}

The following lemma collects some results around common matrices and operators that arise when analyzing energy conservation.
\begin{lemma}\label{lma:EnergyProjection}
Let $\mathbf{M}$ be a symmetric positive definite matrix and define $\mathbf{p} = \mathbf{M}{\dot{\x}}$. Then
\begin{align} \label{eqn:Rpe}
    \mathbf{R}_{\mathbf{p}} = \mathbf{M}^{-1} - \frac{{\dot{\x}}\,{\dot{\x}}^\tr}{{\dot{\x}}^\tr\mathbf{M}{\dot{\x}}}
\end{align}
has null space spanned by $\mathbf{p}$ and
\begin{align}
    \mathbf{R}_{\dot{\x}} = \mathbf{M} - \frac{\mathbf{p}\,\mathbf{p}^\tr}{\mathbf{p}^\tr\mathbf{M}^{-1}\mathbf{p}}
\end{align}
has null space spanned by ${\dot{\x}}$. These matrices are related by
$\mathbf{R}_{\dot{\x}} = \mathbf{M}\mathbf{R}_{\mathbf{p}}\mathbf{M}$ and the matrix $\mathbf{M}^{-1}\mathbf{R}_{\dot{\x}} = \mathbf{R}_{\mathbf{p}}\mathbf{M} = \mathbf{P}_e$ is a projection operator of the form
\begin{align} \label{eqn:EnergyProjector}
    \mathbf{P}_e = \mathbf{M}^{\,\frac{1}{2}}\Big[\mathbf{I} - \hat{\mathbf{v}} \hat{\mathbf{v}}^\tr\Big]\mathbf{M}^{-\frac{1}{2}}
\end{align}
where $\mathbf{v} = \mathbf{M}^{\,\frac{1}{2}}{\dot{\x}}$ and $\hat{\mathbf{v}} = \frac{\mathbf{v}}{\|\mathbf{v}\|}$ is the normalized vector.
Moreover, ${\dot{\x}}^\tr\mathbf{P}_e = \mathbf{0}$.
\end{lemma}
\begin{proof}
Right multiplication of $\mathbf{R}_{\mathbf{p}}$ by $\mathbf{p}$ gives:
\begin{align*}
    \mathbf{R}_{\mathbf{p}} \mathbf{p}
    &= \left(\mathbf{M}^{-1} - \frac{{\dot{\x}}\,{\dot{\x}}^\tr}{{\dot{\x}}^\tr\mathbf{M}{\dot{\x}}}\right) \mathbf{M}{\dot{\x}} \\
    &= {\dot{\x}} - {\dot{\x}} \left(\frac{{\dot{\x}}^\tr\mathbf{M}{\dot{\x}}}{{\dot{\x}}^\tr\mathbf{M}{\dot{\x}}}\right) = \mathbf{0},
\end{align*}
so $\mathbf{p}$ lies in the null space. Moreover, the null space is no larger since each matrix is formed by subtracting off a rank 1 term from a full rank matrix.

The relation between $\mathbf{R}_{\mathbf{p}}$ and $\mathbf{R}_{\dot{\x}}$ can be shown algebraically
\begin{align*}
    \mathbf{M}\mathbf{R}_{\mathbf{p}}\mathbf{M}
    &= \mathbf{M}\left(\mathbf{M}^{-1} - \frac{{\dot{\x}}\,{\dot{\x}}^\tr}{{\dot{\x}}^\tr\mathbf{M}{\dot{\x}}}\right)\mathbf{M} \\
    &= \mathbf{M} - \frac{\mathbf{M}{\dot{\x}}\,{\dot{\x}}^\tr\mathbf{M}}{{\dot{\x}}^\tr\mathbf{M}{\dot{\x}}} \\
    &= \mathbf{M} - \frac{\mathbf{p}\,\mathbf{p}^\tr}{\mathbf{p}^\tr\mathbf{M}^{-1}\mathbf{p}} = \mathbf{R}_{\dot{\x}}.
\end{align*}
Since $\mathbf{M}$ has full rank, $\mathbf{R}_{\dot{\x}}$ has the same rank as $\mathbf{R}_{\mathbf{p}}$ and its null space must be spanned by ${\dot{\x}}$ since $\mathbf{R}_{\dot{\x}}\xd = \mathbf{M}\mathbf{R}_{\mathbf{p}}\mathbf{M}{\dot{\x}} = \mathbf{M}\mathbf{R}_{\mathbf{p}}\mathbf{p} = \mathbf{0}$.

With a slight algebraic manipulation, we get
\begin{align*}
    \mathbf{M}\mathbf{R}_{\mathbf{p}}
    &= \mathbf{M}\left(\mathbf{M}^{-1} - \frac{{\dot{\x}}\,{\dot{\x}}^\tr}{{\dot{\x}}^\tr\mathbf{M}{\dot{\x}}}\right) \\
    &= \mathbf{M}^{\,\frac{1}{2}} \left(\mathbf{I} - \frac{\mathbf{M}^{\,\frac{1}{2}}{\dot{\x}}\,{\dot{\x}}^\tr\mathbf{M}^{\,\frac{1}{2}}}{{\dot{\x}}^\tr\mathbf{M}^{\,\frac{1}{2}}\mathbf{M}^{\,\frac{1}{2}}{\dot{\x}}}\right) \mathbf{M}^{-\frac{1}{2}} \\
    &= \mathbf{M}^{\,\frac{1}{2}} \left(\mathbf{I} - \frac{\mathbf{v}\vv^\tr}{\mathbf{v}^\tr\mathbf{v}}\right) \mathbf{M}^{-\frac{1}{2}} \\
    &= \mathbf{P}_e
\end{align*}
since $\frac{\mathbf{v}\,\mathbf{v}^\tr}{\mathbf{v}^\tr\mathbf{v}} = \hat{\mathbf{v}}\hat{\mathbf{v}}^\tr$. Moreover,
\begin{align*}
    \mathbf{P}_e\mathbf{P}_e
    &=
    \mathbf{M}^{\,\frac{1}{2}} \left(\mathbf{I} - \hat{\mathbf{v}}\hat{\mathbf{v}}^\tr\right) \mathbf{M}^{-\frac{1}{2}}\\
    &\ \ \ \ \ \ \cdot
    \mathbf{M}^{\,\frac{1}{2}} \left(\mathbf{I} - \hat{\mathbf{v}}\hat{\mathbf{v}}^\tr\right) \mathbf{M}^{-\frac{1}{2}} \\
    &= \mathbf{M}^{\,\frac{1}{2}} \left(\mathbf{I} - \hat{\mathbf{v}}\hat{\mathbf{v}}^\tr\right) \left(\mathbf{I} - \hat{\mathbf{v}}\hat{\mathbf{v}}^\tr\right) \mathbf{M}^{-\frac{1}{2}} \\
    &= \mathbf{M}^{\,\frac{1}{2}} \left(\mathbf{I} - \hat{\mathbf{v}}\hat{\mathbf{v}}^\tr\right) \mathbf{M}^{-\frac{1}{2}} = \mathbf{P}_e,
\end{align*}
since $\mathbf{P}_\perp = \mathbf{I} - \hat{\mathbf{v}}\hat{\mathbf{v}}^\tr$ is an orthogonal projection operator. Therefore, $\mathbf{P}_e^2 = \mathbf{P}_e$ showing that it is a projection operator as well.

Finally, we have
\begin{align*}
    {\dot{\x}}^\tr\mathbf{P}_e
    &= {\dot{\x}}^\tr\mathbf{M}\mathbf{R}_{\mathbf{p}}
    = {\dot{\x}}^\tr\mathbf{M}\left(\mathbf{M}^{-1} - \frac{{\dot{\x}}\,{\dot{\x}}^\tr}{{\dot{\x}}^\tr\mathbf{M}{\dot{\x}}}\right)\\
    &= \left[{\dot{\x}}^\tr - \left(\frac{{\dot{\x}}^\tr\mathbf{M}{\dot{\x}}}{{\dot{\x}}^\tr\mathbf{M}{\dot{\x}}}\right){\dot{\x}}^\tr\right]
    = \mathbf{0}.
\end{align*}
\end{proof}

\section{Proofs of lemmas and theorems in the paper}
\label{apx:proofs}

\subsection{Lemma~\ref{lma:EnergyConservation} on energy conservation and dissipation}
\label{apx:LemmaEnergyConservation}
Proof of Lemma~\ref{lma:EnergyConservation} on energy conservation and dissipation of Finsler systems with zero-work modification terms. See Subsection~\ref{sec:BendingFinslerGeometries} for the full statement.
\begin{proof}
This energy is conserved if its time derivative is zero.
Substituting ${\ddot{\q}} = -\mathbf{M}^{-1}\big(\bm{\xi} + \mathbf{f}_f + \partial\psi + \mathbf{B}{\dot{\q}}\big)$ into Equation~\ref{eqn:EnergyTimeDerivative} of Appendix Lemma~\ref{lma:EnergyTimeDerivative} using Lagrangian $\mathcal{L}(\mathbf{q}, {\dot{\q}}) = \mathcal{L}_e - \psi$
gives
\begin{align*}
    \dot{\mathcal{H}}_{\mathcal{L}_e}
    &= {\dot{\q}}^\tr\Big(\mathbf{M} \big(-\mathbf{M}^{-1}(\bm{\xi} + \mathbf{f}_f + \partial\psi + \mathbf{B}{\dot{\q}})\big)+ \bm{\xi} + \partial\psi\Big) \\
    &= {\dot{\q}}^\tr\big(-(\bm{\xi}+\partial\psi)
        +(\bm{\xi} + \partial\psi)
        - (\mathbf{f}_f + \mathbf{B}{\dot{\q}})\big) \\
    &= -{\dot{\q}}^\tr\mathbf{f}_f - {\dot{\q}}^\tr\mathbf{B}{\dot{\q}}.
\end{align*}
Therefore, energy reduces at a rate $-{\dot{\q}}^\tr\mathbf{B}{\dot{\q}}$ (and is conserved when $\mathbf{B}=\mathbf{0}$) if and only if $\mathbf{f}_f$ is zero-work with ${\dot{\q}}^\tr\mathbf{f}_f = 0$ everywhere.

Finally, when $\mathbf{f}_f$ is geometric $\mathbf{M}{\ddot{\q}} + \bm{\xi} + \mathbf{f}_f = \mathbf{0}$ is geometric since $\mathbf{M}$ is HD0 and $\bm{\xi}$ is HD2. Moreover, since $\mathcal{L}_e$ is lower bounded by 0, $\mathcal{H}=\mathcal{L}_e+\psi$ is lower bounded when $\psi$ is.
\end{proof}

\subsection{Lemma~\ref{lma:ConstrainedGeometricFabrics} on constrained fabrics}
\label{apx:LemmaConstrainedGeometricFabrics}

Proof of Lemma~\ref{lma:ConstrainedGeometricFabrics} in  Subsection~\ref{sec:GeometricFabrics} showing that constrained fabrics are also fabrics.
\begin{proof}
The formula derives from the standard algebra of constrained systems used in classical mechanics \cite{ClassicalMechanicsTaylor05}, which we derive here for completeness. The Principle of Virtual Work tells us the constraint force takes the form $\mathbf{f}_c = \mathbf{J}^\tr\bm{\lambda}$ for unknown Lagrange multipliers $\bm{\lambda}$. (Since $\mathbf{J}$'s rows are orthogonal to the constraint's tangent space, $\mathbf{f}_c$ defined this way is always zero work for any valid velocity consistent with the constraints.) The constraint $\Cbb(\mathbf{q}) = \mathbf{0}$ induces two additional constraints by its time derivatives $\dot{\Cbb} = \mathbf{J}{\dot{\q}} = \mathbf{0}$ and $\ddot{\Cbb} = \mathbf{J}{\ddot{\q}} + {\dot{\J}}{\dot{\q}} = \mathbf{0}$. The constraint and its first derivative are satisfied by assumption, so it's the third that enters into the system:
\begin{align*}
    &\mathbf{M}{\ddot{\q}} + \mathbf{f} + \mathbf{J}^\tr\bm{\lambda} = \mathbf{0}\\
    &\mathbf{J}{\ddot{\q}} + {\dot{\J}}{\dot{\q}} = \mathbf{0},
\end{align*}
where we've captured all force terms in $\mathbf{f} = \bm{\xi} + \mathbf{f}_f + \partial\psi + \mathbf{B}{\dot{\q}}$. Multiplying the first equation by $\mathbf{J}\mathbf{M}^{-1}$ and subtracting the second, we can solve for $\bm{\lambda}$:
\begin{align*}
    \bm{\lambda} = (\mathbf{J}\mathbf{M}^{-1}\mathbf{J}^\tr)^{-1}\big[-\mathbf{J}\mathbf{M}^{-1}\mathbf{f}+{\dot{\J}}{\dot{\q}}\big].
\end{align*}
These Lagrange multipliers define the constraint forces; plugging them back into the first equation and collecting terms gives:
\begin{align*}
    \mathbf{M}{\ddot{\q}} + \mathbf{P}_{\!/\mkern-5mu/\!}\mathbf{f} + \mathbf{P}_\perp\mathbf{M}{\dot{\J}}{\dot{\q}} = \mathbf{0}.
\end{align*}
This result gives Equation~\ref{eqn:ConstrainedGeometricFabrics} once $\mathbf{f}$ is split apart. Additionally, adding and subtracting $\bm{\xi}$ gives
\begin{align*}
    &\mathbf{M}{\ddot{\q}} + \bm{\xi} - \bm{\xi} + \mathbf{P}_{\!/\mkern-5mu/\!}\big(\bm{\xi}+\mathbf{f}_f\big) + \mathbf{P}_\perp \mathbf{M}{\dot{\J}}{\dot{\q}} \\
    &\ \ \
    = \mathbf{M}{\ddot{\q}} + \bm{\xi} + \mathbf{P}_\perp\big(\mathbf{M}{\dot{\J}}{\dot{\q}} - \bm{\xi}\big) + \mathbf{P}_{\!/\mkern-5mu/\!}\mathbf{f}_f,
\end{align*}
which derives the expression for $\wt{\mathbf{f}}_f = \mathbf{P}_\perp\big(\mathbf{M}{\dot{\J}}{\dot{\q}} - \bm{\xi}\big) + \mathbf{P}_{\!/\mkern-5mu/\!}\mathbf{f}_f$

Observing that $\mathbf{M}$ is HD0 (see discussion in Section~\ref{sec:FinslerAndHd2Geometries}), we see that both $\mathbf{P}_{\!/\mkern-5mu/\!}$ and $\mathbf{P}_\perp$ are HD0 (i.e. they do not depend on $\|{\dot{\q}}\|$). Since $\mathbf{f}_f$ is HD2 by definition, its projection is as well. Moreover, $\bm{\xi}$ is known to be HD2 and using tensor notation, we see ${\dot{\J}}{\dot{\q}} = \frac{\partial^2\Cbb}{\partial q^i\partial q^j}\dot{q}^i\dot{q}^j$ is HD2, which means means both $\mathbf{M}{\dot{\J}}{\dot{\q}} - \bm{\xi}$ and its projection are HD2 as well. Together these results show that $\wt{\mathbf{f}}_f$ is HD2, and since we know it derives from the original bending term by adding a zero-work constraint force, $\wt{\mathbf{f}}_f$ must be zero work as well and is hence a bending term. This proves that $(\mathcal{L}_e, \wt{\mathbf{f}}_f)$ is a geometric fabric.
\end{proof}

\subsection{Theorem~\ref{thm:GeometricFabricStability} on Geometric Fabric Stability}
\label{apx:TheoremGeometricFabricStability}

Proof of theorem in Subsection~\ref{sec:StabilityAnalysis} characterizing the stability of geometric fabrics and convergence to KKT solutions of the constrained potential.
\begin{proof}
By Lemma~\ref{lma:EnergyConservation}, the forced system will reduce total energy $\mathcal{H}(\mathbf{q}, {\dot{\q}}) = \mathcal{L}_e(\mathbf{q}, {\dot{\q}}) + \psi(\mathbf{q})$, at a rate $\dot{\mathcal{H}} \leq -{\dot{\q}}^\tr\mathbf{B}(\mathbf{q},{\dot{\q}}){\dot{\q}}$ (the inequality is due to the impact model). Since $\mathcal{L}_e\geq 0$ and $\psi(\mathbf{q})$ is lower bounded, $\mathcal{H}$ is lower bounded and it must be that $\dot{\mathcal{H}}\rightarrow 0$. Therefore, ${\dot{\q}}^\tr\mathbf{B}(\mathbf{q},{\dot{\q}}){\dot{\q}}\rightarrow 0$ which implies ${\dot{\q}}\rightarrow\mathbf{0}$ since $\mathbf{B}$ is positive definite. Since velocity approaches zero we must also have ${\ddot{\q}}\rightarrow\mathbf{0}$. At convergence, a subset of inequality constraints will be active resulting in an equality constraint for those dimensions which we denote $\wt{\mathbf{g}} = \mathbf{0}$. Combining those active inequality constraints with the equality constraints $\mathbf{c}$ into $\Cbb = [\mathbf{c};\wt{\mathbf{g}}]$ with Jacobian $\mathbf{J} = \partial_\mathbf{q}\Cbb$, the constrained the equations of motion at convergence can be expressed
\begin{align}
    \mathbf{M}{\ddot{\q}} + \bm{\xi} + \mathbf{f}_f + \mathbf{J}^\tr\bm{\lambda}= -\partial_\mathbf{q}\psi - \mathbf{B}{\dot{\q}}.
\end{align}
Since ${\ddot{\q}}\rightarrow\mathbf{0}$ we have $\mathbf{M}{\ddot{\q}}\rightarrow\mathbf{0}$ ($\mathbf{M}$ is HD0) and since ${\dot{\q}}\rightarrow\mathbf{0}$ we have $\bm{\xi}+\mathbf{f}_f \rightarrow\mathbf{0}$ (both $\bm{\xi}$ and $\mathbf{f}_f$ are HD2). That leaves just $\mathbf{J}^\tr\bm{\lambda}$ on the left hand side. On the right hand side we have $\mathbf{B}{\dot{\q}}\rightarrow\mathbf{0}$ since $\mathbf{B}$ is bounded, so the system approaches
\begin{align}
    \mathbf{J}^\tr\bm{\lambda} = -\partial_\mathbf{q}\psi.
\end{align}
Therefore, it approaches a KKT solution since the system evolution constraints are everywhere satisfied and contact models require $\bm{\lambda}$ satisfy the required complimentarity conditions.
\end{proof}

\subsection{Theorem~\ref{thm:SystemEnergization} on Energized Fabrics}
\label{apx:TheoremOnEnergizedFabrics}

Proof of Theorem~\ref{thm:SystemEnergization} in Subsection~\ref{sec:EnergizingHd2Geometries} showing how to energize arbitrary HD2 geometries and proving the resulting system is a geometric fabric.
\begin{proof}
We derive these results using notation $\mathbf{h} = -\pi(\mathbf{q}, {\dot{\q}})$ so that the HD2 geometry can be expressed ${\ddot{\q}} + \mathbf{h}(\mathbf{q},{\dot{\q}}) = \mathbf{0}$. Using $\mathbf{M} {\ddot{\q}} + \bm{\xi} = \mathbf{0}$ to denote the unforced equations of motion,
Equation~\ref{eqn:EnergyTimeDerivative} of Appendix Lemma~\ref{lma:EnergyTimeDerivative} gives the time derivative of the energy $\dot{\mathcal{H}}_e = {\dot{\q}}^\tr\big(\mathbf{M}{\ddot{\q}} + \bm{\xi}\big)$. Substituting a system of the form ${\ddot{\q}} = -\mathbf{h}(\mathbf{q},{\dot{\q}}) - \alpha_{\mathcal{H}_e}{\dot{\q}}$, setting it to zero, and solving for $\alpha_{\mathcal{H}_e}$ gives
\begin{align*}
&\dot{\mathcal{H}}_e = {\dot{\q}}^\tr\Big[\mathbf{M}\big(-\mathbf{h} - \alpha_{\mathcal{H}_e}{\dot{\q}}\big) + \bm{\xi}\Big] = \mathbf{0} \\
&\Rightarrow -{\dot{\q}}^\tr\mathbf{M}\mathbf{h} - \alpha_{\mathcal{H}_e} {\dot{\q}}^\tr\mathbf{M}{\dot{\q}} + {\dot{\q}}^\tr\bm{\xi} = \mathbf{0} \\
&\Rightarrow \alpha_{\mathcal{H}_e} = -\frac{{\dot{\q}}^\tr\mathbf{M}\mathbf{h} - {\dot{\q}}^\tr\bm{\xi}}{{\dot{\q}}^\tr\mathbf{M}{\dot{\q}}}\\
&\ \ \ \ \ \ \ \ \ \ =
-({\dot{\q}}^\tr\mathbf{M}{\dot{\q}})^{-1}{\dot{\q}}^\tr\big[\mathbf{M}\mathbf{h} - \bm{\xi}\big].
\end{align*}
This result gives the formula for Equation~\ref{eqn:EnergizationTransformAlpha}. Substituting this solution for $\alpha_{\mathcal{H}_e}$ back in gives
\begin{align*}
    {\ddot{\q}}
    &= -\mathbf{h} + \left(\frac{{\dot{\q}}^\tr}{{\dot{\q}}^\tr\mathbf{M}{\dot{\q}}}\big[\mathbf{M}\mathbf{h} - \bm{\xi}\big]\right) {\dot{\q}} \\
    &= -\mathbf{h} + \left[\frac{{\dot{\q}}\,{\dot{\q}}^\tr}{{\dot{\q}}^\tr\mathbf{M}{\dot{\q}}}\right]\big(\mathbf{M}\mathbf{h} - \bm{\xi}\big).
\end{align*}
Algebraically, it helps to introduce $\mathbf{h} = \mathbf{M}^{-1}\big[\mathbf{f}_f + \bm{\xi}\big]$ to first express the result as a difference away from $\bm{\xi}$; in the end we will convert back to $\mathbf{h}$. Doing so and moving all the terms to the left hand side of the equation gives
\begin{align*}
    &\ \ \ \ {\ddot{\q}} + \mathbf{h} - \left[\frac{{\dot{\q}}\,{\dot{\q}}^\tr}{{\dot{\q}}^\tr\mathbf{M}{\dot{\q}}}\right]\big(\mathbf{M}\mathbf{h} - \bm{\xi}\big) = \mathbf{0} \\
    &\Rightarrow {\ddot{\q}} + \mathbf{M}^{-1}\big[\mathbf{f}_f + \bm{\xi}\big] \\
    &\ \ \ \ \ - \left[\frac{{\dot{\q}}\,{\dot{\q}}^\tr}{{\dot{\q}}^\tr\mathbf{M}{\dot{\q}}}\right]\Big(\mathbf{M}\M^{-1}\big[\mathbf{f}_f + \bm{\xi}\big] - \bm{\xi}\Big) = \mathbf{0} \\
    &\Rightarrow \mathbf{M} {\ddot{\q}} + \mathbf{f}_f + \bm{\xi} \\
    &\ \ \ \ \ - \mathbf{M}\left[\frac{{\dot{\q}}\,{\dot{\q}}^\tr}{{\dot{\q}}^\tr\mathbf{M}{\dot{\q}}}\right]\Big(\mathbf{f}_f + \big(\bm{\xi} - \bm{\xi}\big)\Big) = \mathbf{0} \\
    &\Rightarrow \mathbf{M} {\ddot{\q}} + \bm{\xi} + \mathbf{M}\left[\mathbf{M}^{-1} - \frac{{\dot{\q}}\,{\dot{\q}}^\tr}{{\dot{\q}}^\tr\mathbf{M}{\dot{\q}}}\right] \mathbf{f}_f = \mathbf{0} \\
    &\Rightarrow \mathbf{M} {\ddot{\q}} + \bm{\xi} + \mathbf{M}\mathbf{R}_{\mathbf{p}}\big(\mathbf{M}\mathbf{h} - \bm{\xi}\big) = \mathbf{0},
\end{align*}
where $\mathbf{R}_{\mathbf{p}} = \mathbf{M}^{-1} - \frac{{\dot{\q}}\,{\dot{\q}}^\tr}{{\dot{\q}}^\tr\mathbf{M}{\dot{\q}}}$ matching Equation~\ref{eqn:Rpe} of Lemma~\ref{lma:EnergyProjection}, and we substitute $\mathbf{f}_f = \mathbf{M}\mathbf{h} - \bm{\xi}$ back in. By Lemma~\ref{lma:EnergyProjection} $\mathbf{M}\mathbf{R}_{\mathbf{p}} = \mathbf{P}_e$, so we get Equation~\ref{eqn:ZeroWorkEnergizationForm}.
Moreover, again by Lemma~\ref{lma:EnergyProjection}, ${\dot{\q}}^\tr\mathbf{P}_e = \mathbf{0}$, so ${\dot{\q}}^\tr\mathbf{f}_f = -{\dot{\q}}^\tr\mathbf{P}_e\big[\mathbf{M}\pi + \bm{\xi}\big] = 0$ showing that the new term is a zero work modification. Finally, both $\mathbf{h}$ and $\bm{\xi}$ are HD2 by definition, and since $\mathbf{P}_e$ and $\mathbf{M}$ are HD0, the modifying term is HD2 making it a bending term. Therefore, the energized system is a bent Finsler system and hence a geometric fabric.
\end{proof}

\subsection{Theorem~\ref{thm:FabricsInGeneralizedCoords} on Fabrics in Generalized Coordinates}
\label{apx:TheoremFabricsInGeneralizedCoords}

Proof of Theorem~\ref{thm:FabricsInGeneralizedCoords} showing how to express in generalized coordinates energized systems that are subsequently constrained.
See Subsection~\ref{sec:GeneralizedCoordinates}.

\begin{proof}
We will show the equivalence by calculation. Let $\mathbf{h}(\mathbf{x},{\dot{\x}}) = -\pi(\mathbf{x},{\dot{\x}})$ so that the geometry's equation becomes ${\ddot{\x}} + \mathbf{h}(\mathbf{x},{\dot{\x}}) = \mathbf{0}$. Multiplying through by $\mathbf{M}$ gives its metric-weighted (force) form $\mathbf{M}{\ddot{\x}} + \mathbf{f} = \mathbf{0}$ where $\mathbf{f} = \mathbf{M}\mathbf{h}$.
Using Equation~\ref{eqn:EnergizationTransformAlpha} of Theorem~\ref{thm:SystemEnergization}, this system's energization with respect to energy $\mathcal{L}_e$ with equations of motion $\mathbf{M}{\ddot{\x}} + \bm{\xi} = \mathbf{0}$ can be expressed
$\mathbf{M}{\ddot{\x}} + \bm{\xi}^{\mathbf{h}} = \mathbf{0}$ where
\begin{align}
    \bm{\xi}^{\mathbf{h}} = \bm{\xi} + \mathbf{M}\left[\mathbf{M}^{-1} - \frac{{\dot{\x}}\xd^\tr}{{\dot{\x}}^\tr\mathbf{M}{\dot{\x}}}\right] \big(\mathbf{f} - \bm{\xi}\big).
\end{align}
Let $\mathbf{J} = \partial_\mathbf{q}\phi$. The pullback of the energized geometry generator is
\begin{align}\nonumber
    &\mathbf{J}^\tr\mathbf{M}\left(\mathbf{J}{\ddot{\q}} + {\dot{\J}}{\dot{\q}}\right) + \mathbf{J}^\tr\bm{\xi}^{\mathbf{h}} = \mathbf{0} \\\nonumber
    &\Rightarrow \big(\mathbf{J}^\tr\mathbf{M}\mathbf{J}\big){\ddot{\q}} + \mathbf{J}^\tr\big(\bm{\xi}^{\mathbf{h}} + \mathbf{M}{\dot{\J}}{\dot{\q}}\big) = \mathbf{0}\\\nonumber
    &\Rightarrow \big(\mathbf{J}^\tr\mathbf{M}\mathbf{J}\big){\ddot{\q}} + \mathbf{J}^\tr \bm{\xi} \\\nonumber
    &\ \ \ \ + \mathbf{J}^\tr\mathbf{M}\left[\mathbf{M}^{-1} - \frac{{\dot{\x}}\xd^\tr}{{\dot{\x}}^\tr\mathbf{M}{\dot{\x}}}\right]
            \big(\mathbf{f} - \bm{\xi}\big)
        + \mathbf{J}^\tr\mathbf{M}{\dot{\J}}{\dot{\q}} = \mathbf{0} \\\label{eqn:EnergizationPullback}
    &\Rightarrow \wt{\mathbf{M}}{\ddot{\q}} + \wt{\bm{\xi}}
        + \mathbf{J}^\tr\mathbf{M}\left[\mathbf{M}^{-1} - \frac{{\dot{\x}}\xd^\tr}{{\dot{\x}}^\tr\mathbf{M}{\dot{\x}}}\right]
            \big(\mathbf{f} - \bm{\xi}\big),
\end{align}
where $\wt{\mathbf{M}} = \mathbf{J}^\tr\mathbf{M}\mathbf{J}$ and $\wt{\bm{\xi}} = \mathbf{J}^\tr\big(\bm{\xi} + \mathbf{M}{\dot{\J}}{\dot{\q}}\big)$ form the standard pullback of $(\mathbf{M}, \bm{\xi})$.

Now we show that pulling back the geometry first then energizing it in $\mathbf{q}$ results in a system matching Equation~\ref{eqn:EnergizationPullback}. We can calculate the geometry pullback with respect to the energy metric $\mathbf{M}$ by pulling back the metric weighted force form of the geometry $\mathbf{M}{\ddot{\x}} + \mathbf{f} = \mathbf{0}$, where again $\mathbf{f} = \mathbf{M}\mathbf{h}$. The pullback is
\begin{align}\nonumber
    &\mathbf{J}^\tr\mathbf{M}\big(\mathbf{J}{\ddot{\q}} + {\dot{\J}}{\dot{\q}}\big) + \mathbf{J}^\tr\mathbf{f} = \mathbf{0} \\\nonumber
    &\Rightarrow \big(\mathbf{J}^\tr\mathbf{M}\mathbf{J}\big) {\ddot{\q}} + \mathbf{J}^\tr\big(\mathbf{f} + \mathbf{M}{\dot{\J}}{\dot{\q}}\big) \\\label{eqn:PullbackGeometry}
    &\Leftrightarrow \wt{\mathbf{M}}{\ddot{\q}} + \wt{\mathbf{f}} = \mathbf{0}
\end{align}
where $\wt{\mathbf{M}} = \mathbf{J}^\tr\mathbf{M}\mathbf{J}$ as before and $\wt{\mathbf{f}} = \mathbf{J}^\tr \big(\mathbf{f} + \mathbf{M}{\dot{\J}}{\dot{\q}}\big)$.

Let $\wt{\mathcal{L}}_e = \mathcal{L}_e\big(\phi(\mathbf{q}), \mathbf{J}{\dot{\q}}\big)$ be the pullback of the energy function $\mathcal{L}_e$. We know that the Euler-Lagrange equation commutes with the pullback, so applying the Euler-Lagrange equation to this pullback energy $\wt{\mathcal{L}}_e$ is equivalent to pulling back the Euler-Lagrange equation of $\mathcal{L}_e$. This means we can calculate the Euler-Lagrange equation of $\mathcal{L}_e$ as
\begin{align}
    &\big(\mathbf{J}^\tr\mathbf{M}\mathbf{J}\big) {\ddot{\q}} + \mathbf{J}^\tr\big(\bm{\xi} + \mathbf{M}{\dot{\J}}{\dot{\q}}\big) = \mathbf{0}\\
    &\Leftrightarrow \wt{\mathbf{M}}{\ddot{\q}} + \wt{\bm{\xi}} = \mathbf{0},
\end{align}
with $\wt{\mathbf{M}} = \mathbf{J}^\tr\mathbf{M}\mathbf{J}$ and $\wt{\bm{\xi}} = \mathbf{J}^\tr\big(\bm{\xi} + \mathbf{M}{\dot{\J}}{\dot{\q}}\big)$ (both as previously defined). Therefore, energizing \ref{eqn:PullbackGeometry} with $\wt{\mathcal{L}}_e$ gives
\begin{align}\nonumber
    &\wt{\mathbf{M}}{\ddot{\q}} + \wt{\bm{\xi}}
      + \wt{\mathbf{M}}\left[
        \wt{\mathbf{M}}^{-1} - \frac{{\dot{\q}}\qd^\tr}{{\dot{\q}}^T\wt{\mathbf{M}}{\dot{\q}}}\right]
        \big(\wt{\mathbf{f}} - \wt{\bm{\xi}}\big) = \mathbf{0} \\\nonumber
    &\Rightarrow
    \wt{\mathbf{M}}{\ddot{\q}} + \wt{\bm{\xi}} \\\nonumber
    &\ \ \ \ + \big(\mathbf{J}^\tr\mathbf{M}\mathbf{J}\big)
          \left[
            \wt{\mathbf{M}}^{-1} - \frac{{\dot{\q}}\qd^\tr}{{\dot{\q}}^\tr \mathbf{J}^\tr\mathbf{M}\mathbf{J}{\dot{\q}}}
          \right] \cdot\\\nonumber
    &\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
        \Big(\mathbf{J}^\tr\big(\mathbf{f} + \mathbf{M}{\dot{\J}}{\dot{\q}}\big)
          - \mathbf{J}^\tr\big(\bm{\xi} + \mathbf{M}{\dot{\J}}{\dot{\q}}\big)
        \Big)
      = \mathbf{0} \\\nonumber
    &\Rightarrow
    \wt{\mathbf{M}}{\ddot{\q}} + \wt{\bm{\xi}} \\\nonumber
    &\ \ \ \ + \big(\mathbf{J}^\tr\mathbf{M}\big)
          \mathbf{J}\left[
            \big(\mathbf{J}^\tr\mathbf{M}\mathbf{J}\big)^{-1} - \frac{{\dot{\q}}\qd^\tr}{{\dot{\x}}^\tr\mathbf{M}{\dot{\x}}}
          \right]\mathbf{J}^\tr
        \big(\mathbf{f}- \bm{\xi}\big) = \mathbf{0}\\\label{eqn:PullbackEnergizationAlmost}
    &\Rightarrow
    \wt{\mathbf{M}}{\ddot{\q}} + \wt{\bm{\xi}} \\\nonumber
    &\ \ \ \ + \mathbf{J}^\tr\mathbf{M}
          \left[
            \mathbf{J}\big(\mathbf{J}^\tr\mathbf{M}\mathbf{J}\big)^{-1}\mathbf{J}^\tr - \frac{{\dot{\x}}\xd^\tr}{{\dot{\x}}^\tr\mathbf{M}{\dot{\x}}}
          \right]
        \big(\mathbf{f}- \bm{\xi}\big) = \mathbf{0}.
\end{align}
Since $\phi$ is full rank $\mathbf{J}^\tr\mathbf{M}\mathbf{J}$ is invertible, so
\begin{align*}
    \mathbf{J}^\tr\mathbf{M}\mathbf{J}\big(\mathbf{J}^\tr\mathbf{M}\mathbf{J}\big)^{-1}\mathbf{J}^\tr
    = \mathbf{J}^\tr = \mathbf{J}^\tr\mathbf{M}\big(\mathbf{M}^{-1}\big),
\end{align*}
which means we can rewrite Equation~\ref{eqn:PullbackEnergizationAlmost} as
\begin{align}\label{eqn:PullbackEnergized}
    \wt{\mathbf{M}}{\ddot{\q}} + \wt{\bm{\xi}}
      + \mathbf{J}^\tr\mathbf{M}
          \left[
            \mathbf{M}^{-1} - \frac{{\dot{\x}}\xd^\tr}{{\dot{\x}}^\tr\mathbf{M}{\dot{\x}}}
          \right]
        \big(\mathbf{f}- \bm{\xi}\big) = \mathbf{0}.
\end{align}
This final expression in Equation~\ref{eqn:PullbackEnergized} matches the expression for the energized geometry pullback in Equation~\ref{eqn:EnergizationPullback}.
\end{proof}

\section{Details on Designing with Fabrics}

\subsection{Derivations of Fabric Component Algebra}
\label{apx:DerivationsOfFabricComponentAlgebra}

A fabric component was defined in Subsection~\ref{sec:AlgebraOfComponents} as a Finsler energy paired with a force term $\{\mathcal{L}_e, \mathbf{f}\}$. We often refer to this pair as a Finsler energy paired with a {\em HD2 geometry} $\pi(\mathbf{x},{\dot{\x}})$ knowing that the actual component would be $\mathbf{f} = -\mathbf{M}\pi$ where $\mathbf{M}$ is the mass from the energy's equations of motion $\mathbf{M}{\ddot{\x}} + \bm{\xi} = \mathbf{0}$. Since the energization formulas in Theorem~\ref{thm:SystemEnergization} are functions of the Finsler energy $\mathcal{L}_e$ only through its equations of motion parameterized by $(\mathbf{M},\bm{\xi})$, we can think of the fabric component as effectively a triple $(\mathbf{M}, \bm{\xi},\mathbf{f})$ representing the energy equations alongside an additional forcing term, which itself can be viewed as two pairs $\big((\mathbf{M}, \bm{\xi}), (\mathbf{M},\mathbf{f})\big)$. The algebra defined in Subsection~\ref{sec:AlgebraOfComponents} is equivalent to the reduction of sums and pullback to the sums and pullback of these individual pairs. For concreteness, and in alignment with the terminology of \cite{optimizationFabricsForBehavioralDesignArXiv2020}, we call these pairs {\em specs}, short for {\em spectral semi-sprays}, which is descriptive of the differential equation they represent emphasizing the spectral role of the metric $\mathbf{M}$.

The one operation on these specs not discussed in Subsection~\ref{sec:AlgebraOfComponents} is summation and the linearity of summation. This property can be easily seen by
\begin{align*}
    &\alpha(\mathbf{M}_1,\mathbf{f}_1) + \beta(\mathbf{M}_2,\mathbf{f}_2) \\
    &\ \ \Leftrightarrow
    \alpha \big(\mathbf{M}_1{\ddot{\x}} + \mathbf{f}_1\big) + \beta \big(\mathbf{M}_2{\ddot{\x}} + \mathbf{f}_2\big) = \mathbf{0} \\
    &\Rightarrow
    \big(\alpha\mathbf{M}_1 + \beta\mathbf{M}_2\big){\ddot{\x}} + \big(\alpha\mathbf{f}_1 + \beta\mathbf{f}_2\big) = \mathbf{0} \\
    &\ \ \Leftrightarrow \Big(\big(\alpha\mathbf{M}_1 + \beta\mathbf{M}_2\big), \big(\alpha\mathbf{f}_1 + \beta\mathbf{f}_2\big)\Big).
\end{align*}
As mentioned in the discussion leading up to Equation~\ref{eqn:PullbackCommutesWithEL}, the pullback of energy equations is equivalent to applying the Euler-Lagrange equation on the pullback energy. Since application of the Euler-Lagrange equation is linear in the sense $\mathrm{EL}\big[\alpha \mathcal{L}_1 + \beta \mathcal{L}_2\big] = \alpha\,\mathrm{EL}\big[\mathcal{L}_1\big] + \beta\,\mathrm{EL}\big[\mathcal{L}_2\big]$,
if we treat any two energies that differ by only a constant as equivalent, the mapping defined by the $\mathrm{EL}[\cdot]$ operator defines an isomorphism between energies and the subspace of specs defined by those energies. That means, algebraically, we can treat the two representations as equivalent.

That algebra of specs derives the summation operation on fabric components. For the pullback operation, we can simply note that by Equation~\ref{eqn:PullbackCommutesWithEL}, the first (energy) spec of the component has a pullback matching the application of the Euler-Lagrange equation applied to the pullback energy, so the pullback of the energy captures that part. It additionally captures the pullback of the mass matrix of the second spec, so we need only additionally account for how that second spec's force term pulls back. Looking at the definition of component pullback in Equation~\ref{eqn:SpecPullback} we see that that's exactly what the force term pullback expression computes.

Finally, we note that while the definition of fabric components and its associated algebra as presented in Subsection~\ref{sec:AlgebraOfComponents} is compact and gives insight into what fabric components are and how they interact, it may be more efficient in practice to represent them in compact spec form explicitly as triples $(\mathbf{M},\bm{\xi},\mathbf{f})$ in code. That said, leveraging automatic differentiation frameworks such as those used in \cite{LiRmp2RobotLearning2021} may lead to alternative efficient representations.

We note that some fabric components may represent potential terms with dampers as their $\mathbf{f}$ term. These components intermix with the other fabric components using the same algebra since a component is defined independent of $\mathbf{f}$ semantics.

\subsection{Formulas for speed control}
\label{apx:FormulasSpeedControl}

This section gives explicit formulas for calculating the damping $\beta(\mathbf{q},{\dot{\q}})$ described in Subsection~\ref{sec:SpeedControl} for regulating a specific measure of speed.

Once a geometry is energized by a Finsler energy $\mathcal{L}_e$, its speed profile is defined by what is necessary to conserve $\mathcal{L}_e$. In practice, we usually want to regulate a different measure of energy (speed), an {\em execution energy} $\mathcal{L}_e^{\mathrm{ex}}$, such as Euclidean energy in either C-space or end-effector space (or, more commonly, some combination of both).
Denoting the fabric by $(\mathcal{L}_e, -\mathbf{M}\pi - \bm{\xi})$ which defines nominal fabric equations
\begin{align*}
    &\mathbf{M} {\ddot{\q}} + \bm{\xi} + \mathbf{f}_f \\
    &\ \ = \mathbf{M} {\ddot{\q}} + \bm{\xi} - \mathbf{M}\pi - \bm{\xi} \\
    &\ \ = \mathbf{M} ({\ddot{\q}} - \pi) = \mathbf{0}
\end{align*}
so that the resolved behavior matches the geometric policy ${\ddot{\q}} = \pi$,
we can regulate the speed by adjusting $\alpha_\mathrm{reg}$ in a forced and damped equation of the form
\begin{align} \label{eqn:FundamentalsOfSpeedControl}
    {\ddot{\q}} = -\mathbf{M}^{-1}\partial_\mathbf{q}\psi(\mathbf{q}) + \pi(\mathbf{q}, {\dot{\q}}) + \alpha_\mathrm{reg} {\dot{\q}}.
\end{align}
As long as $\alpha_\mathrm{reg} < \alpha_{\mathcal{L}_e}$, where $\alpha_{\mathcal{L}_e}$ is the energization coefficient (Equation~\ref{eqn:EnergizationTransformAlpha} of Theorem~\ref{thm:SystemEnergization}), the system will remain stable per Theorem~\ref{thm:GeometricFabricStability} since $\alpha_\mathrm{reg}$ always act as a nonzero damper with respect to the energized system. Note that $\alpha_{\mathcal{L}_e}$ need not always be smaller than $0$, so we have the flexibility with $\alpha_\mathrm{reg}$ to both remove and add energy as needed (although, we don't necessarily have to use that flexibility; the method outlined below restricts it to being positive to aid intuition). Below we calculate explicit formulas for $\alpha_\mathrm{reg}$ to regulate a given measure of speed (energy).

Let $\mathcal{L}_e^\mathrm{ex}$ be an {\em execution} energy, which may differ from the fabric's energy $\mathcal{L}_e$, and let $\alpha_\mathrm{ex}^0$ and $\alpha_\mathrm{ex}^\psi$ be the energization coefficients, respectively, for $\mathrm{energize}_{\mathcal{L}_e^\mathrm{ex}}\big[\pi\big]$ and $\mathrm{energize}_{\mathcal{L}_e^\mathrm{ex}}\big[-\mathbf{M}^{-1}\partial_\mathbf{q}\psi + \pi\big]$. The coefficient $\alpha_\mathrm{ex}^0$ ensures that system
\begin{align}
    {\ddot{\q}} = \pi(\mathbf{q},{\dot{\q}}) + \alpha_\mathrm{ex}^0{\dot{\q}}
\end{align}
conserves $\mathcal{L}_e^{\mathrm{ex}}$ and the coefficient $\alpha_\mathrm{ex}^\psi$ ensures that the {\em forced} system\footnote{As noted above, the energization coefficient calculation in Theorem~\ref{thm:SystemEnergization} does not require that the underlying differential equation being energized be a geometry (in this case it's a geometry forced by a potential). It simply states that if the underlying differential equation is indeed a geometry, then the resulting energized equation is a geometric fabric. Here we're just using the result to calculate a reference coefficient.}
\begin{align}
    {\ddot{\q}} = -\mathbf{M}^{-1}\partial_\mathbf{q}\psi + \pi(\mathbf{q},{\dot{\q}}) + \alpha_\mathrm{ex}^\psi{\dot{\q}}
\end{align}
conserves $\mathcal{L}_e^{\mathrm{ex}}$.
Likewise, let $\alpha_{\mathcal{L}_e}$ denote the energy coefficient of the actual system energy $\mathrm{energize}_{\mathcal{L}_e}\big[\pi\big]$ (designed to conserve $\mathcal{L}_e$).

The forced and damped energized system
\begin{align}
    {\ddot{\q}} = -\mathbf{M}^{-1}\partial_\mathbf{q}\psi + \pi(\mathbf{q},{\dot{\q}}) + \alpha_{\mathcal{L}_e}{\dot{\q}} - \wt{\beta} {\dot{\q}},
\end{align}
would be stable for $\wt{\beta}>0$. Adding and subtracting $\alpha_\mathrm{ex}^0$ gives
\begin{align*}
    {\ddot{\q}}
    &= -\mathbf{M}^{-1}\partial_\mathbf{q}\psi + \pi(\mathbf{q},{\dot{\q}}) + \big(\alpha_{\mathcal{L}_e} + \alpha_\mathrm{ex}^0 - \alpha_\mathrm{ex}^0\big){\dot{\q}} - \wt{\beta} {\dot{\q}} \\
    &= -\mathbf{M}^{-1}\partial_\mathbf{q}\psi + \pi(\mathbf{q},{\dot{\q}}) + \alpha_\mathrm{ex}^0{\dot{\q}} - \beta {\dot{\q}}\\
    &= -\mathbf{M}^{-1}\partial_\mathbf{q}\psi + \mathrm{energize}_{\mathcal{L}_\mathrm{ex}}\big[\pi(\mathbf{q},{\dot{\q}})\big] - \beta {\dot{\q}},
\end{align*}
where $\beta = \wt{\beta} + \alpha_\mathrm{ex}^0 - \alpha_{\mathcal{L}_e}$. The constraint $\wt{\beta}>0$ written this way manifests as $\beta > \alpha_\mathrm{ex}^0 - \alpha_{\mathcal{L}_e}$.

This calculation suggests we can regulate speed by choosing $\alpha_\mathrm{reg}$ in Equation~\ref{eqn:FundamentalsOfSpeedControl} as
\begin{align} \label{eqn:AlphaRegForm}
    \alpha_\mathrm{reg}
    = \alpha_\mathrm{ex}^\eta - \beta_\mathrm{reg}(\mathbf{q}, {\dot{\q}}) + \alpha_\mathrm{boost}
\end{align}
where $\alpha_\mathrm{ex}^\eta = \eta \alpha_\mathrm{ex}^0 + (1-\eta)\alpha_\mathrm{ex}^\psi$ for $\eta\in[0,1]$ and $\alpha_\mathrm{boost} \leq 0$ as long as $\beta_\mathrm{reg} > \alpha_\mathrm{ex}^0 - \alpha_{\mathcal{L}_e}$.
The term $\alpha_\mathrm{boost}$ can be used temporarily in the beginning to boost the system up to speed. Since it is transient, we drop $\alpha_\mathrm{boost}$ momentarily for simplicity, but return to it at the end of the section to describe a good boosting policy. In our experiments, we use that term only for the cubby navigation problem in order to reduce the strength of the potential function near the beginnings of motions to get a clear signal on what the geometry itself is encoding.

Under this choice of $\alpha_\mathrm{reg}$ given in Equation~\ref{eqn:AlphaRegForm}, using $\mathbf{v}_{\!/\mkern-5mu/\!} = \big(\alpha_\mathrm{ex}^\psi - \alpha_\mathrm{ex}^0\big){\dot{\q}}$ to denote the component of $-\mathbf{M}^{-1}\partial_\mathbf{q}\psi$ additionally removed by including it within the energization operation (and dropping $\alpha_\mathrm{boost}$), we can express the system as
\begin{align*}
    {\ddot{\q}} &= \eta \mathbf{v}_{\!/\mkern-5mu/\!} + \mathrm{energize}_{\mathcal{L}_e^\mathrm{ex}}\big[-\mathbf{M}^{-1}\partial_\mathbf{q}\psi + \pi\big] - \beta_\mathrm{reg}(\mathbf{q}, {\dot{\q}}){\dot{\q}},
\end{align*}
using a similar analysis to the above. Intuitively, the energized system will conserve execution energy, $\eta \mathbf{v}_{\!/\mkern-5mu/\!}$ can inject energy into the system to speed it up (by an amount parameterized by $\eta\in[0,1]$), and $\beta_\mathrm{reg}$ can bleed off energy. Note that the lower bound on the damper isn't strictly zero $\beta_\mathrm{reg} > \alpha_\mathrm{ex}^0 - \alpha_{\mathcal{L}_e}$, so there may be both cases where that term is allowed to inject extra energy and when that term is forced to be strictly a margin above zero for stability.

In practice, we enforce that $\beta_\mathrm{reg}$ is strictly positive by choosing $\beta_\mathrm{reg} = s_\beta(\mathbf{q}) B + \underline{B} + \max\{0, \alpha_\mathrm{ex}^\eta - \alpha_{\mathcal{L}_e}\} > 0$ with both $B$ and $\underline{B}$ positive. This gives us a more intuitive interface to speed control with $\eta$ adjusting energy in and $\beta_\mathrm{reg}$ adjusting energy out.
In this form, $\underline{B}$ acts as a constant baseline damping coefficient, and $s_\beta(\mathbf{q}) B$ provides additional damping near the convergence point with $B>0$ and $s_\beta(\mathbf{q})$ acting as a switch transitioning from 0 to 1 as the system approaches the target. $\max\{0, \alpha_\mathrm{ex}^\eta - \alpha_{\mathcal{L}_e}\}$ ensures the stability bound $\alpha_\mathrm{reg} = \beta_\mathrm{reg} \geq \alpha_\mathrm{ex} - \alpha_{\mathcal{L}_e}$ is satisfied with the baseline coefficients making it a strict inequality.

Specifically, in our experiments we use
\begin{align}
\label{eqn:damping_gate}
 s_\beta(\mathbf{q}) = \frac{1}{2} \Big(\tanh\big(-\alpha_\beta (\|\mathbf{q}\| - r) \big) + 1\Big)
\end{align}
where $\alpha_\beta \in \mathbb{R}^+$ is a gain defining the switching rate, and $r \in \mathbb{R}^+$ is the radius where the switch is half-way engaged. Denoting the desired execution energy as $\mathcal{L}_e^\mathrm{ex,d}$, we use the following policy for $\eta$
\begin{align}
\label{eqn:energy_gate}
    \eta = \frac{1}{2} \Big(\tanh\big(-\alpha_\eta (\mathcal{L}_e^\mathrm{ex} - \mathcal{L}_e^\mathrm{ex,d}) - \alpha_\mathrm{shift} \big) + 1\Big)
\end{align}
where $\alpha_\eta,\alpha_\mathrm{shift} \in \mathbb{R}^+$ adjust the rate and offset, respectively, of the switch as an affine function of the speed (execution energy) error.

When used, $\alpha_\mathrm{boost}$ acts to explicitly inject energy along the direction of motion to quickly boost the system up to speed. When ${\dot{\q}} = \mathbf{0}$, this term has no effect, so the initial direction of motion is chosen by the potential $\-\partial_\mathbf{q}\psi$. Once ${\dot{\q}} \neq \mathbf{0}$, $\alpha_\mathrm{boost}$ quickly accelerates the system to the desired speed. When $\pi$ is a geometry, quickly reaching a high speed means that the influence of the non-geometric potential is diminished, promoting path consistency. We reduce $\alpha_\mathrm{boost}$ to zero permanently once the system is close to convergence so it does not affect convergence stability.

In these cases, $\alpha_\mathrm{boost}$ is modeled as $\alpha_\mathrm{boost} = k \: \eta \big(1-s_\beta(\mathbf{q})\big)\frac{1}{\|{\dot{\q}}\| + \epsilon}$, where $k \in \mathbb{R}^+$ is a gain that directly sets the desired level of acceleration, $\eta$ (from above) sets $\alpha_\mathrm{boost}=0$ when the desired speed is achieved, and $1-s_\beta(\mathbf{q})$ sets $\alpha_\mathrm{boost}=0$ when the system is within the region of higher damping. The normalization by $\|{\dot{\q}}\|+\epsilon$ ensures that $\alpha_\mathrm{boost}$ is directly applied along $\hat{{\dot{\q}}}$ with a very small positive value for $\epsilon$ to ensure numerical stability. This overall design injects more energy into the system when $-\alpha_\mathrm{boost} < \alpha_\mathrm{ex}^\eta - \alpha_{\mathcal{L}_e}$. Since this injection occurs for finite time, the total system energy is still bounded. The additional switches ensure that the system is still subject to positive damping, guaranteeing convergence.

\subsection{Designing fabric components}
\label{apx:DesigningHD2Terms}

This section reviews a number of techniques for designing HD2 functions and maps for the design of Finsler energies and HD2 geometries referenced in Subsection~\ref{sec:ToolsForFabricDesign} on tools for fabric design.

The design of fabric components largely involves the construction of a geometric policy expressing a desired component behavior and a Finsler energy whose metric will act as that policy's priority metric. Both of these parts are defined by HD2 terms, so the question of design largely reduces to the question of how to design HD2 functions and maps.

As a reminder, a Homogeneous of Degree $k$  (HD$k$) function is a function $f(\mathbf{x})$ for which $f(\alpha \mathbf{x}) = \alpha^k f(\mathbf{x})$. In our case, our functions are typically HD0 through HD2, specifically in velocity. For instance, if $\pi(\mathbf{x},{\dot{\x}})$ is a geometric policy, the {\em geometric} modifier means it is HD2 in velocity, so $\pi(\mathbf{x},\alpha{\dot{\x}}) = \alpha^2 \pi(\mathbf{x},{\dot{\x}})$, and usually we assuming this holds just for $\alpha\geq 0$. It places no restrictions on how the function should scale with $\mathbf{x}$.

Often we design HD2 terms by multiplying an HD0 term by an HD2 scaling factor. That HD0 term might be entirely independent of velocity, or it may be dependent on the direction of the velocity, but not its norm (speed). Note that in the 1D case, the velocity direction is a switch, the sign of the scalar one-dimensional velocity, switching from 1 to -1 as velocity passes from positive to negative. This discrete switch is valid with our theory (systems have to come to a rest (zero velocity) before experiencing the switch) and is used commonly in our designs, especially for one-dimensional barrier geometries and energies.

In higher dimensions, switches depending on the norm of a velocity may not result in positive definite mass matrices. The presentation in \cite{ratliff2021geometry} allows that case by relaxing the third requirement of Definition~\ref{def:FinslerStructure} to require only invertibility, which is internally consistent with the geometric theory. But here, as a generalization of classical mechanics, we want mass matrices to remain positive definite, so building such switching terms in 1D spaces (where is naturally positive definite) and pulling them back to higher-dimensional spaces becomes an important technique for Finsler energy design.

A common pattern is to start with a scalar function $\psi(\mathbf{x})$ (often similar in structure and form to standard potential functions), and either use it directly as a position only base function, or compute a base policy as $\pi_0(\mathbf{x}) = -\partial_\mathbf{x} \psi$. In this discussion we will use the base policy as an example and construct a geometric policy. This base policy is HD0 simply because it doesn't depend on velocity at all. It's often useful to depend on velocity direction $\wh{{\dot{\x}}} = {\dot{\x}}/\|{\dot{\x}}\|$, though. Since $\wh{{\dot{\x}}}$ is already independent of speed $s = \|{\dot{\x}}\|$, any function of velocity direction is automatically HD0 as well. That gives $\pi_0^\sigma(\mathbf{x},{\dot{\x}}) = \sigma(\wh{{\dot{\x}}})\pi_0(\mathbf{x})$ as a new HD0 base. Finally, any Finsler energy $\mathcal{L}_e(\mathbf{x},{\dot{\x}}) = \frac{1}{2}\mathcal{L}_g^2(\mathbf{x},{\dot{\x}})$ is HD2 by definition, so scaling $\pi_0^\sigma$ by $\mathcal{L}_e$ gives us an HD2 geometry
\begin{align}
    \pi(\mathbf{x},{\dot{\x}}) = \mathcal{L}_e(\mathbf{x},{\dot{\x}})\,\pi_0^\sigma(\mathbf{x},{\dot{\x}}).
\end{align}
For instance, we might choose
\begin{align}
    \psi(\mathbf{x}) &= \frac{1}{2}\|\mathbf{x}_0 - \mathbf{x}\|^2
\end{align}
as a quadratic attractor to a target $\mathbf{x}_0$. The negative gradient of the potential defines $\pi_0(\mathbf{x}) = -\partial\psi = \mathbf{x}_0 - \mathbf{x}$. Then with a slight abuse of notation, we might switch on the velocity direction using a scalar $\sigma$ of the form $\sigma(s) = 1/(1+\exp(-\beta s))$ for $\beta \in\Rbb_+$ with $s = \mathbf{v}^\tr\wh{{\dot{\x}}}$ for some constant vector $\mathbf{v}$. The resulting HD0 policy becomes
\begin{align}
    \pi_0^\sigma(\mathbf{x},{\dot{\x}})
    &= \sigma(\mathbf{v}^\tr\wh{{\dot{\x}}}) \pi_0(\mathbf{x}) \\\nonumber
    &= \frac{\mathbf{x}_0 - \mathbf{x}}{1+\exp\big(-\beta \mathbf{v}^\tr\wh{{\dot{\x}}}\big)}.
\end{align}
Finally, we can choose to scale that using a Euclidean Finsler energy $\mathcal{L}_e = \frac{1}{2}\|{\dot{\x}}\|^2$ to get a final geometry of the form
\begin{align}
    \pi(\mathbf{x},{\dot{\x}}) = \frac{\frac{1}{2}\|{\dot{\x}}\|^2\big(\mathbf{x}_0 - \mathbf{x}\big)}{1+\exp\big(-\beta \mathbf{v}^\tr\wh{{\dot{\x}}}\big)}
\end{align}

Note that another way to design an HD2 geometry is by designing a Finsler geometry derived from a Finsler energy. We can use one Finsler energy to derive the geometry, and a second Finsler energy to define the priority metric. Energization would then bend the geometry associated with the second (priority) energy to align with the first energy's geometry.

A common pattern for designing Finsler energies is to start with a Riemannian (kinetic) energy in 1D of the form $\mathcal{K}(x,\dot{x}) = \frac{1}{2}g(x)\dot{x}^2$ where $g(x)$ is everywhere positive and dependent on position only. We then scale it by a scaling function that depends on velocity directionality. Since we are in 1D, that normalized velocity becomes a switch on the sign of the velocity, so we can design this scaling function in two parts
\begin{align}
    \sigma(\dot{x}) = \left\{
    \begin{array}{cc}
        \sigma_-(\dot{x}) & \mbox{for $\dot{x} \leq 0$}\\
        \sigma_+(\dot{x}) & \mbox{for $\dot{x} > 0$}
    \end{array}
    \right..
\end{align}
The energy then becomes $\mathcal{L}_e(x,\dot{x}) = \frac{1}{2}\wt{g}(x,\dot{x})\dot{x}^2$ where $\wt{g}(x,\dot{x}) = \sigma(\dot{x})g(x)$, which is effectively two different metrics depending on directionality of the velocity.

Note that in the 1D case, the derived $m(x,\dot{x}) = \partial^2_{\dot{x}\dot{x}}\wt{g}$ must have the property that
\begin{align}
\mathcal{L}_e(x,\dot{x}) = \frac{1}{2}\wt{g}(x,\dot{x})\dot{x}^2 = \frac{1}{2}m(x,\dot{x})\dot{x}^2,
\end{align}
which is only true globally when $m = \wt{g}$. This in general is not necessarily true in higher dimensions.

Finally, we describe a common pattern in designing potential functions. Potential functions define forces that push against the fabric system's mass. The mass defines the fabric priority; the stronger the mass the higher the priority the fabric takes. Thought another way, these masses define the fabric system's inertia; a given force has less effect pushing a sizable mass (high priority) than pushing a small mass (low priority). It is often important that a potential function adjust effectively to the fabric's priority. However, it can be challenging to design potential forces directly to scale in an intuitive way.

For that reason, it's often easier to design a potential function by designing its gradient in terms of a priority metric of its own and a desired {\em acceleration}. Specifically, if $\wt{\psi}(\mathbf{x})$ is a function (similar to the potential function in form) defining an acceleration policy $\pi(\mathbf{x}) = -\partial\wt{\psi}$, and if $\mathbf{M}$ is its associated priority, the resulting potential function $\psi(\mathbf{x})$ would be one for which $-\partial\psi = -\mathbf{M}\,\partial\wt{\psi}$ (force equals mass times acceleration) if such a potential function exists. Given an expression for a vector field, it's easy to check whether it is a gradient of a potential by taking its Jacobian and verifying the Jacobian is symmetric (it's a well-known result in multivariate calculus that a vector field is the gradient of a potential if and only if its Jacobian is everywhere symmetric).

One-dimensional vector fields are always symmetric so we can always design one-dimensional potentials using an acceleration policy and an associated priority mass. For higher-dimensional potentials, one can use the results discussed in \cite{cheng2018rmpflow} (Appendix D.4) that show that if the acceleration policy is radially symmetric and the metric is compatibly symmetric along the same radial lines, the resulting vector field defined by scaling the acceleration policy by the metric has symmetric Jacobian and is therefore the gradient of some potential.

When designing potential functions in this way, it is necessary only to show that the potential exists. Knowing that, we can implement the system using only its gradient given by the metric-scaled acceleration policy. The resulting system then fits the theory and we have convergence and stability guarantees.

\section{Additional Experiments: Equality Constraints}

This supplementary experiment is a specific realization of constrained fabrics. In principle, constrained fabrics can be realized in many different ways, e.g., by explicitly calculating the constraint Lagrange multiplier and adding the resulting constraint force to the unconstrained fabric as discussed in Appendix \ref{apx:LemmaConstrainedGeometricFabrics}. This method would additionally need to employ Baumgarte stabilization \cite{flores2011parametric} to combat numerical drift. Other approaches may be taken as well such as solving a nonlinear program subject to a discretized momentum-impulse form of the fabric and additional equality constraints as done by Stewart and Trinkle \cite{stewart2000implicit}. There are also many other different discretization and nonlinear program schemes that could be created to solve constrained fabrics. In this particular case, we create the following discretization,

\begin{align}
    \mathbf{q} &= \mathbf{q}_k \\
    {\dot{\q}} &= \frac{\mathbf{q}_{k+1} - \mathbf{q}_{k-1}}{2 \Delta t} \\
    {\ddot{\q}} &= \frac{\mathbf{q}_{k+1} - 2 \mathbf{q}_k + \mathbf{q}_{k-1}}{\Delta t^2}
\end{align}
where $k$ indicates the discrete time index and $\mathbf{q}_{k-1}$ and $\mathbf{q}_k$ are known. We then minimize the following cost function,

\begin{align}
    L(\mathbf{q}, {\dot{\q}}) = \frac{1}{2} ({\ddot{\q}}_d - {\ddot{\q}})^\tr \mathbf{M}(\mathbf{q}, {\dot{\q}}) ({\ddot{\q}}_d - {\ddot{\q}}) + \frac{1}{2} \lambda \Cbb(\mathbf{q}_{k+1})^2
\end{align}
where ${\ddot{\q}}_d$ is the desired unconstrained acceleration of the fabric calculated as ${\ddot{\q}}_d = -\mathbf{M}(\mathbf{q}, {\dot{\q}})^{-1} \mathbf{f}(\mathbf{q}, {\dot{\q}})$, $\Cbb(\mathbf{q}_{k+1})$ is a scalar constraint function of interest, and $\lambda$ is a weight. The first term solves for system acceleration that naturally follows the fabric while the second term is a penalty on the constraint, where $\lambda$ can be set to a large value to enforce the constraint. Altogether, this formulation implicitly solves for a new acceleration that adds a constraint force on the unconstrained fabric theoretically realizing (\ref{eqn:ConstrainedGeometricFabrics}).

We implemented the above for the geometric fabric designed in Section \ref{sec:SystemExperiments} subject to the constraint that the end-effector must lie on a plane. We repeatedly solved the program with $\Delta t = 0.04 s$ and a Gauss-Newton optimizer with four iterations per time step to forward integrate the system. The result is shown in Fig. \ref{fig:constrained_fabric}. The robot can maintain the constraint at the micrometer level and the end-effector can settle to within 1 mm of the closest point on the plane to the target.

\begin{figure}[!t]
  \centering
  \includegraphics[width=.8\linewidth]{figures/constrained_fabric.png}
  \caption{Constrained geometric fabric moving the robot end-effector to orange target points while the end-effector is constrained to the blue plane.}
  \label{fig:constrained_fabric}
\end{figure}

\section{End-effector Pose Control}
There are many different ways geometric fabrics can be used to control a robot, including end-effector pose control. This is a straight-forward extension to design outlined in Section \ref{sec:SystemExperiments} that attracted a point on the end-effector towards a target point. In this case, we calculate three additional task spaces that map to each axis of a orthonormal coordinate system placed at the end-effector. In each of these axis spaces, we place the very same attractor components discussed in Section \ref{subsubsec:attraction}. With this extension to the original fabric, we can exhibit full pose control over the end-effector.

\section{Additional experiments: Cubby navigation heuristics}
An important use of geometric fabrics is to heuristically bias the system without affecting convergence to local minima of the potential. This feature is especially useful for designing global navigation heuristics. A geometric fabric system as described in Subsection~\ref{sec:SystemExperiments} is subject to local minima since the target attractor is greedy. In many environments, though, simple heuristics can coarsely shape the behavior en route so that these local avoidance behaviors are sufficient for global navigation as well (similar to the heuristics described in \cite{ratliff2018rmps} Section V-D).

Section XI of \cite{geometricFabricsForMotionArXiv2020}  describes such an experimental setup in detail. In those experiments, geometric fabrics are used to navigate between a wall of cubbies and a set of floor boxes in a standard pick and place industrial setting. The heuristic geometries encoded simple geometric fields that push the system away from the interior of any cubby and box except for the one it wanted to enter into. The final system successfully navigated between any pair of cubby and surrounding boxes even with perturbations to the environment.

The simplicity and success of these heuristic terms suggest that learning such heuristics conditioned on environmental features over the top of an existing obstacle navigation fabric system would be straightforward, especially given automatic differentiation tools such as \cite{LiRmp2RobotLearning2021} and leveraging imitation learning methods such as \cite{xie2021imitation}. Fabrics constitute a compact encoding of behavioral components that span multiple tasks which could make it a powerful medium for learning of strongly generalizing policies. It is an area of future work to explore these learning applications in detail, especially in the context of meta learning shared nominal fabric layers.

\subsection{Multivariate calculus notation}
\label{apx:CalculusNotation}

We avoid tensor notation here whenever possible in favor of the more compact and familiar matrix-vector notation of advanced calculus.

We assume each partial derivative generates a new index ranging over the partials. The partial $\partial_\mathbf{x} f(\mathbf{x})$ of a function $f:\Rbb^n\rightarrow\Rbb$ is treated as a column vector. And if $\mathbf{g}:\Rbb^n\rightarrow\Rbb^m$ is a map, $\mathbf{g}(\mathbf{x})$ itself is treated as a column vector and $\partial_\mathbf{x} \mathbf{g}(\mathbf{x}) \in \Rbb^{m\times n}$ denotes a matrix with partials ranging over the columns. When there are multiple partials as in $\partial_{\mathbf{x}\mathbf{y}}h(\mathbf{x},\mathbf{y})$, we apply the partials in the order listed. In this case, we apply $\mathbf{x}$ first to create a column vector $\partial_\mathbf{x} h(\mathbf{x},\mathbf{y})$ then apply $\mathbf{y}$ to create a matrix $\partial_\mathbf{y}\big(\partial_\mathbf{x} h\big)$.
When more partials are involved, the resulting object becomes a higher-order tensor, but we use the convention that multiplication on the right by a matrix or vector contracts the indices in the reverse order they were generated. For instance, if $b(\mathbf{x},\mathbf{y},\mathbf{z})$ is a function, $\partial_{\mathbf{x}\mathbf{y}\mathbf{z}}b$ is a three index tensor, but $\partial_{\mathbf{x}\mathbf{y}\mathbf{z}}b \mathbf{v}$ will unambiguously contract along the $\mathbf{z}$ partial first leaving a matrix whose rows align with the partials along $\mathbf{x}$ and whose columns align with the partials along $\mathbf{y}$. Specifically, $\partial_{\mathbf{x}\mathbf{y}\mathbf{z}}b \mathbf{v} = \sum_k \frac{\partial}{\partial z^k}\big(\partial_{\mathbf{x}\mathbf{y}}b\big) v^k$. This convention is applied recursively when more than three variables are involved.

\section{Extended Details on Experiments}

\begin{table}[!htb]
	\caption{End-effector target error (ETE) and end-effector wall distance (EWD) as a function of policy class (GF or RMP) and target wall distance (TWD).}
	\label{table:targets_wall}
	\begin{center}
	    \vspace{-4mm}
		\begin{tabular}{| p{1cm} | c | c | c | c |}
		    \hline
			& ETE (m) & ETE (m) & EWD (m)  & EWD (m)\\
			\hline
			TWD (m)
			& GF & RMP & GF & RMP\\
			\hline
			0.3 & 3.041e-5 & 8.09e-05 & 3.0e-1 & 3.001e-1 \\
			\hline
			0.275 & 2.455e-5 & 7.89e-05 & 2.75e-1 & 2.751e-1 \\
			\hline
			0.25 & 1.919e-5 & 7.714e-05 & 2.5e-1 & 2.501e-1 \\
			\hline
			0.225 & 1.432e-5 & 7.568e-05 & 2.25e-1 & 2.251e-1 \\
			\hline
            0.2 & 1.151e-5 & 7.455e-05 & 2.0e-1 & 2.001e-1 \\
			\hline
			0.175 & 1.549e-5 & 7.389e-05 & 1.75e-1 & 1.751e-1 \\
			\hline
			0.15 & 2.88e-5 & 7.564e-05 & 1.5e-1 & 1.501e-1 \\
			\hline
			0.125 & 5.834e-5 & 1.152e-4 & 1.251e-1 & 1.251e-1\\
			\hline
			0.1 & 1.564e-4 & 8.482e-4 & 1.002e-1 & 1.009e-1 \\
			\hline
			0.075 & 1.733e-3 & 8.84e-3 & 7.673e-2 & 8.38e-2 \\
			\hline
			0.05 & 1.939e-2 & 1.677e-2 & 6.939e-2 & 6.672e-2 \\
			\hline
			0.025 & 4.315e-2 & 4.351e-2 & 6.816e-2 & 6.846e-2 \\
			\hline
			0.0 & 6.761e-2 & 9.625e-2 & 6.761e-2 & 9.62e-2 \\
			\hline
			-0.025 & 9.225e-2 & 1.21e-1 & 6.725e-2 & 9.602e-2 \\
		    \hline
			-0.05 & 1.17e-1 & 1.459e-1 & 6.698e-2 & 9.595e-2 \\
			\hline
			-0.075 & 1.418e-1 & 1.715e-1 & 6.683e-2 & 9.651e-2 \\
			\hline
			-0.1 & 1.668e-1 & 1.957e-1 & 6.681e-2 & 9.571e-2 \\
			\hline
			-0.125 & 1.919e-1 & 2.214e-1 & 6.691e-2 & 9.644e-2 \\
			\hline
			-0.15 & 2.171e-1 & 2.457e-1 & 6.709e-2 & 9.581e-2 \\
			\hline
			-0.175 & 2.422e-1 & 2.714e-1 & 6.724e-2 & 9.647e-2 \\
			\hline
			-0.2 & 2.673e-1 & 2.957e-1 & 6.732e-2 & 9.578e-2 \\
			\hline
			-0.225 & 2.924e-1 & 3.213e-1 & 6.735e-2 & 9.638e-2 \\
			\hline
			-0.25 & 3.173e-1 & 3.455e-1 & 6.737e-2 & 9.564e-2 \\
			\hline
			-0.275 & 3.424e-1 & 3.712e-1 & 6.737e-2 & 9.636e-2 \\
			\hline
			-0.3 & 3.674e-1 & 3.961e-1 & 6.738e-2 & 9.632e-2 \\
			\hline
		\end{tabular}
	\end{center}
	\vspace{-4mm}
\end{table}

\subsection{Derivations on manifolds}
\label{apx:Manifolds}

Nonlinear geometry is most commonly developed in the mathematics literature in terms of smooth manifolds, and often in an abstract coordinate-free notation which can be challenging for practitioners. Even tensor notation can be a barrier to entry for readers unfamiliar with the algebra. To make our presentation more accessible, we will stick to coordinate descriptions and simpler vector notations from advanced calculus (see Appendix~\ref{apx:CalculusNotation}) similar to many books on classical mechanics aimed at an engineering audience.

Manifolds will be formally defined in the traditional way (see \cite{LeeSmoothManifolds2012} for a complete introduction), but for our purposes, we will consider them simply as $d$-dimensional spaces $\mathcal{Q}$ with elements identified with $\mathbf{q}\in\Rbb^d$ in $d$ generalized coordinates, and when we say $\mathbf{q}\in\mathcal{Q}$ we mean $\mathbf{q}$ is the coordinate representation of the manifold point described in our chosen coordinates. Often we implicitly assume a system evolves over time $t$ in a trajectory $\mathbf{q}(t)$ with velocity ${\dot{\q}} = \frac{d\mathbf{q}}{dt}$, and say that coordinate velocity vector ${\dot{\q}}$ is an element of the tangent space $\mathcal{T}_\mathbf{x}\mathcal{Q}$, which can be thought of as the space of all velocity vectors that can be generated in that way.
When discussing general tangent space vectors (elements of the tangent space without reference to a specific trajectory), we often use the notation $\mathbf{v}\in\mathcal{T}_\mathbf{q}\mathcal{Q}$.

\subsection{A note on covariance and transform trees}
\label{apx:NoteOnCovariance}

Geometric consistency across changes of coordinates on a manifold is an important concept since many manifolds cannot be described globally by a single set of coordinates. Algorithms and systems should not change their fundamental behavior on the manifold based on the coordinate system they're described in. When an algorithm or system has such a coordinate independence, it's called {\em covariant}.

It is well known that the Euler-Lagrange equation is covariant to changes of coordinate \cite{LeonardSusskindClassicalMechanics}, so anything derived here from the Euler-Lagrange equation is automatically covariant. For everything else, such as HD2 geometries and bending terms, we are not rigorous with developing specific intrinsically covariant subclasses. Such covariant subclasses exist and can be derived, but the mathematics behind that can be complicated and unnecessarily restrictive. In practice, we instead use transform trees to derive covariant systems.

A transform tree is simply a tree whose nodes are different spaces and whose edges are the maps linking those spaces. For instance, the simplest tree is a single branch $\phi:\mathcal{Q}\rightarrow\mathcal{X}$ mapping generalized coordinates $\mathcal{Q}\subset\Rbb^d$ to a single space $\mathcal{X}\subset\Rbb^n$. In fact, all trees can be expressed in this form since paths from the root to any given node are unique and the maps encountered along that path can be combined into a single composed map.
The single branch tree is constructed by stacking all nodes into a single space and defining a map by stacking these unique path-composed maps into a single map. In practice, more expressive tree structures are important for computational performance, but conceptually it suffices to think of simply a single link $\phi$ as described here. (Note that $n\geq d$ when the tree is well-formed so that policies on the nodes can fully constrain the space.)

Equation~\ref{eqn:PullbackCommutesWithEL} shows that the pullback operation defined on a transform tree is compatible with the pullback implicit in the covariant Euler-Lagrange equation. This is because pullback is compatible with the definition of covariance. Calculationally, we prove that a system is covariant if, under a change of coordinates, the resulting system is simply scaled by the Jacobian (transpose) of the change of coordinates map. For instance, if $\phi$ denotes a change of coordinates (i.e. when $n = d$ and $\phi$ is smooth and invertible making it what is known as a {\em diffeomorphism}), a term $\mathbf{h}(\mathbf{x},{\dot{\x}}, {\ddot{\x}})$ is covariance if its expression in $\mathcal{Q}$ amounts to $\wt{\mathbf{h}}(\mathbf{q},{\dot{\q}},{\ddot{\q}}) = \mathbf{J}^\tr \mathbf{h}(\mathbf{x},{\dot{\x}},{\ddot{\x}})$ where $\mathbf{J} = \partial\phi$ is the Jacobian of $\phi$ and ${\dot{\J}}$ is its time derivative (and $\mathbf{x} = \phi(\mathbf{q}), {\dot{\x}} = \mathbf{J}{\dot{\q}}, {\ddot{\x}} = \mathbf{J}{\ddot{\q}} + {\dot{\J}}{\dot{\q}}$). The utility of this definition of covariance is that if $\wt{\mathbf{h}}(\mathbf{q},{\dot{\q}}, {\ddot{\q}}) = \mathbf{0}$ then in $\mathcal{Q}$ we get the same
\begin{align}
    &\wt{\mathbf{h}}(\mathbf{q},{\dot{\q}},{\ddot{\q}}) =\mathbf{0} \\\nonumber
    &\ \ \Rightarrow  \mathbf{J}^\tr\mathbf{h}(\mathbf{x},{\dot{\x}},{\ddot{\x}}) = \mathbf{0} \\\nonumber
    &\ \ \Rightarrow \mathbf{h}(\mathbf{x},{\dot{\x}},{\ddot{\x}}) = \mathbf{0}
\end{align}
since $\mathbf{J}$ is invertible when $\phi$ is a change of coordinates. In other words, covariant equations hold in all coordinate systems.

Now, looking at Equation~\ref{eqn:SpecPullback} and its derivation leading to Equation~\ref{eqn:PullbackDerived}, we see that when $\phi$ is a change of coordinates, the system transforms covariantly by definition. Specifically, when the term is $\mathbf{h}(\mathbf{x},{\dot{\x}},{\ddot{\x}}) = \mathbf{M}(\mathbf{x},{\dot{\x}}){\ddot{\x}} + \mathbf{f}(\mathbf{x},{\dot{\x}})$, the pullback operation computes $(\mathbf{J}^\tr\mathbf{M}\mathbf{J}\big){\ddot{\q}} + \big(\mathbf{f} + \mathbf{M}{\dot{\J}}{\dot{\q}}\big) = \mathbf{J}^\tr\big(\mathbf{M}{\ddot{\x}} + \mathbf{f}\big) = \mathbf{J}\,\mathbf{h}(\mathbf{x},{\dot{\x}},{\ddot{\x}})$ (we're just reversing the argument used to derive it).
Therefore, if we have any arbitrary transform tree $\phi:\mathcal{Q}\rightarrow\mathcal{X}$ (viewed as a single branch) and we want to change the coordinates using $\zeta:\wt{\mathcal{Q}}\rightarrow\mathcal{Q}$, we can simply build a transform tree linking these transforms
\begin{align}
    \wt{\mathcal{Q}}\stackrel{\zeta}{\longrightarrow} \mathcal{Q}
    \stackrel{\phi}{\longrightarrow}\mathcal{X}
\end{align}
and use standard pullback operations to get covariance automatically. For instance, if $\mathcal{Q}$ is the generalized coordinates of the robot, we can define a transform tree with $\mathcal{Q}$ at the root and a collection of task spaces extending from it such as the end-effector space, body points, distances to obstacles, etc. Then, if we change to a different coordinate representation of the generalized coordinates we simply express that as part of the transform tree and the fundamental behavior of our system does not change because of the natural covariance of the transform tree itself. An important example is when the generalized coordinates includes a free flying SE3 base such as in a drone with an SO3 rotation component. SO3 cannot be globally expressed in a single coordinate system, so we inevitably need to change coordinates at times. (We can also change representations, between Euler angles, quaternions, axis-angle and matrix representations, as long as the transform is represented as part of the tree.) Covariance ensures that the system behavior remains consistent under these representational transformations.

Another way to understand the transform tree is that it represents an optimization problem. Pullback can equivalently be viewed as the solution to a least squares optimization problem
\begin{align}
    &\min_{\ddot{\q}} \frac{1}{2}\|{\ddot{\x}}^d - {\ddot{\x}}\|_\mathbf{M}^2 \\\nonumber
    &\mathrm{s.t.}\ {\ddot{\x}} = \mathbf{J}{\ddot{\q}} + {\dot{\J}}{\dot{\q}},
\end{align}
analogous to the least squares formulation of classical mechanics under Gauss's Principle of Least Constraint \cite{udwadia1996analytical}. Because it's an optimality principle, the fundamental minimum does not change under changes of coordinates.

redpajama_set_name: RedPajamaArXiv