This paper was converted on www.awesomepapers.org from LaTeX by an anonymous user.
Want to know more? Visit the Converter page.

A variational approach to geometric mechanics
for undulating robotic locomotion

Sean Even1, Patrick S. Martinez2, Cora Keogh3, Oliver Gross4,∗, Yasemin Özkan-Aydın1,∗ and Peter Schröder2 1University of Notre Dame, 2California Institute of Technology, 3Dublin City University, 4EPFL (Swiss Federal Institute of Technology, Lausanne)
∗Corresponding authors.
Abstract

Limbless organisms of all sizes use undulating patterns of self-deformation to locomote. Geometric mechanics, which maps deformations to motions, provides a powerful framework to formalize and investigate the theoretical properties and limitations of such modes of locomotion. However, the inherent level of abstraction poses a challenge when bridging the gap between theory or simulations and laboratory experiments. We investigate the challenges of modeling motion trajectories of an undulating robotic locomotor by comparing experiments and simulations performed with a variational integrator. Despite the extensive simplifications that the model based on a geometric variation principle entails, the simulations show good agreement on average. Notably, our approach merely requires the knowledge of the dissipation metric—a Riemannian metric on the configuration space, which can in practice be approximated by means closely resembling resistive force theory.

Index Terms:
Geometric mechanics, biomechanics, bioinspired robotics, undulating locomotion, limbless locomotion, kinematics

I Introduction

Understanding how animals navigate and propel themselves through diverse environments is a fundamental question in biomechanics [1]. The dynamics of locomotion are intertwined with environmental factors—complex interactions between an organism’s morphology and the surrounding medium critically influence performance metrics such as velocity, maneuverability, and energy expenditure. By studying these interactions, we gain valuable insights into how animals have evolved to optimize their movement, which can, in turn, inspire innovative designs in robotics and other fields seeking to mimic or leverage these natural strategies.

Among all organisms, snakes stand out for their exceptional versatility in locomotion. Navigating diverse environments—from dense underbrush and rocky crevices to smooth, sandy surfaces—with remarkable adaptability [2]. They employ various locomotion strategies, including undulating, sidewinding, and concertina movements [3, 4] with each gait finely tuned to overcome specific environmental challenges. This adaptability not only allows them to efficiently exploit a wide range of ecological niches but also highlights the sophistication of their biomechanical systems, which contribute to their widespread success as a species across varied habitats [5, 6, 7, 8]. Their unique locomotion strategies have inspired numerous robotic studies, aiming to replicate these natural mechanisms for advanced robotic mobility[9, 10, 11, 12, 13, 14].

Robotic snake locomotion has been heavily influenced by the application of geometric mechanics, which provides a unified framework for controlling these dynamic systems [15, 16, 17, 18, 19, 20]. One of the foundational models, proposed by Hirose [21], utilizes a serpenoid shape to represent the continuous curve of a snake’s body, allowing for precise control of its undulatory motion. This geometric approach has been extended by Branyan et al. [22], who implemented a gait optimization strategy to enhance locomotion efficiency, through coordinated shape changes. Similarly, [23] applied soap–bubble optimization to identify optimal gaits based on the curvature of the shape space, while [24] optimized gaits for their power consumption based on a power dissipation metric. Expanding these principles to more complex terrains, Dai et al. [25] optimized robotic snake gaits on granular surfaces, leveraging a geometric framework to simplify the dynamic interaction between the robot and the shifting terrain. These control approaches demonstrate the flexibility and effectiveness of geometric mechanics in enabling precise, efficient locomotion across various environments.

The inherent level of abstraction, however, can pose a challenge when geometric mechanics-based simulations are to be adapted to reproduce or predict real-world experiments.

Refer to caption
Figure 1: Geometric mechanics maps gaits of undulatory robotic locomotors to motion trajectories (blue) in world coordinates. The accuracy of this map compared to laboratory experiments (green) depends on the choice of model parameters, for which the Riemannian metric on the configuration space provides a natural description.

This paper addresses the challenges of bridging the gap between abstract theoretical models and empirical validation using a practical example. Restricting our attention to the case of undulating locomotion, we

  • •

    calibrate a general simulation framework [26] to match motion trajectories of a snake-like robot.

  • •

    employ an approximation of the dissipation metric which closely resembles Resistive Force Theory (RFT).

  • •

    employ a variational integrator, derived from a geometric variational principle, to integrate motion trajectories from given shape sequences.

  • •

    compare the locomotive performance of simulations and laboratory experiments for various gaits. In particular, we validate if performance differences of the gaits carry over to the real world.

This allows us to test whether a highly simplified model can still accurately describe the system’s behavior, aiming to replace costly detailed simulations with more efficient calculations.

A convenient calibration of the geometric mechanics model for a specific laboratory setup enables a more efficient approximation of the trajectories. In particular, by demonstrating the effectiveness of an intuitive approach to geometric mechanics that only requires the specification of a Riemannian metric in the form of RFT-like calibrations, we show the practical accessibility to the multitude of benefits this abstract setup provides. Ultimately, this paves the way to new areas of application such as inverse design or reinforcement learning [27].

II Preliminaries

Geometric mechanics analyzes dynamical systems based on a special structure of the configuration space ℳ{\mathcal{M}}. For any state of a locomotive system in ℝ3\mathbb{R}^{3} the locomotors shape γ^∈𝒮{\hat{{\gamma}}}\in{\mathcal{S}} is treated distinguished from its respective positioning γ∈ℳ{\gamma}\in{\mathcal{M}}. Two positioned shapes have the same shape if they differ only by a rigid body motion g∈SE⁡(3)g\in\operatorname{SE}(3) in world space111Although the locomotion that we consider in the present paper is limited to two dimensions, we treat it as a special case of the three-dimensional case, which allows us to establish a more general set of equations that was used for the numerical implementation.. That is, g:ℝ3→ℝ3,x↦A​x+bg\colon\mathbb{R}^{3}\to\mathbb{R}^{3},\,x\mapsto Ax+b for some A∈SO⁡(3)A\in\operatorname{SO}(3) and b∈ℝ3b\in\mathbb{R}^{3}.

As a result, the configuration space ℳ{\mathcal{M}} decomposes into six-dimensional fibers,

{γ=g​(γ^)∣g∈SE⁡(3)}⊂ℳ,\{\gamma=g({\hat{{\gamma}}})\mid g\in\operatorname{SE}(3)\}\subset{\mathcal{M}}, (1)

consisting of all possible positions of a body posture γ^∈𝒮{\hat{{\gamma}}}\in{\mathcal{S}}.

Therefore, the dynamics of an undulating robotic locomotor can be described by a smooth curve γ^:[0,T]→𝒮{\hat{{\gamma}}}\colon[0,T]\to{\mathcal{S}} into the shape space 𝒮{\mathcal{S}} together with a smooth map g:[0,T]→SE⁡(3)g\colon[0,T]\to\operatorname{SE}(3) positioning each shape in world space. A path γ:[0,T]→ℳ{\gamma}\colon[0,T]\to{\mathcal{M}} in the configuration space is said to be a lift of γ^{\hat{{\gamma}}} if π​(γ)=γ^\pi({\gamma})={\hat{{\gamma}}}, where π:ℳ→𝒮\pi\colon{\mathcal{M}}\to{\mathcal{S}} denotes the projection from the configuration space onto the shape space. By decoupling shapes from their position, optimization and control problems in robotics can often be reduced to a problem on the shape space 𝒮{\mathcal{S}}.

II-A Geometric locomotion

For scenarios involving, e.g., highly damped environments or granular media one can exploit a linear relationship between a system’s shape changes γ^′\hat{\gamma}^{\prime} and its position to integrate the motion trajectory of a shape changing body [16, 23]. Typically, a local connection form ϖ​(γ^)\varpi(\hat{\gamma}) is used to map infinitesimal shape changes to body velocities in world coordinates by

γ′=−ϖ​(γ^)​γ^′.\displaystyle\gamma^{\prime}=-\varpi(\hat{\gamma})\,\hat{\gamma}^{\prime}. (2)

Eq. (2) assures that any 11-parameter family t↦γ^t∈𝒮t\mapsto{\hat{{\gamma}}}_{t}\in{\mathcal{S}} of body poses can, up to a global rigid body motion, be uniquely lifted to a 11-parameter family t↦γt∈ℳt\mapsto{\gamma}_{t}\in{\mathcal{M}} describing the motion trajectory of the shape changing body.

When we consider a closed loop in the shape space 𝒮{\mathcal{S}}, it describes a periodic sequence of shapes, which is referred to as a gait. Despite the periodic nature of gaits, the resulting lift γ{\gamma} does in general not close up (Fig. 1). This aperiodicity is precisely why the geometric description of the dynamical system can lead to a net displacement of shape-changing bodies. The resulting net displacement of the locomotor after one gait cycle is given by g​(0)−1​g​(T)∈SE⁡(n)g(0)^{-1}g(T)\in\operatorname{SE}(n) (Fig. 1) and the extent of this displacement depends on the geometry of the fiber bundle [28].

II-B A variational approach

Instead of numerically integrating Eq. (2), variational integrators provide an alternative approach to integrating motion trajectories from shape sequences. They are derived following the guiding principle to first discretize and then optimize [29]. By construction, they exhibit several advantageous properties as they are automatically symplectic and momentum preserving, and exhibit good energy behavior for exponentially long times [29, 30].

A unified framework that makes variational integrators straightforward to use for a variety of scenarios—including the one we examine—has been proposed by Gross et al. [26]. In this section, we will briefly outline the key aspects of this approach, while for a more detailed exposition and derivations, we refer to the original reference. Notably, this approach does not require explicit prior knowledge of the local connection forms. Instead, it merely requires knowledge of the dissipation metric on ℳ{\mathcal{M}} for which even a rough approximation in the spirit of RFT [31, 32] has proven sufficient in practice.

II-B1 Discretization and local dissipation metric

We discretize the snake-like undulating locomotor at a given time tt as a polygonal curve γt{\gamma}^{t} consisting of NN vertices in ℝ3\mathbb{R}^{3} (Fig. 2). Moreover, we restrict our attention to time-discrete sequences of discrete positioned shapes, i.e., an indexed sequence γ0,…,γT∈ℳ=(ℝ3)N.{\gamma}^{0},\ldots,{\gamma}^{T}\in{\mathcal{M}}=(\mathbb{R}^{3})^{N}. To each vertex γkt{\gamma}^{t}_{k} of a curve γt{\gamma}^{t} we assign a unit tangent vector Tkt∈S2T^{t}_{k}\in S^{2}. When a positioned shape is transformed by a rigid motion γ↦A​γ+b{\gamma}\mapsto A{\gamma}+b, the tangent vectors at the vertices are transformed as well by T↦A​TT\mapsto AT.

Refer to caption
Figure 2: We represent the robotic system comprised of linked elements as a polygonal curve. Displacements of vertices anisotropically dissipate energy to the environment which is represented by local dissipation tensors (see Sec. II-B)

Displacing a body in a viscous medium causes energy dissipation due to, e.g., viscous drag. In the spirit of RFT, we approximate this energy dissipation by summing up the energy dissipated from displacing individual vertices, thus neglecting any effects from interactions. Treating vertices as local dissipation elements (Fig. 2) we associate them with anisotropic local dissipation tensors, i.e., symmetric and positive definite blocks of the form

Dkt≔wk​(I+(ϵ−1)​Tkt⊗Tkt)∈ℝ3,3.\displaystyle D_{k}^{t}\coloneqq w_{k}(I+(\epsilon-1)T^{t}_{k}\otimes T^{t}_{k})\in\mathbb{R}^{3,3}. (3)

Here, wkw_{k} are integration weights, while the anisotropy ratio ϵ∈(0,1]\epsilon\in(0,1] controls the ease of tangential motion compared to normal motion (Fig. 2). For ϵ≈0\epsilon\approx 0 tangential displacements cause close to no energy dissipation, while for ϵ=1\epsilon=1 the tensor becomes isotropic and there are no preferred directions and no net displacement can be achieved [31].

II-B2 Variational energy

The discrete variational energy we consider is of the form

ℰ​(γ0,…,γT)=12​∑t=0T−1⟨𝒟(t,t+1)​Δ​p(t,t+1),Δ​p(t,t+1)⟩.\displaystyle{\mathcal{E}}({\gamma}^{0},\ldots,{\gamma}^{T})=\tfrac{1}{2}\sum_{t=0}^{T-1}\langle{\mathcal{D}}^{(t,t+1)}\Delta p^{(t,t+1)},\Delta p^{(t,t+1)}\rangle. (4)

Here, denoting the concatenation of the vertex positions of γt{\gamma}^{t} by ptp^{t}, Δ​p(t,t+1)≔pt+1−pt\Delta p^{(t,t+1)}\coloneqq p^{t+1}-p^{t} and 𝒟(t,t+1){\mathcal{D}}^{(t,t+1)} is symmetric, positive definite and models the dissipation metric. We follow the definition in [26] and choose 𝒟(t,t+1)≔12​(𝒟t+𝒟t+1){\mathcal{D}}^{(t,t+1)}\coloneqq\tfrac{1}{2}({\mathcal{D}}^{t}+{\mathcal{D}}^{t+1}) where 𝒟t∈ℝ3​N,3​N{\mathcal{D}}^{t}\in\mathbb{R}^{3N,3N} are block-diagonal matrices with blocks of the form given in Eq. (3). This energy measures the total energy dissipation affected by the displacements Δ​p(t,t+1)\Delta p^{(t,t+1)}.

II-B3 Variational integrator

A physical motion extremizes Eq. (4) under variations by rigid body transformations [26]. The resulting Euler-Lagrange equation corresponds to the non-holonomic constraints of the system [29]. For bodies initially at rest, they are given by

μ​(γt,γt+1)=0.\mu({\gamma}^{t},{\gamma}^{t+1})=0. (5)

Here, μ\mu is the geometric momentum, which is defined as

μ​(γt,γt+1)≔−12​∑k(μrot​(γt,γt+1)kμtran​(γt,γt+1)k),\displaystyle\mu({\gamma}^{t},{\gamma}^{t+1})\coloneqq-\tfrac{1}{2}\!\sum_{k}\begin{pmatrix}\mu_{\rm rot}({\gamma}^{t},{\gamma}^{t+1})_{k}\\ \mu_{\rm tran}({\gamma}^{t},{\gamma}^{t+1})_{k}\end{pmatrix}, (6)

where

μrot​(γt,γt+1)k≔\displaystyle\mu_{\rm rot}({\gamma}^{t},{\gamma}^{t+1})_{k}\coloneqq pkt+1×(Dt​Δ​p(t,t+1))k\displaystyle\ p^{t+1}_{k}\times(D^{t}\Delta p^{(t,t+1)})_{k}
+pkt×(Dt+1​Δ​p(t,t+1))k\displaystyle+p^{t}_{k}\times(D^{t+1}\Delta p^{(t,t+1)})_{k}
μtran​(γt,γt+1)k≔\displaystyle\mu_{\rm tran}({\gamma}^{t},{\gamma}^{t+1})_{k}\coloneqq (D(t,t+1)​Δ​p(t,t+1))k.\displaystyle\ (D^{(t,t+1)}\Delta p^{(t,t+1)})_{k}.

Since the shapes are considered to be given, for every timestep, the six constraints in Eq. (5) match the six degrees of freedom for their respective position (Eq. (1)). Therefore, Eq. (5) can be solved numerically exactly at each time step (Alg. 1), thus avoiding the accumulation of integration errors. In practice, given two consecutive shapes together with their respective material parameters, Eq. (5) is satisfied for roots of

SE⁡(3)→ℝ6,gt↦μ​(γt,gt​(γ^t+1)).\operatorname{SE}(3)\to\mathbb{R}^{6},\ g_{t}\mapsto\mu({\gamma}^{t},g_{t}({\hat{{\gamma}}}^{t+1})). (7)

That is, γ^t+1{\hat{{\gamma}}}_{t+1} is positioned by gt∈SE⁡(3)g_{t}\in\operatorname{SE}(3) such that Eq. (5) is satisfied.

Algorithm 1 — IntegrateMotionTrajectory [26]
1:shapes (γ^0,…,γ^T)({\hat{{\gamma}}}^{0},\ldots,{\hat{{\gamma}}}^{T}), parameters (w0,…,wN,ϵ)(w_{0},\ldots,w_{N},\epsilon)
2:positioned shapes γ0,…,γT{\gamma}^{0},\ldots,{\gamma}^{T}
3:γ0←γ^0{\gamma}^{0}\leftarrow{\hat{{\gamma}}}^{0}
4:for t=1,…,Tt=1,\ldots,T do
5:     gt←g_{t}\leftarrow solve μ​(γt−1,gt​(γ^t))=0\mu({\gamma}^{t-1},g_{t}({\hat{{\gamma}}}^{t}))=0 ▷\triangleright Eqs. (5, 6, 7)
6:     γt←gt​(γ^t){\gamma}^{t}\leftarrow g_{t}({\hat{{\gamma}}}^{t})
7:end for

For three spacial dimensions, determining such a positioning amounts to solving a non-linear system in six variables (Eq. (7)), while in two spacial dimensions only three variables are left.

III Shape space and gait design

We describe the examined gaits with help of the space 𝒮serp{{\mathcal{S}}_{\rm serp}} of serpenoid curves, which is commonly used to describe undulations of non-anthropomorphic organisms in low-Reynolds number environments [21, 20] and highly damped environments [20, 33].

Serpenoid curves are composed from a superposition of a standing and a traveling wave. Typically, planar curves are described in terms of their curvature function

κ​(s,t)=w1​(t)​sin⁡(2​π​ξ​s)+w2​(t)​cos⁡(2​π​ξ​s),\kappa(s,t)=w_{1}(t)\sin(2\pi\xi s)+w_{2}(t)\cos(2\pi\xi s), (8)

which, by the fundamental theorem of plane curves [34], uniquely determines the curves’ shape up to a rigid body transformation. Here, w1​(t)w_{1}(t) and w2​(t)w_{2}(t) are time-dependent coefficients that correspond to the coordinates in the shape space, and ξ\xi is the spatial frequency of body undulation (Fig. 3). Therefore, the serpenoid shape space 𝒮serp{{\mathcal{S}}_{\rm serp}} is two-dimensional and can be identified with ℝ2\mathbb{R}^{2} (Figs. 1, 3).

Notably, in contrast to other common choices of shape spaces—such as the space of actuation angles—the dimension of the shape space is independent of the number of actuators on the robot. This is favorable since it lowers the complexity for eventual gait optimization problems on the shape space.

Refer to caption
Figure 3: The ellipses corresponding to two prototypical gaits within the shape space, accompanied by the shapes of the sequence and motion trajectory.

III-A Gaits

Rieser et al. [20] find that the gaits of various undulating living systems tend to follow approximately circular, closed-loop trajectories. Therefore, for the gait design we restrict ourselves to a low-dimensional representation and consider gaits

[0,1]∋t↦(w1​(t)w2​(t))≔Rθ​(a​(cos⁡(2​π​t)σ​sin⁡(2​π​t)))+(xcyc)[0,1]\ni t\mapsto\begin{pmatrix}w_{1}(t)\\ w_{2}(t)\end{pmatrix}\coloneqq R_{\theta}\left(a\begin{pmatrix}\cos(2\pi t)\\ \sigma\sin(2\pi t)\end{pmatrix}\!\right)+\begin{pmatrix}x_{\rm c}\\ y_{\rm c}\end{pmatrix}

determined by ellipses embedded in 𝒮serp{{\mathcal{S}}_{\rm serp}} (Fig. 3). Here, RθR_{\theta} denotes a rotation in ℝ2\mathbb{R}^{2} by the angle θ\theta, (xc,yc)T(x_{\rm c},y_{\rm c})^{T} the center of the ellipse, aa the length of the larger principal axis and σ∈[0,1]\sigma\in[0,1] the ellipse’s flatness. Moreover, we consider the spatial frequency ξ\xi as a variable. In practice, we discretize those gait ellipses as polygonal curves with their vertices corresponding to time steps (Fig. 3).

For the purposes of this paper, we consider two kinds of gaits: We either draw uniform samples from the 66-dimensional parameter domain within “reasonable” bounds, or we seek for more “optimal” gaits. For the latter, we minimize a naive loss function

(σ,xc,yc,θ,a,ξ)↦−ΔCoM,(\sigma,x_{\rm c},y_{\rm c},\theta,a,\xi)\mapsto-\Delta_{\mathrm{CoM}},

where ΔCoM\Delta_{\mathrm{CoM}} denotes the magnitude of the translational part of the net displacement of the center of mass resulting from a forward simulation of one full gait cycle with Alg. 1. A prototypical gait pair is shown before and after optimization for an experiment, simulation, and the elliptical gait profile (SI Movie 0:31-0:44).

We investigate the accuracy of the framework underlying Alg. 1 to describe the motion trajectory of our robotic locomotor in three different ways. First, we qualitatively compare the resulting motion trajectories (Sec. V, Fig. 7). Second, we examine how accurately performance differences are reproduced (Sec. V-A, Fig. 8). Last, we investigate the possibility of regularizing objective functions to optimize, e.g., the power consumption of gaits (Sec. V-B, Fig. 9).

IV Experimental protocol

IV-A Snake robot

For the experimental setup, we developed a ten-degree-of-freedom planar snake robot (Fig. 4). The robot (length=0.92​m\text{length}=0.92\,\text{m}, mass=1.38​kg\text{mass}=1.38\,\text{kg}) consists of a chain of Dynamixel XL430-W250-T servo motors connected by 3D-printed ABS plastic connectors printed on a Stratasys F170 3D printer, a Robotis Open-RB 150 control board, and two 11.111.1\,V 15001500\,mAh 120C Drone batteries attached to the head and tail. This type of robot design is well-studied and intentionally simple to avoid introducing unknown sources of error associated with more complex designs.

Refer to caption
Figure 4: Description of our snake robot which includes ten planar actuators (Dynamixel XL430-W250-T) and rubber wheels that introduce frictional anisotropy. The robot also contains reflective markers at the rotational axis of the motors (marked in red) as well as on the head and tail of the robot which allows each segment to be tracked in real time using the Optitrack Motion Capture system.

To account for the frictional anisotropy of the scales, which allow snakes to slither, we attached a pair of passive rubber Lego wheels (diameter=56​mm\text{diameter}=56\,\text{mm}, thickness=11.9​mm\text{thickness}=11.9\,\text{mm}, connected by an axle) to each of the connector parts linking the motors (18.2318.23\,mm from the servo’s rotational axis) [9]. Additional wheels were placed beneath the center of mass of both the head and tail of the system (a total of 12 wheel pairs). Functionally, displacements of the wheels in the directions perpendicular to the rolling direction of the wheels generate significantly higher friction than in the tangential directions along the rolling direction.

To investigate whether the underlying geometric model can also be employed to enhance the robot’s gait efficiency, we have integrated a power sensor (Adafruit INA260) into the robot. By connecting it in series with the onboard battery, the sensor logs power consumption during operation, enabling us to calculate the Cost of Transport and quantitatively assess efficiency (see Sec. V-B).

IV-B Experimental setup

We utilized an Optitrack Motion Capture System comprised of eight PrimeX-13 cameras (120​frames/sec120\,\nicefrac{{\rm frames}}{{\rm sec}} sample rate) to track the trajectory of the reflective markers that are attached to the robot at the midpoint of the head and tail segments and the axis of rotation of each of the servo motors. Experiments are conducted in a 3​m×3​m3\,\text{m}\times 3\,\text{m} arena covered with foam mats.

At the start of each experiment, the robot was placed in the arena with its motors configured to initiate the first configuration of its gait. The robot is equipped with a power located on its head which is used to initiate trials. The simulated gait cycle was discretized and input into the system as a 50×1050\times 10 matrix. Each gait cycle was executed iteratively three times, with a 22-second pause between cycles to distinguish the multiple cycles within a single trial. We conducted three experiments for each gait cycle and calculated the mean and standard deviation (STD) of the results to assess the variability and performance consistency across trials (Fig. 5).

Refer to caption
Figure 5: Overlay of three lab trials given the same input data which correspond to the experiment in the center of Fig. 7.

IV-C Fitting of material parameters

Traditional RFT states that motion is predicated on velocity and tangential drag on individual body elements [31, 32]. Our simplified robotic segments can be modeled as point-like contacts on a chosen surface. That is, the tracking points and hence centers–of–mass are at the joints and pivot over the wheelbase (Fig. 2). We can then use actual weight measurements of the robotic components to determine wjw_{j} (Eq. (3)) and fit the drag–like anisotropy ratio ϵ\epsilon as a proxy for tangential drag. The coefficient is inversely proportional to displacement as seen in Fig. 6, meaning we can employ, e.g., a binary search [35] to find the best–RMS fit between the simulated and laboratory data. For our robot, the average of the fitted values was ϵavg=0.1865\epsilon_{\rm avg}=0.1865, which we used for all the experiments presented.

Refer to caption
Figure 6: Displacement (body lengths) vs. time graphs for various anisotropy ratios, ϵ\epsilon. Red surface shading corresponds with smaller root–mean–square (RMS) error between the simulated and lab data (blue). The yellow-highlighted plot shows the simulation using the best–fit ϵ\epsilon.

V Evaluation

To frame how close our highly simplified model is to reality, we generated twelve random gaits and subsequently optimized them according to Sec. III-A. We refer to this set of 24 gaits as the Simulation (Sim) data.

Additionally, when simulating the Sim-data set in laboratory experiments, we collect the corresponding Experimental (Exp) data with the help of the motion capture system. A third dataset is obtained by resimulating (Resim) the data of the experimental data set. This is because, at each time step, the true angles of the gait were passed to the system as its goal configuration. However, the motors did not always reach their exact goal configuration. We must therefore consider the slightly different shape sequences that have been realized in practice.

We remark that there is a bijection between the corresponding gaits of each dataset—the simulated gait produces the experimental, which in turn produces the resimulated data.

Fig. 7 shows that our highly simplified model fails to reliably reproduce individual trials. Nonetheless, it effectively captures the overall nature of the motion trajectories observed during our laboratory experiments. This is further demonstrated in the SI Movie where all three data classes are shown visually for the gait that achieves maximum CoM displacement (SI Movie 0:03-0:17).

Refer to caption
Figure 7: Motion trajectories of ten laboratory experiments are shown in green, with simulations corresponding to motion-captured data in blue. The dashed red lines indicate the respective center of mass trajectories.

V-A Mapping performance differences

To examine the Sim2Real correspondence of the proposed framework, we analyze how accurately differences in the performances of gaits are mapped. In theory, gaits that excel or underperform in simulations should exhibit similar performance trends in real-world experiments. We compare each of the gaits—Sim, Exp and Resim—to every other gait in the same class, and consider the performance ratios

δ​((i,j),X)≔[ΔCoM​(gaiti)ΔCoM​(gaitj)]i,j∈class​X\delta((i,j),X)\coloneqq\left[\tfrac{\Delta_{\text{CoM}}(\text{gait}_{i})}{\Delta_{\text{CoM}}(\text{gait}_{j})}\right]_{i,j\in\mathrm{class}\,X}

of their net displacements, which is a unitless metric correlating the gait performances within each class. For the experiments, each gait was tested through nine trials and the average CoM displacement is considered the experimental ground truth.

Refer to caption
Figure 8: These heat maps show the pairwise quotients of the gait performance ratios across the Experimental (Exp), Simulated (Sim), and Resimulated (Resim) classes. The statistical summary below the maps provides the mean value and standard deviation of cross-class comparisons, focusing on the consistency between experimental data and simulations.

Then, for corresponding gait pairings (i,j)(i,j) of two classes XX and YY, we consider the quotients

Ξ​((i,j),X,Y)≔δ​((i,j),X)δ​((i,j),Y)\Xi((i,j),X,Y)\coloneqq\tfrac{\delta((i,j),X)}{\delta((i,j),Y)}

of the performance ratios to determine how well differences in gaits were mapped from one class to another. This quotient naturally captures all disparities in the performance ratios caused by the class transitions. The data is displayed in a Heat Map in Fig. 8, with the mean and the STD being computed without the data of the major diagonal where gaits are compared against themselves.

The comparison between experimental, simulated, and resimulated gaits provides insight into the effectiveness of our simulation framework. The Exp/Sim comparison shows strong alignment (mean=1.0410\text{mean}=1.0410, STD=0.2974\text{STD}=0.2974), indicating that our initial simulations closely follow experimental data. However, the Exp/Resim comparison reveals slightly more variability (mean=1.0609\text{mean}=1.0609, STD=0.3718\text{STD}=0.3718), suggesting that resimulating based on tracked shapes introduces some divergence. Finally, the Sim/Resim comparison (mean=1.0932\text{mean}=1.0932, STD=0.4801\text{STD}=0.4801) suggests that while the two simulated datasets share general trends. The resimulation leads to a large variability, which is another indication that the robot does not achieve the desired shapes accurately.

Since the STD of all the between-class comparisons is significant, it is unlikely that individual trials will be accurately reproduced. But since the average of all trials was approximately 11, we conclude that Alg. 1 is effective at mapping displacement from simulation to reality in the aggregate, if the amount of trials is sufficiently large.

V-B Regularization

Another consideration is power efficiency. In robotic systems, battery power is at a premium and we tested the ability to improve gait efficiency.

Experimentally we know that energy is dissipated into the environment through the frictional resistance of the wheels. Therefore, we penalize the total energy dissipation (see Eq. (4)) in the simulation by introducing a penalty term to the objective of the optimization. That is, we minimize

−ΔCoM+c​ℰ,-\Delta_{\mathrm{CoM}}+c\,{\mathcal{E}},

where the Dissipation Coefficient c>0c>0 is a positive constant. This punishes energy loss while still attempting to maximize the CoM displacement. Additionally, we introduce the possibility for a constraint on the spatial frequency to generate gaits with a prescribed wavelength.

In a robotic system, a key metric for evaluating gait efficiency is the Cost of Transport

CoT≔P/m​g​v,\mathrm{CoT}\coloneqq\nicefrac{{P}}{{mgv}},

where PP denotes the power [J/s][\nicefrac{{\rm J}}{{\rm s}}], mm is the mass [kg][{\rm kg}], gg is gravitational acceleration [m/s2][\nicefrac{{\rm m}}{{{\rm s}^{2}}}] and vv is the velocity [m/s][\nicefrac{{\rm m}}{{\rm s}}]. To minimize the cost of transport, a gait needs to balance minimizing power consumption while achieving significant displacement. To measure this, we ran three trials with three gaits each and measured average power consumption. Simultaneously, the CoM displacement was recorded using the motion capture system. The data for each point displayed in Fig. 9 is a summary of information of 9 power and displacement trials.

Refer to caption
Figure 9: The effectiveness of penalizing dissipation in the experiment for three representative gaits (A, B, C). For each gait, the mean Cost of Transport is a function of the Dissipation Coefficient cc, and the standard deviation of the Cost of Transport is represented by error bars.

To determine if the optimization with an additional penalty term leads to efficiency improvements, we took a range of values for cc, from 0 to 2.52.5, for three representative gaits. For Gaits A and B, only the dissipation penalty is introduced to the optimization increase in weight. While for Gait C, a spatial frequency constraint was added. The experimental results indicate that the inclusion of penalty terms can improve the robot’s performance in terms of CoT.

When cc is low, the gaits obtained from optimization will barely change. On the other extreme, when the penalty on dissipation is increased beyond a threshold, the optimal actuation of the robot becomes very subtle, and thus the CoM displacement drops significantly, increasing the CoT. Additionally, as the motor movements become more subtle, static friction is encountered at a higher rate—increasing the required torque power consumption. However, when cc is between 1.0−2.01.0-2.0, a balance is struck between power consumption and displacement. The strongest evidence of this is in Gait A at c=1.7c=1.7 where we see a 23% improvement in the CoT as displayed (SI Movie 0:17-0:31). We see then, that penalizing dissipation is reflected in simulated results.

V-C Limitations and Challenges

There are two major limitations of the proposed approach in its current form: Slipping and reliance on motion capture data. Slipping typically occurs when the passive rubber wheels fail to maintain consistent traction with the surface, leading to deviations in the robot’s trajectory compared to the predicted one. This issue is primarily due to frictional inconsistencies on the experimental surface, as well as motor control limitations that prevent them from precisely reaching desired states. Slipping can introduce compounding errors in CoM displacement measurements, making it difficult to accurately resimulate the experimental trials. To mitigate this, future work could explore refining contact surface properties for consistent traction or upgrading the wheel design to something like a kirigami skin [36, 37]. Additionally, feedback control systems that dynamically adjust the robot’s movement in response to slipping as well as experimentation on various surfaces would provide valuable insights into gait efficiency and accuracy.

We note that the relatively high anisotropy ratio fitted to our experimental values betrays this slipping. That is, with a wheeled system that reinforces tangential motion, we expect a low (ϵ≈0\epsilon\approx 0) value. The lack of normal friction is then reflected in the simulation framework as higher tangential dissipation and can be effectively seen as a metric for the slipping frequency.

Moreover, our approach currently relies on external motion capture measurements for calibration. While our gait generation and improvement pipeline is powerful because it relies on only one physical parameter, ϵ\epsilon, this could be problematic in dynamic environments. We conducted experimental trials using a motion capture system to analyze the behavior of the current environment and evaluate its dissipative properties. Ideally, real-time environmental feedback would enable the system to compute the anisotropy ratio onboard and dynamically adjust its gait in response to changing conditions.

VI Conclusion and Outlook

We successfully demonstrated the implementation of a variational integrator based on geometric mechanics to model the undulatory motion of a snake-like robot. By characterizing the dissipative losses to the environment, the proposed method effectively characterizes the motion resulting from a sequence of robot configurations. It greatly simplifies the computational cost of modeling complex robot dynamics, particularly when compared to current physics engines [38]. We validated the accuracy of the simulation through the mean comparison of various gaits against experimental data and demonstrated strong correspondence of CoM displacement. Additionally, we demonstrated the ability to regularize gaits and improve efficiency as measured by the Cost of Transport by up to 23.2%23.2\%.

Although our results indicate an effective ability to qualitatively reproduce real-world results, there are many ways in which the gait generation pipeline could be improved. In its current state, Alg. 1 can be used to evaluate the effectiveness of a gait through a full forward simulation. This is only tractable due to the algorithm’s computational efficiency. For more reliable gait optimizations, more sophisticated optimization methods relying on, e.g., back-propagation should be explored. Nonetheless, we envision that the gaits we generate could be used as training data to convert a sequence of motor positions directly to a sequence of CoM displacement.

Neural network architectures such as Recurrent Neural Networks (RNNs) [39], Long Short-Term Memory networks (LSTMs)[40], and Transformers [41] have proven highly effective in sequence-to-sequence optimization tasks, making them well-suited for mapping motor position sequences to CoM displacement in robotic systems, enabling efficient and accurate gait prediction without the need for full forward simulations. In this way, we could infer the neural network to detect candidate gaits and evaluate performance before a robot is built. The model may not replicate every trial precisely; however, the experimental behavior can be accurately captured through a large sample size, as discussed in Sec. V-A.

This research opens new avenues for efficient and adaptable robotic locomotion, with potential applications in reinforcement learning, inverse design, and real-world robotic deployment across diverse environments.

Acknowledgment

We would like to thank the Naughton Undergraduate Fellowship program for supporting Cora Keogh. Additional support was provided by SideFX software.

References