Multi-legged locomotion with slipping

# Modeling multi-legged robot locomotion with slipping and its experimental
validation

Ziyou<EMAIL_ADDRESS>Robotics, University of Michigan  Dan Zhao
Mechanical Engineering, University of Michigan  Now at XPENG Robotics  Shai
Revzen  Robotics, University of Michigan  Electrical Engineering and Computer
Science, University of Michigan  Ecology and Evolutionary Biology, University
of Michigan  Applied Physics Program, University of Michigan

###### Abstract

Multi-legged robots with six or more legs are not in common use, despite
designs with superior stability, maneuverability, and a low numbers of
actuators being available for over 20 years. This may be in part due to the
difficulty in modeling multi-legged motion with slipping and producing
reliable predictions of body velocity. Here we present a detailed measurement
of the foot contact forces in a hexapedal robot with multiple sliding
contacts, and provide an algorithm for predicting these contact forces and the
body velocity. The algorithm relies on the recently published observation that
even while slipping, multi-legged robots are principally kinematic, and
employs a friction law ansatz that allows us to compute the shape-change to
body-velocity connection and the foot contact forces. This results in the
ability to simulate motion plans for a large number of potentially slipping
legs. In homogeneous environments this can run in (parallel) logarithmic time
of the planning horizon.

Key words: Legged locomotion, friction, modeling, ground contact force

###### Contents

  1. 1 Introduction
    1. 1.1 Multi-legged contact forces
  2. 2 Results
    1. 2.1 BigANT with force-torque sensors
      1. 2.1.1 BigANT: tripod gait
      2. 2.1.2 BigANT: metachronal gait
    2. 2.2 Demonstration with other legged systems
      1. 2.2.1 Multipods with 6, 8, 10, and 12 legs
      2. 2.2.2 Ghost Robotics Spirit
    3. 2.3 Analysis of running time
  3. 3 Discussion
  4. 4 Materials and methods : algorithm
    1. 4.1 Spring Support Model : finding the contacts
    2. 4.2 Local connection model : traction forces
      1. 4.2.1 Friction forces
    3. 4.3 Solving for planar body velocity
  5. 5 Materials and methods : robot experiments
    1. 5.1 Motion capture
      1. 5.1.1 Force/Torque measuring BigANT
      2. 5.1.2 Multipod robots
      3. 5.1.3 Ghost Robotics Spirit
    2. 5.2 Parameter fitting of $K$, $\mu$ and $\mathrm{H}$
  6. 6 Appendix
    1. 6.1 Contact detection
      1. 6.1.1 Normal force and planar moment balance solution
      2. 6.1.2 Special case: one leg in contact
      3. 6.1.3 Special case: two legs in contact
    2. 6.2 Planar motion

##  1 Introduction

Most recent research in legged robots has focused on bipedal or quadrupedal
robots, yet the vast majority of legged animal species use six or more legs,
are smaller, and therefore navigate a relatively much rougher terrain. Such
“multi-legged” systems – a term we refrain from using for quadrupeds and
bipeds – can exhibit complex tradeoffs between loads on the legs, and move
with substantial slipping at the feet. We are not aware of any multi-legged
robot which: (1) simultaneously measured all the ground contact forces
produced while moving and possibly slipping; (2) accurately modeled these
forces. Here we present such a robot, together with a fast algorithm for
modeling multi-legged system that predicts foot forces, and also provides (3)
the body velocity needed for producing motion plans. Recent work [1]
demonstrated that multi-legged systems, even while slipping, move in a
“geometric” or “principally kinematic” way, where the body velocity arises
through shape change rather than momentum buildup. That work showed that there
exists a “local connection”, in the sense of [2, 3], describing a shape
dependent linear relationship between the instantaneous time derivative of
shape and the instantaneous body velocity. While conventional Coulomb-friction
based modeling approaches cannot produce such a linear relationship, our
algorithm gives a local connection by construction. Furthermore, it scales
very favorably with the number of legs, and admits easy parallelization.

###  1.1 Multi-legged contact forces

Much of the work on multi-legged robots has focused on hexpedal robots,
starting with [4, 5]. Hexapods are appealing because they can have a tripod of
supporting legs while moving another tripod into place for support on the next
step. This static stability allows for the possibility of easier control, and
usable motion even when the robot is slipping. Significant prior work was done
on the hexpedal robots of the RHex family [6, 7], many variations of which
have been built over the past 20 years. Additional families of multi-legged
robots include the RoACH robots [8], the Sprawl robots [9], various multi-
legged robots that used “whegs” [10], and several studies of multi-legged
robots with large numbers of legs [11, 12, 13].

Unlike more simply supported robots, in a multi-legged system every closed
kinematic loop between body and ground can support an internal force that
produces no net body acceleration. This implies a sizable space of contact
forces which are unobservable from motion tracking, and cannot be measured
using the commonly used approach of placing a force plate for the robot to
move over [14]. This is because such a force plate only measures the total
wrench applied to it – not the interplay of all individual foot contact
forces, which may trade off in various ways step-to-step. The authors of [15]
measured the individual ground contact forces of RHex by installing each
individual leg with a 3d force sensor.

We set out to understand the results of [1] – why do multiple legs with
(presumably) Coulomb friction interactions with the substrate produce a local
connection between shape velocity and body velocity which is, algebraically at
least, incompatible with such a model. To address the question through
experiment we assembled a hexapedal robot with 6-Degree of Freedom(DoF) force-
torque sensors at the base of each leg, enabling the contact forces to be
measured directly [16]. To address the question through theory, we constructed
a highly simplified quasi-static model of body-ground interaction, and
replaced the Coulomb friction term which is linear in normal force but non-
smooth and non-linear in contact velocity with a friction ansatz that is
bilinear in normal force and contact velocity [17, Chapter 4].

Here we present a refined version of this model and show that it correctly
predicts the interplay of forces measured with the robot. We resolve the
seeming contradiction of having the “wrong” (ansatz) friction model produce
the “correct” forces, by showing that Coulomb friction and our ansatz produce
very similar motion predictions, for reasons we partially demonstrate in
section 3 and in [18]. From a computational perspective, we present a
numerical study demonstrating that our computation time is almost independent
of the number of legs, unlike the behavior of popular state-of-the-art robot
simulation tools.

##  2 Results

Our algorithm takes as inputs: (1) the positions $q_{j}$ and velocities
$\dot{q}_{j}$ of the robot’s feet; (2) the spring stiffness $k_{j}$ of each
leg; (3) the friction coefficients $\mu_{j}$ and friction anisotropy
${w}_{{xy,{j}}}$. As outputs it provides: (1) body height ${p}_{{z,{0}}}$,
pitch $\alpha_{x}$, and roll $\alpha_{y}$ slopes; (2) body velocities
${\dot{p}}_{{xy,{0}}}$ and $\dot{\theta}$; (3) 3D forces at the feet $F_{j}$.

###  2.1 BigANT with force-torque sensors

Figure 1: A tripod gait trajectory of BigANT. We plotted the trajectory
measured from motion capture (red), allowing it to be compared with the
integrated body velocity prediction from out simulation with classical Coulomb
friction (blue), and with our viscous-Coulomb friction law (orange). We
indicated the body location and orientation (rectangles) and indicated the
location of the robot body frame origin (crosses) at the beginning, half-way
point and the end. We plotted three Euler angles, three body frame velocities.
We calculated roll, pitch, supporting force and body velocity, and integrated
turning velocity to get yaw angle (viscous-Coulomb: orange, Coulomb friction:
blue). We used motion capture data as ground truth for kinematics (red). We
plotted the ground contact force in $x,y$-axis, and supporting force ratio in
$z$-axis of each leg (leg names are: HL hind left, ML mid left, FL front left,
HR hind right, MR mid right, FR front right). We used force torque sensor
measurements as ground truth for forces (red). We plotted the estimated
$z$-axis force ratio (purple), estimated $xy$-axis friction forces (viscous-
Coulomb: orange, Coulomb friction: blue).

To verify our algorithm’s accuracy, we built a version of the BigANT robot
with a 6-DoF force-torque sensor attached to each leg, and calibrated the
sensors to report the contact forces at the feet (see [16]). BigANT is a
hexpedal robot which has only one motor per leg. The leg drivetrain is a four-
bar linkage designed to provide advantageous foot clearance and gearing [17,
Chapter 2.2]. We recorded the robot motion using a marker-based motion
tracking system running at 100 frames per second(fps). In particular, we
measured the BigANT’s foot positions ($q_{j}$) in body frame, and we estimated
its body position and orientation from the markers attached to the robot
chassis. We obtained foot velocities $\dot{q}_{j}$ by differentiating $q_{j}$
using a Savitzky-Golay filter.

The remaining inputs to the algorithm were not so easy to determine. Because
our model is quasi-static, the mass plays no direct role, except for its
appearance in $Mg$ as the sum-total of normal forces at the feet. The force
and moment balance equations remain unchanged regardless of the units selected
for force, and these affect only $Mg$, the stiffnesses $K$, and the friction
coefficients $\mathrm{H}$. We therefore chose $Mg=1$. Using marker positions,
we estimated the robot body’s height, pitch and roll according to §4.1. We
estimated the spring constants $K$ and two anisotropic friction model
coefficients per leg ($\mu_{k}$ and ${w}_{{xy,{k}}}$ of eqn. 15) using least-
squares minimization of suitable penalty function. In total, we fitted 18
constant model parameters to predict a time series of six 3-dimensional leg
forces and 6-DoF body velocity measurements, thus there is little risk of
over-fitting.

####  2.1.1 BigANT: tripod gait

We first ran the robot through an alternating tripod gait driven with a
“Buehler clock” [6], with the steering strategy described in [19]. We
collected 21 tripod gait trials [20], with 4-5 cycles in each trial, and a
total of 102082 frames consisting of $84\pm 1$ cycles. The motions of the
shaft angles were scheduled to have a slow ground contact phase and a fast leg
reset phase. We show in figure 1 a comparison of forces and kinematics modeled
by our multi-contact algorithm with viscous-Coulomb friction, our algorithm
with classical Coulomb friction, and the experimental measurements. We
integrated body velocity and showed the robot trajectory in figure 1.

Because our physical modeling assumptions only define contact forces up to a
positive scale factor, we chose a single a positive scalar $\sigma(t)$ for
every time-step, such that

$\sigma:=\arg\min_{c}\sum_{k}(c|{\hat{F}}_{k}|-|F_{k}|)^{2}$

between the 12 dimensional prediction $\hat{F}$ and the measured horizontal
forces $F$ was minimized.

We reported prediction error statistics in figure 2. The mean and 1st and 3rd
quantile of the run time per frame for viscous-Coulomb friction is 0.19ms
(0.18, 0.24). When running a single approximation with the choice of
$\epsilon=10^{-5}$, the Coulomb friction solver takes 3.7ms (3.1, 3.9) per
frame. For a full set of iterations to convergence, Coulomb friction took
10.4ms (3.25,15.0) per frame.

Figure 2: Prediction error distributions for BigANT robot with tripod gait. We
plotted the distribution of measurement residuals from the mean (red) to
compare with residuals of the predictions of the spring supported model
(purple), or residuals of predictions from viscous-Coulomb (orange) and
Coulomb friction (blue).

####  2.1.2 BigANT: metachronal gait

Figure 3: We plotted: (A) Metachronal gait phase vs. motor shaft angle for all
six legs. (B) Magnitude of slipping velocity vs. magnitude of planar force
divided by normal force, overlaying points from all six feet. (C) Ground
contact forces. (D) Planar forces with mean subtracted (red) and model
prediction errors for planar forces (viscous-Coulomb: orange, Coulomb:blue).

We wanted to further study why the viscous-Coulomb friction model gave similar
body velocity and forces predictions as classic Coulomb friction model. Since
non-slip motions provide little insight into the question of which friction
force model to use, we developed a metachronal gait with exacerbated multi-
legged slipping events. Each foot contacted the ground for $2/3$ of a cycle
leading to four feet, two from each side, being in contact with the ground at
any time. To ensure that feet slip, we needed to ensure that the distances
between contacting feet change while in stance. We facilitated this by
ensuring that the contact feet have vastly incompatible velocities by making
the shaft angle a cubic function of time during stance. We collected 12
metachronal slipping gait trials, with the robot moving forward 4-6 cycles in
each. In total, the data consited of 43934 frames and $60\pm 1$ cycles. The
resulting gait produced much more slipping than the tripod gait, with slipping
velocities ranging in $(-51.8,111.8)$[mm/s] (5% quantile, 95% quantile). We
showed the gait in figure 3.

To determine whether viscous-Coulomb or classical Coulomb friction was
indicated in these data, we examined the force measurements from the slipping
gait. Plotting ${F}_{{x,{j}}}/{F}_{{z,{j}}}$ against ${v}_{{x,{j}}}$ (see
figure 3[B]) shows the expected structure for classical Coulomb friction,
namely a step function.

###  2.2 Demonstration with other legged systems

To test whether the proposed model generalizes to other legged systems, we
further examined using our model on Multipod robots with 6-12 leaf-spring legs
and an undulation gait, and on the commercially available quadruped Ghost
Robotics Spirit 40.

####  2.2.1 Multipods with 6, 8, 10, and 12 legs

We used the publicly available Multipod dataset [21, 22, 23] used in [1]. Each
contralateral pair of legs in a multipod has two DoF – yaw and roll – and the
roll (vertical) motion is further compounded with the spring flexing of the
leg itself (see figure 4). Here we used a slice of these data with the robot
running at frequency 0.3Hz and phase offset $1.35\pi$ to demonstrate our
algorithm. We assumed the mass of the robot is linear in the number of legs –
an explicit design feature of these robots – and set $mg=N$. We used $K=1$ as
the spring constant and isotropic friction model $\mu_{k}=1$ on all legs.

Figure 4: Modeling of Multipod with undulation gait at frequency 0.3Hz and
phase offset 1.35. We showed Multipod with 6-12 legs. We plotted the estimated
velocity and trajectory (blue), compared with motion tracking (red). Side
velocity plots have the same unit and scale as the heading velocity plots.

####  2.2.2 Ghost Robotics Spirit

Our physical modeling approach built upon the assumption that friction
dissipates the robot’s body momentum quickly in comparison to the time-scale
of gait motions. We intentionally selected a commercial quadruped, the Ghost
Robotics Spirit, where this assumption brakes down to test how well the
connection-based model could approximate the motion of such a quadruped. We
collected 921 frames comprising about 9 cycles of motion (see figure 5).
Because our model has no inertia, it tends to produce spurious high frequency
changes in its predictions. To obtain a more realistic time series, we added a
simple model of robot inertia in the form of a first-order IIR lowpass filter

$y_{n}=\gamma y_{n-1}+(1-\gamma)x_{n},$

where $x_{n}$ is our raw model prediction and $y_{n}$ is the filtered
prediction. We manually selected $\gamma=0.15$ to bring the power spectral
density (PSD, computed using scipy.signal.welch) of the estimated body
velocities close to that of the motion tracking derived velocities.

Figure 5: We showed the commercial quadrepedal robot, Spirit, used in the
experiment. We plotted the estimated velocity and trajectory (orange),
filtered estimation (blue) and motion tracking (red).

###  2.3 Analysis of running time

We compared the computation speed between our algorithm using the viscous-
Coulomb friction model and a widely-used physics simulation engine MuJoCo [24,
v2.2.1] using robot models with 3 to 50 legs. Since our focus is on multi-
legged contacts, our models consisted of a round disk with 3 to 50 legs
equally spaced on its circumference. We gave each leg two rotational DoF, a
vertical translation DoF, and limited leg motions so that their workspaces did
not overlap. We tested the execution time of both MuJoCo and our algorithm at
1000 randomly chosen poses and velocities for each number of legs, and re-
normalized the running time by dividing by the median execution time of the 3
legged case, to reveal how each simulation approach scaled with the number of
legs (see figure 6). While both algorithms reveal an increase in execution
times, our algorithm slows down by less than a factor of 3 with 50 legs,
compared with a factor of 13 for MuJoCo. This suggests that an optimized
implementation of our algorithm could be used for multi-legged motion planning
for any practical number of contacts.

Because we are using an inertia-free model of physics in the form of a local
connection, the body velocity at any instant is only a function of the shape
change and shape velocity at that instant. Hence, in a homogeneous environment
all time-steps of a motion plan can be computed in parallel. To demonstrate
the performance gains, we simulated 10,000 random poses and velocities of a
hexapod robot. We used $P=1,\cdots,4$ processors to compute the body velocity
matrices in parallel, then integrated them in a single linear process (note:
this over-estimates the parallelization overhead, since the product of $N$
matrices can be parallelized to take $\log_{2}N$ time, but was linear here).
In figure 6 we show that the algorithm parallelizes well, with the overhead at
four processors falling below 1.5, i.e. a net speedup of $4/1.5$.

Figure 6: (a) Plot of normalized run time of the multi-contact algorithm and
MuJoCo simulation versus number of legs. We plotted the distribution of time-
step computation times on $1000$ randomly initialized configurations for each
number of legs from $3$ to $50$. Plot indicates distribution percentiles $2.5$
to $97.5$ (lightly shaded); $25$ to $75$ (shaded); and median (dotted line).
The execution times are normalized relative to the median execution time of
each simulation on the $3$ leg case. The robot configurations consisted of a
disk with $N$ equally spaced legs on the rim as illustrated by examples with
$N=3$, $21$ and $42$. (b) Plot of parallelization overhead splitting the
algorithm over $M$ threads. The overhead is execution time times $M$ the
number of threads, in units of the median execution time on a single thread.
Perfectly parallelizable workloads give $1$ whereas unparalellizable workloads
give $M$. We plotted the workload distributions at $M=1\ldots 4$ for a
hexapod, running $100$ randomized trajectories each $10000$ time-steps long
(ribbon with same quantile as in (a)).

##  3 Discussion

Multi-legged robots (with six or more legs) are not widely studied in robotics
community. One reason might be that the complexity of modeling the multi-
contact ground interaction constrains both motion planning and simulations for
design. Motivated by a previous discovery [1] – that multi-legged robots move
as if they are governed by a local connection, i.e. quasi-statically – we
developed a simplified ground-interaction model and validated it
experimentally. Our algorithm consists of simulating a spring supported body
in a small-angle approximation for pitch and roll to obtain the vertical foot
loadings. We then introduced the viscous-Coulomb ansatz to replace classical
Coulomb friction in generating the horizontal forces to produce a linear set
of equations which can be solved to give rise to the local connection.

Our experimental verification demonstrated that while the actual contact
forces were, as expected, governed by classical Coulomb friction, our viscous-
Coulomb friction model gave equally good predictions of both contact forces
and body velocities, while computing 50 times faster for a hexapod. Our
algorithm scales to large numbers of contacts with virtually no change in
execution time, and parallelizes with very low overhead.

Figure 7: Contour of error between viscous-Coulomb approximation to Coulomb
friction around equilibrium velocity.

To understand how a system governed by sliding Coulomb friction can be modeled
by a viscous-Coulomb friction model, one may examine compare the relative
error of a viscous friction model to that of the “true” Coulomb friction.
Because both models are isotropic, we can assume without loss of generality
that the velocity is in the $x$ direction. Because both models are
homogeneous, we can assume without loss of generality that the speed is $1$.
What remains is to study the relative error of predicting the Coulomb friction
force for contact velocities close to $(1,0)$ and the prediction obtained by
using viscous drag model instead. In figure 7 we present the contours of
relative error when using a single viscous friction model instead of Coulomb
friction over the specified range of velocities. The plot demonstrates that
with $|\delta v|<0.2|v|$, the viscous-Coulomb force prediction for velocity
$v+\delta v$ will be within $2\%$ of the classical Coulomb friction force
prediction. The linearity between slipping velocity and friction forces was
also observed as an average relationship in numerical simulation [18] and
experiments [12].

We are thus left with the conclusion that a viscous-Coulomb ansatz model for
friction produces very similar predictions to those produced by the classical,
tribologically accurate Coulomb friction model. Comparing the motion
predictions obtained from both models, they are far more similar to each other
than either is to the measured motion, suggesting that the dominant error in
these models was not the use of an incorrect friction model. However, the
viscous-Coulomb model, in the context of our multi-contact algorithm provides
a significant performance boost. It is faster to compute; it scales better
with the number of contacts; and it is easier to parallelize.

From the perspective of physics, that our ansatz produces motion plans as
accurate as those produced by Coulomb friction, but also provably produces a
local connection and principally kinematic motion in the geometric mechanics
sense, provides further justification for the observation of [1] that local
connection models provide a framework that includes multi-legged locomotion.
While the local connection models of [1] were data-driven, here we have shown
that such models can be obtained using a principled modeling approach.

The algorithm we presented here provides merely a starting point – it is a
means for rapidly and accurately estimating multi-contact robot-environment
interactions. Such estimates are building blocks for motion planning, model
predictive control, design optimization, and many other potential
applications. The algorithm itself can be extended to include contacts with
non-flat ground, and the various quantities we estimated by fitting could be
converted to online, adaptively tuned quantities. We hope that such advances
will stimulate the adoption of multi-legged robots in field robotics, and
provide reliable and adaptable bio-inspired locomotion platforms.

### Acknowledgements

We would like to thank Andy Ruina for his key insight in explaining the
relationship between our viscous-Coulomb ansatz and Coulomb friction. We also
thank the many students who have worked in the BIRDS Lab at the University of
Michigan on collecting the large robot motion datasets used herein.

##  4 Materials and methods : algorithm

We propose an algorithm to estimate world frame body velocity from its body
shape and shape changing velocity at current time frame. The algorithm is
composed of two steps: (A) find which feet are in contact with the ground and
estimate their gravity induced loading using a spring support model; (B)
construct instantaneous local connection model to estimate the planar body
velocity.

The inputs to the spring support model (A) are: (1) the 3D positions of the
feet $q_{j}$ in the body frame; (2) the spring stiffness $k_{j}$ of each leg.
The outputs of the spring support model are: (1) body height ${p}_{{z,{0}}}$,
pitch $\alpha_{x}$, and roll $\alpha_{y}$ slopes; (2) gravity induced loading
on each foot ${F}_{{z,{j}}}$ and, implicit in that, which feet are in contact
with the ground.

Once the contacting feet are known, we solve for force and moment equilibrium
using a viscous-Coulomb friction ansatz which is bi-linear in ${F}_{{z,{j}}}$
and the foot sliding velocities in the world frame ${\dot{p}}_{{xy,{j}}}$,
providing an local connection model (B). The inputs to the connection model
are: (1) the 2D positions ${q}_{{xy,{j}}}$ and velocities
${\dot{q}}_{{xy,{j}}}$ of the feet in the body frame; (2) the friction
coefficients $\mu_{j}$ and friction anisotropy ${w}_{{xy,{j}}}$; (3) the
gravity induced loading ${F}_{{z,{j}}}$ computed in (A). The outputs of this
local connection model are body velocities: (1) body velocities
${p}_{{xy,{0}}}$ and $\dot{\theta}$; (2) 2D traction forces at the feet
${F}_{{xy,{j}}}$.

Suppose we are given a system with $N$ legs (or other contacts), indexed by
$j=1\ldots N$. The time varying foot positions in the body frame of reference
are given by $q_{j}\in\mathbb{R}^{3}$, $j=1\ldots N$. We assume the
transformation from body frame to world frame is given by a time varying rigid
body transformation $\Omega\in SE(3)$. The world frame foot positions $p_{j}$
and velocities $\dot{p}_{j}$ are

$\displaystyle p_{j}$ $\displaystyle:=\Omega q_{j}$ (1)
$\displaystyle\dot{p}_{j}$
$\displaystyle=\dot{\Omega}\,q_{j}+\Omega\dot{q}_{j}=\Omega\left[\Omega^{-1}\dot{\Omega}q_{j}+\dot{q}_{j}\right]$
(2)

Let $p_{0}$ represent the origin of the body frame. We assume a simplified
form for the rigid body transformation approximation $\Omega^{\prime}$, where
pitch $\alpha_{y}$ and roll $\alpha_{x}$ angles are small, so they can be
approximated by their first order Taylor approximation. We also assume the
rigid body motion is only time varying in the horizontal plane, i.e.
$\alpha_{x}$ $\alpha_{y}$ and ${p}_{{z,{0}}}$ vary so slowly that their
derivatives can be approximated by $0$. The detailed construction of
$\Omega^{\prime}$ and $\Omega^{\prime-1}\dot{\Omega}^{\prime}$ are in 4.

$\displaystyle\Omega^{\prime}:$
$\displaystyle=\begin{bmatrix}\cos\theta&-\sin\theta&\alpha_{x}&{p}_{{x,{0}}}\\\
\sin\theta&\cos\theta&-\alpha_{y}&{p}_{{y,{0}}}\\\
-\alpha_{x}&\alpha_{y}&1&{p}_{{z,{0}}}\\\ 0&0&0&1\end{bmatrix}$
$\displaystyle\Omega^{\prime-1}\dot{\Omega}^{\prime}:$
$\displaystyle=\begin{bmatrix}0&-\dot{\theta}&0&{\dot{p}}_{{x,{0}}}\\\
\dot{\theta}&0&0&{\dot{p}}_{{y,{0}}}\\\ 0&0&0&0\\\ 0&0&0&0\end{bmatrix}$ (3)

Because of these simplifying assumptions, we can decouple the movements in the
$xy$ plane, and the physical units of vertical and horizontal length are
decoupled. We use the planar rotation
${\mathrm{R}_{\theta}}:=\left[\begin{smallmatrix}\cos\theta&-\sin\theta\\\
\sin\theta&\cos\theta\end{smallmatrix}\right]$, and
$\mathrm{S}:=\left[\begin{smallmatrix}0&-1\\\ 1&0\end{smallmatrix}\right]$ to
represent foot position and velocity in world frame $xy$ plane with:

$\displaystyle{p}_{{xy,{j}}}$
$\displaystyle={\mathrm{R}_{\theta}}{q}_{{xy,{j}}}+{p}_{{xy,{0}}}$ (4)
$\displaystyle{\dot{p}}_{{xy,{j}}}$
$\displaystyle={\mathrm{R}_{\theta}}\left(\dot{\theta}{\mathrm{R}_{\theta}}^{-1}\mathrm{S}{\mathrm{R}_{\theta}}{q}_{{xy,{j}}}+{\dot{q}}_{{xy,{j}}}\right)$
$\displaystyle\qquad+{\dot{p}}_{{xy,{0}}}+[\alpha_{y},-\alpha_{x}]^{T}{q}_{{z,{j}}}$
(5)

###  4.1 Spring Support Model : finding the contacts

Figure 8: Visualization of search for contact state in a 2D “robot”. We
indicated the height and pitch (${p}_{{z,{0}}}$,$\alpha$) states searched
(labels 0-3 of (a)), and visualized the pose and contacts of the “robot” in
each of these states (corresponding labels of plots in (b)). Each “robot” leg
(zigzag lines in (b)) defines a corresponding codimension 1 plane (line here)
in ($z_{0}$,$\alpha$) at which it contacts the ground (colored lines in (a)
with color same as the leg in (b), (c)). At a ${p}_{{z,{0}}}$ above the plane,
the leg is in the air; below it the leg will be in contact and generate normal
forces. With each state being searched (number label in (a)), there is a
closed-form solution of the force equilibrium, which we connect to that state
with a line interval (black in (a)). If the equilibrium lies in the same
contact state the algorithm terminates (star; step 3). Otherwise the portion
of the line segment in another contact state is counter-factual (black dashed
in (a)). Instead, we switch to the new contact state and solve again. Each
such transition between contact states lies on a plane correspoding to the leg
that made contact (black dot in (a); circled leg in (b)).

In this section, we show how to decouple the roll, pitch and vertical (z-axis)
motion of the robot, and determine which legs are in contact and what
supporting force each leg generates. We model the robot as a “body plane”,
with each leg assumed to be a vertical spring attached to this plane. We
assume the system is at force and moment balance. A simplified version of this
model, without accounting for roll and pitch, can be found in [25, 17]. A
similar spring-leg model was used to study legged animals and robots [26, 13],
but they did not specify how to determine which legs are in contact.

Consider a pitch, roll, and height state
$\mathbf{s}=(\alpha_{x},\alpha_{y},{p}_{{z,{0}}})$. From (1), we have

${p}_{{z,{j}}}=-\alpha_{x}{q}_{{x,{j}}}+\alpha_{y}{q}_{{y,{j}}}+{q}_{{z,{j}}}+{p}_{{z,{0}}}={p}_{{z,{j}}}(\mathbf{s}).$
(6)

Taking $0$ to be the ground level, and up being the positive z-axis direction,
those ${p}_{{z,{j}}}<0$ are in contact with the ground. Assuming the normal
supporting force ${F}_{{z,{j}}}(\mathbf{s})$ is linearly dependent on
${p}_{{z,{j}}}$, we define the individual leg normal force and the resulting
planar moment function by,

$\displaystyle{F}_{{z,{j}}}(\mathbf{s})$
$\displaystyle:=\left\\{\begin{matrix}&-K_{j}\
{p}_{{z,{j}}}(\mathbf{s})\quad\text{if }{p}_{{z,{j}}}(\mathbf{s})<0\\\ \\\
&0\quad\text{otherwise}\end{matrix}\right.$
$\displaystyle{M}_{{x,{j}}}(\mathbf{s})$
$\displaystyle:=-{q}_{{y,{j}}}{F}_{{z,{j}}}(\mathbf{s})\qquad{M}_{{y,{j}}}(\mathbf{s}):={q}_{{x,{j}}}{F}_{{z,{j}}}(\mathbf{s}),$
(7)

and we denote the total force and moment by,

$\displaystyle F_{z}(\mathbf{s})=\sum_{i=1}^{N}{F}_{{z,{j}}}(\mathbf{s}),\quad
M_{x}(\mathbf{s})=\sum_{i=1}^{N}{M}_{{x,{j}}}(\mathbf{s}),$ $\displaystyle
M_{y}(\mathbf{s})=\sum_{i=1}^{N}{M}_{{y,{j}}}(\mathbf{s}).$ (8)

When $\alpha_{x}=\alpha_{y}=0$, the total normal force at height $z$ such that
${p}_{{z,{N_{k}}}}\leq-z<{p}_{{z,{N_{k}+1}}}$ is
$F_{z}([0,0,z])=\sum_{j=1}^{N_{k}}K_{j}(z+{p}_{{z,{j}}})$. WLOG, we can sort
${p}_{{z,{j}}}$ in non-increasing order, and we let $z=-{p}_{{z,{N_{k}}}}$,
starting with $N_{k}=1$, where only the lowest foot is in contact. We increase
$N_{k}$ until $F_{z}([0,0,-{p}_{{z,{N_{k}}}}])\leq
Mg<F_{z}(0,0,-{p}_{{z,{N_{k}+1}}})$, and then linearly interpolate with slope
$K_{N_{k}+1}$ to find $z^{*}$ such that force balance is achieved. For that
$z^{*}$, legs $k=1\ldots N_{k}$ are in contact with the ground. Throughout the
paper, we use the index $k$ to vary only over legs which are in contact with
the ground based on this criterion, and by ${F}_{{z,{k}}}$ the normal force of
those legs.

Next, we solve for the full state, $\mathbf{s}$, containing the small pitch
and roll angles, and the body height, maintaining vertical force balance, and
moment balance of the moments generated by the normal forces, i.e.
$F_{z}-Mg=M_{x}=M_{y}=0$. We start with an initial condition
$\mathbf{s}_{0}=(0,0,z^{*})$, with $F_{z}=Mg$. Taking
$\alpha_{x}$,$\alpha_{y}$ and $z^{*}$ as unknowns and holding the legs in
contact constant, these values are a solution $\mathbf{s}^{*}_{0}$ for a
3-dimensional linear system. We check whether the legs in contact at
$\mathbf{s}^{*}_{0}$ are the same as in $\mathbf{s}_{0}$; if so, then
$\mathbf{s}^{*}_{0}$ is the result from our model. If not, then we search
along the line segment starting at $\mathbf{s}_{0}$ and ending at
$\mathbf{s}^{*}_{0}$ for the first change in contacts, which must therefore
occur on a plane describing the contact condition for the first leg which
would change contact state going along this line segment. This transition
point is taken as $\mathbf{s}_{1}$, and the process repeats for the new legs
in contact. Because contact forces are zero on the corresponding contact
condition plane, $F_{z}$, $M_{x}$ and $M_{y}$ are continuous through the
change in contacting legs. The detailed expression of the equations is in
§6.1.1.

As the search iterates, we may encounter a state where only one or two legs
are in contact, and the linear force torque balance equation becomes under-
determined. To resolve these states we include an additional assumption – that
the origin of the body plane is the center of mass of the robot body. Under
this assumption, when there are fewer than three legs in contact, the COM
generates a moment around the contact point(s) and we tilt the body plane,
i.e. change $\alpha_{x}$ and $\alpha_{y}$, approximating the rotation this
moment would induce, until an additional leg contacts the ground.

We proceed to describe the tilting directions as if they were rotations with
an axis. However, the actual linear map they describe is a shearing operation
whose neutral plane intersects the $xy$ plane on a line containing the
rotation axis. When only one leg is in contact, the rotation is in the plane
containing the leg and the COM, around the contact point. When two legs are in
contact, the rotation is around the line connecting their contact points, and
in the direction of the moment the net $F_{z}$ generates around this line. The
detailed solution is in §6.1.2 and §6.1.3.

We used a 2D “robot” in $xz$-plane to visually illustrate our search algorithm
in figure 8. In this 2D case, the algorithm searches for robot height and
pitch using foot position in $x,z$ coordinates. The 3D model extends the
contact switching lines in 2D searching space to planes in 3D, and its
visualization can be found in figure 9.

###  4.2 Local connection model : traction forces

After knowing which legs are in contact and their gravity loading we solve for
the body planar velocity $({\dot{p}}_{{xy,{0}}},\dot{\theta})$, obtained by
imposing force and moment balance. While classical approaches suggest Coulomb
friction is the correct tribological model for sliding dry contacts, we show
that a viscous-Coulomb ansatz which is bilinear in both loading force and
sliding velocity makes for a linear system of equations that leads to a local
connection model.

####  4.2.1 Friction forces

The classical approaches to mechanics suggest that the contact between foot
and ground should be modeled by Coulomb friction (middle term below).

$\displaystyle{F}_{{xy,{k}}}$
$\displaystyle=-\frac{{\dot{p}}_{{xy,{k}}}}{\|{\dot{p}}_{{xy,{k}}}\|}\mu_{k}{F}_{{z,{k}}}=\mathrm{H}_{k}{\dot{p}}_{{xy,{k}}}$
(9)

The choice of
$\mathrm{H}_{k}=-\mu_{k}{F}_{{z,{k}}}/{\|{\dot{p}}_{{xy,{k}}}\|}$ would
provide equality, but this would produce the well-known problem of singularity
at ${\dot{p}}_{{xy,{k}}}=0$. Define $v_{k}:=\|{\dot{p}}_{{xy,{k}}}\|$. We
explore the tractability of alternative friction models using

$\displaystyle\mathrm{H}_{k}:=-\mu_{k}{F}_{{z,{k}}}\frac{\varepsilon+v_{k}}{\varepsilon+v_{k}^{2}}$
(10)

When $\varepsilon\to 0$, $\mathrm{H}_{k}\to$ Coulomb friction model; when
$\varepsilon\to\infty$, $\mathrm{H}_{k}\to-\mu_{k}{F}_{{z,{k}}}$, the friction
force becomes ${F}_{{xy,{k}}}=-\mu_{k}{F}_{{z,{k}}}{\dot{p}}_{{xy,{k}}}$, a
combination of viscous and Coulomb friction, depending on both slipping rate
and normal force.

We further deconstruct (9) in terms of $\dot{\theta}$, ${\dot{p}}_{{x,{0}}}$,
and ${\dot{p}}_{{y,{0}}}$:

$\displaystyle{F}_{{xy,{k}}}$
$\displaystyle=\mathrm{H}_{k}\left({\mathrm{R}_{\theta}}\left[\dot{\theta}{\mathrm{R}_{\theta}}^{-1}\mathrm{S}{\mathrm{R}_{\theta}}{q}_{{xy,{k}}}+{\dot{q}}_{{xy,{k}}}\right]+{\dot{p}}_{{xy,{0}}}\right)$
$\displaystyle=\left(\mathrm{H}_{k}\mathrm{S}{\mathrm{R}_{\theta}}{q}_{{xy,{k}}}\right)\dot{\theta}+\mathrm{H}_{k}\,{\dot{p}}_{{xy,{0}}}+\left(\mathrm{H}_{k}{\mathrm{R}_{\theta}}{\dot{q}}_{{xy,{k}}}\right)$
(11)

###  4.3 Solving for planar body velocity

From our quasi-static assumption, we have horizontal plane force and moment
balance, i.e. $\sum{F}_{{x,{k}}}=\sum{F}_{{y,{k}}}=\sum{M}_{{z,{k}}}=0$. From
horizontal force balance, using (11), we obtain two equations in
$\dot{\theta}$, ${\dot{p}}_{{x,{0}}}$, and ${\dot{p}}_{{y,{0}}}$

$\displaystyle 0=\sum_{k=1}^{N_{k}}{F}_{{xy,{k}}}$
$\displaystyle=\left(\sum_{k=1}^{N_{k}}\mathrm{H}_{k}\mathrm{S}{\mathrm{R}_{\theta}}{q}_{{xy,{k}}}\right)\dot{\theta}+\left(\sum_{k=1}^{N_{k}}\mathrm{H}_{k}\right){\dot{p}}_{{xy,{0}}}$
$\displaystyle\qquad+\left(\sum_{k=1}^{N_{k}}\mathrm{H}_{k}{\mathrm{R}_{\theta}}{\dot{q}}_{{xy,{k}}}\right)$
(12)

The moment exerted by a leg is given by:

$\displaystyle{M}_{{z,{k}}}$
$\displaystyle={p}_{{xy,{k}}}^{\mathsf{T}}\mathrm{S}{F}_{{xy,{k}}}=\left({p}_{{xy,{k}}}^{\mathsf{T}}\mathrm{S}\mathrm{H}_{k}\mathrm{S}{\mathrm{R}_{\theta}}{q}_{{xy,{k}}}\right)\dot{\theta}+$
$\displaystyle\left({p}_{{xy,{k}}}^{\mathsf{T}}\mathrm{S}\mathrm{H}_{k}\right){\dot{p}}_{{xy,{0}}}+\left({p}_{{xy,{k}}}^{\mathsf{T}}\mathrm{S}\mathrm{H}_{k}{\mathrm{R}_{\theta}}{\dot{q}}_{{xy,{k}}}\right)$
(13)

Giving the obvious third and final equation:

$\displaystyle 0$
$\displaystyle=\sum_{k=1}^{N_{k}}{M}_{{z,{k}}}=\left(\sum_{k=1}^{N_{k}}{p}_{{xy,{k}}}^{\mathsf{T}}\mathrm{S}\mathrm{H}_{k}\mathrm{S}{\mathrm{R}_{\theta}}{q}_{{xy,{k}}}\right)\dot{\theta}$
$\displaystyle+\left(\sum_{k=1}^{N_{k}}{p}_{{xy,{k}}}^{\mathsf{T}}\mathrm{S}\mathrm{H}_{k}\right){\dot{p}}_{{xy,{0}}}+\left(\sum_{k=1}^{N_{k}}{p}_{{xy,{k}}}^{\mathsf{T}}\mathrm{S}\mathrm{H}_{k}{\mathrm{R}_{\theta}}{\dot{q}}_{{xy,{k}}}\right)$
(14)

In the case when $\varepsilon\to\infty$, $\mathrm{H}_{k}$ being rate
${\dot{p}}_{{xy,{k}}}$ independent, the three force and moment balance
equations are linear in the body velocity ${\dot{p}}_{{xy,{0}}}$,
$\dot{\theta}$ and foot velocity in body frame ${\dot{q}}_{{xy,{k}}}$. One
could solve the system by 3d matrix inversion. The detailed expression of the
solution is derived in §6.2.

In addition to classic Coulomb friction and viscous friction, we consider the
possibility that $\mathrm{H}_{k}$ can be dependent on slipping direction,
modeling forces generated by a wheel, skate, claw, or otherwise non-isotropic
frictional contact. We consider an anisotropic viscous friction model, where
$\mathrm{H}_{k}$ is a symmetric positive semidefinite matrix,
$\mathrm{H}_{k}(q):={\mathrm{R}_{\theta}}\mathrm{H}_{q,k}(q){\mathrm{R}_{\theta}}^{-1}$
taken to be independent of ${\dot{p}}_{{xy,{k}}}$, but (possibly non-linearly)
dependent on all elements of $q$. We assume that each contact is associated
with an enhanced traction direction and associated magnitude, expressed in
body coordinates as a vector ${w}_{{xy,{k}}}$, defined as:

$\mathrm{H}_{k}:=-\mu_{k}{F}_{{z,{k}}}{\mathrm{R}_{\theta}}(\mathrm{I}_{2}+{w}_{{xy,{k}}}{w}_{{xy,{k}}}^{\mathsf{T}}){\mathrm{R}_{\theta}}^{-1}$
(15)

This changes the circular cross-section of the friction cone into an
ellipsoidal one. Even with this dependence, the equations (14) and (12) are
still linear in the velocities ${\dot{p}}_{{xy,{0}}}$, $\dot{\theta}$ and
${\dot{q}}_{{xy,{k}}}$. Similar to §4.3, body velocity ${\dot{p}}_{{xy,{0}}}$,
$\dot{\theta}$ can still be solved linearly with respect to shape changing
velocity ${\dot{q}}_{{xy,{k}}}$, giving a general form:

$\displaystyle{\mathrm{R}_{\theta}}^{-1}{\dot{p}}_{{xy,{0}}}$
$\displaystyle=:\sum_{k}{A}_{{xy,{k}}}(q){\dot{q}}_{{xy,{k}}},$
$\displaystyle\dot{\theta}$
$\displaystyle=:-\sum_{k}A_{\theta,k}(q){\dot{q}}_{{xy,{k}}},$

Where the $A_{\cdot,\cdot}(q)$ matrices is kinematic term in the
reconstruction equation of geometric mechanics.

##  5 Materials and methods : robot experiments

###  5.1 Motion capture

We used 10 Qualisys Oqus-310+ cameras as marker-based tracking system, running
at 100 fps with software QTM 2.17 build 4000.

####  5.1.1 Force/Torque measuring BigANT

The BigANT robot has six 1 DoF legs, each actuated by a servo motor (Robotis
Dynamixel MX106), via a 4-bar mechanism we chose for both foot clearance and
advantageous instantaneous gearing ratios. The detailed design and
characteristics can be found in [17, Chapter 2.2]. We manufactured the legs
from foamcore (Elmer’s Products Inc. 3/8” foam board) and fiber reinforced
tape (3M Scotch #8959) using the “plates and reinforced flexures” (PARF)
technique of [27]. We laser-cut the base plate for the BigANT chassis from a
1/4” ABS plastic plate. We then installed an 6 DoF Force/Torque sensor (ATI
Gamma F/T sensor) underneath each leg, and used the vendor’s wireless F/T
sensor system to communicate the measurements to the controlling host computer
at 100 Hz. We calibrated the F/T sensors according to [16].

####  5.1.2 Multipod robots

We constructed the Multipod robots 3 to 6 pairs of 2-DoF segments with spring
steel cantilevers as legs. The mechanical design and characteristics can be
found in [17, Chapter 2.3]. We provide detailed documentation of the motion
capture data in [1], and the dataset itself can be found in [21, 22, 23]. To
model the body motion of unloaded spring legs, we computed the location of
each foot relative to the rigid segment to which its spring was connected
using motion tracking frames in which the leg was unloaded. We then
extrapolated this unloaded position to the frames where the leg was loaded.

####  5.1.3 Ghost Robotics Spirit

We used a Spirit 40 v1.0 Proto Q-UGV Robot from Ghost robotics, operated
through a Samsung Galaxy A50 with on board firmware version 0.12.6.

###  5.2 Parameter fitting of $K$, $\mu$ and $\mathrm{H}$

We estimated the spring constants $K$ by minimizing the L2-norm difference
between modeled and measured ${F}_{{z,{}}}$ distribution among contacting
legs, while adding the coefficient of variation of $K$ between legs as a
regularization penalty. We assumed anisotropic friction coefficients
$\mathrm{H}_{k}$ (see eqn. 15), and inferred the parameters $\mu_{k}$ and
${w}_{{xy,{k}}}$ for each leg by minimizing L2-norm error between measured
forces and forces calculated from slipping velocity measured by motion
capture. We used scipy.optimize.least_squares for both of these parameter
estimation minimizations.

We solved for classical Coulomb friction with scipy.optimize.root using the LM
algorithm, with the solution of the previous time-step as an initial guess.
Because Coulomb friction is non-smooth, we approximated the solution starting
from $\varepsilon=10^{-5}$ (from eqn. 10), using each solution as an initial
condition for another solution with smaller $\varepsilon$ until the relative
change in L2-norm of two consecutive solutions differed less than $10^{-3}$ –
a threshold less than $5\%$ of the median ground speed measured. In the very
rare cases (0.12% of the BigANT tripod gait dataset) where the Coulomb
friction solver failed to converge, we used the ground truth velocity as the
initial condition to obtain the Coulomb friction solution, and these
converged.

## References

## References

  * [1] D. Zhao, B. Bittner, G. Clifton, N. Gravish, and S. Revzen, “Walking is like slithering: A unifying, data-driven view of locomotion,” _Proceedings of the National Academy of Sciences_ , vol. 119, no. 37, p. e2113222119, 2022. [Online]. Available: https://www.pnas.org/doi/abs/10.1073/pnas.2113222119
  * [2] R. L. Hatton and H. Choset, “Geometric motion planning: The local connection, stokes’ theorem, and the importance of coordinate choice,” _The International Journal of Robotics Research_ , vol. 30, no. 8, pp. 988–1014, 2011\. [Online]. Available: https://doi.org/10.1177/0278364910394392
  * [3] J. Ostrowski and J. Burdick, “The geometric mechanics of undulatory robotic locomotion,” _The International Journal of Robotics Research_ , vol. 17, no. 7, pp. 683–701, 1998. [Online]. Available: https://doi.org/10.1177/027836499801700701
  * [4] D. Gorinevsky and A. Y. Shneider, “Force control in locomotion of legged vehicles over rigid and soft surfaces,” _The International Journal of Robotics Research_ , vol. 9, no. 2, pp. 4–23, 1990.
  * [5] R. Carlton and S. Bartholet, “The evolution of the application of mobile robotics to nuclear facility operations and maintenance,” vol. 4, pp. 720–726, 1987.
  * [6] U. Saranli, M. Buehler, and D. E. Koditschek, “Rhex: A simple and highly mobile hexapod robot,” _The International Journal of Robotics Research_ , vol. 20, no. 7, pp. 616–631, 2001. [Online]. Available: https://doi.org/10.1177/02783640122067570
  * [7] I.-C. Chang, C.-H. Hsu, H.-S. Wu, and P.-C. Lin, “An analysis of the rolling dynamics of a hexapod robot using a three-dimensional rolling template,” _NONLINEAR DYNAMICS_ , vol. 109, no. 2, pp. 631–655, JUL 2022.
  * [8] D. Zarrouk, D. W. Haldane, and R. S. Fearing, “Dynamic legged locomotion for palm-size robots,” in _Micro- and Nanotechnology Sensors, Systems, and Applications VII_ , T. George, A. K. Dutta, and M. S. Islam, Eds., vol. 9467, International Society for Optics and Photonics. SPIE, 2015, p. 94671S. [Online]. Available: https://doi.org/10.1117/12.2176780
  * [9] S. Kim, J. E. Clark, and M. R. Cutkosky, “isprawl: Design and tuning for high-speed autonomous open-loop running,” _The International Journal of Robotics Research_ , vol. 25, no. 9, pp. 903–912, 2006. [Online]. Available: https://doi.org/10.1177/0278364906069150
  * [10] D. Zarrouk and R. S. Fearing, “Controlled in-plane locomotion of a hexapod using a single actuator,” _IEEE Transactions on Robotics_ , vol. 31, no. 1, pp. 157–167, 2015.
  * [11] S. Aoi, Y. Egi, and K. Tsuchiya, “Instability-based mechanism for body undulations in centipede locomotion,” _Phys. Rev. E_ , vol. 87, p. 012717, Jan 2013. [Online]. Available: https://link.aps.org/doi/10.1103/PhysRevE.87.012717
  * [12] B. Chong, J. He, S. Li, E. Erickson, K. Diaz, T. Wang, D. Soto, and D. Goldman, “Self propulsion via slipping: frictional resistive force theory for multi-legged locomotors,” 07 2022.
  * [13] B. Chong, Y. O. Aydin, J. M. Rieser, G. Sartoretti, T. Wang, J. Whitman, A. Kaba, E. Aydin, C. McFarland, H. Choset, and D. I. Goldman, “A general locomotion control framework for serially connected multi-legged robots,” _CoRR_ , vol. abs/2112.00662, 2021. [Online]. Available: https://arxiv.org/abs/2112.00662
  * [14] H. Komsuoḡlu, K. Sohn, R. J. Full, and D. E. Koditschek, “A physical model for dynamical arthropod running on level ground,” in _Experimental Robotics_ , O. Khatib, V. Kumar, and G. J. Pappas, Eds. Berlin, Heidelberg: Springer Berlin Heidelberg, 2009, pp. 303–317.
  * [15] C.-J. Kao, C.-S. Chen, and P.-C. Lin, “Reactive force analysis and modulation of the individual legs in a running hexapod robot,” in _2019 IEEE/ASME International Conference on Advanced Intelligent Mechatronics (AIM)_ , 2019, pp. 370–375.
  * [16] Z. Wu and S. Revzen, “In-situ calibration of six-axis force/torque transducers on legged robot,” University of Michigan, Tech. Rep., May 2023.
  * [17] D. Zhao, “Locomotion of low-dof multi-legged robots,” Ph.D. dissertation, University of Michigan, 2021.
  * [18] Z. Wu, D. Zhao, and S. Revzen, “Coulomb Friction Crawling Model Yields Linear Force–Velocity Profile,” _Journal of Applied Mechanics_ , vol. 86, no. 5, 03 2019, 054501. [Online]. Available: https://doi.org/10.1115/1.4042696
  * [19] D. Zhao and S. Revzen, “Multi-legged steering and slipping with low dof hexapod robots,” _BIOINSPIRATION & BIOMIMETICS_, vol. 15, no. 4, JUL 2020\.
  * [20] BIRDS-Lab, “Force/torque bigant processed data,” 2023.
  * [21] ——, “Birds lab multipod robot motion tracking data - processed data and code,” 2021.
  * [22] ——, “Birds lab multipod robot motion tracking data - raw dataset,” 2021.
  * [23] ——, “Birds lab multipod robot motion tracking data - videos,” 2021.
  * [24] E. Todorov, T. Erez, and Y. Tassa, “Mujoco: A physics engine for model-based control,” in _2012 IEEE/RSJ International Conference on Intelligent Robots and Systems_ , 2012, pp. 5026–5033.
  * [25] D. Zhao, V. Sachdva, and S. Revzen, “Modeling multilegged locomotion: the friction dominated limit,” in _Integrative and Comparative Biology_ , 2018\.
  * [26] J. R. Usherwood and B. J. H. Smith, “The grazing gait, and implications of toppling table geometry for primate footfall sequences,” _Biology Letters_ , vol. 14, no. 5, p. 20180137, 2018. [Online]. Available: https://royalsocietypublishing.org/doi/abs/10.1098/rsbl.2018.0137
  * [27] I. Fitzner, Y. Sun, V. Sachdeva, and S. Revzen, “Rapidly prototyping robots: Using plates and reinforced flexures,” _IEEE Robotics & Automation Magazine_, vol. 24, no. 1, pp. 41–47, 2017.

##  6 Appendix

###  6.1 Contact detection

The normal force exerted by a contacting leg:

$\displaystyle{F}_{{z,{k}}}=K_{k}(-\alpha_{x}{q}_{{x,{k}}}+\alpha_{y}{q}_{{y,{k}}}+{q}_{{z,{k}}}+{p}_{{z,{0}}})$

The $x,y$ moment resulted from the normal force:

${M}_{{x,{k}}}=-{q}_{{y,{k}}}K_{k}(-\alpha_{x}{q}_{{x,{k}}}+\alpha_{y}{q}_{{y,{k}}}+{q}_{{z,{k}}}+{p}_{{z,{0}}})\qquad{M}_{{y,{k}}}={q}_{{x,{k}}}K_{k}(-\alpha_{x}{q}_{{x,{k}}}+\alpha_{y}{q}_{{y,{k}}}+{q}_{{z,{k}}}+{p}_{{z,{0}}})$

####  6.1.1 Normal force and planar moment balance solution

When there are three or more legs in contact, the linear system to solve for
normal force and planar moment balance is:

$\displaystyle\left[\begin{array}[]{c}\sum_{k=1}^{N_{k}}{F}_{{z,{k}}}\\\
\sum_{k=1}^{N_{k}}{M}_{{xy,{k}}}\end{array}\right]$
$\displaystyle=\left[\begin{array}[]{c}-Mg\\\ \mathbf{0}\end{array}\right]$
$\displaystyle\sum_{k=1}^{N_{k}}\left[\begin{array}[]{c}-K_{k}{q}_{{x,{k}}}\alpha_{x}+K_{k}{q}_{{y,{k}}}\alpha_{y}+K_{k}{p}_{{z,{0}}}+K_{k}{q}_{{z,{k}}}\\\
K_{k}{q}_{{x,{k}}}{q}_{{y,{k}}}\alpha_{x}-K_{k}{q}_{{y,{k}}}^{2}\alpha_{y}-K_{k}{q}_{{y,{k}}}{p}_{{z,{0}}}-K_{k}{q}_{{y,{k}}}{q}_{{z,{k}}}\\\
-K_{k}{q}_{{x,{k}}}^{2}\alpha_{x}+K_{k}{q}_{{x,{k}}}{q}_{{y,{k}}}\alpha_{y}+K_{k}{q}_{{x,{k}}}{p}_{{z,{0}}}+K_{k}{q}_{{x,{k}}}{q}_{{z,{k}}}\end{array}\right]$
$\displaystyle=\left[\begin{array}[]{c}-Mg\\\ 0\\\ 0\end{array}\right]$
$\displaystyle\sum_{k=1}^{N_{k}}\left[\begin{array}[]{ccc}-K_{k}{q}_{{x,{k}}}&K_{k}{q}_{{y,{k}}}&K_{k}\\\
-K_{k}{q}_{{x,{k}}}{q}_{{y,{k}}}&K_{k}{q}_{{y,{k}}}^{2}&-K_{k}{q}_{{y,{k}}}\\\
-K_{k}{q}_{{x,{k}}}^{2}&K_{k}{q}_{{x,{k}}}{q}_{{y,{k}}}&K_{k}{q}_{{x,{k}}}\end{array}\right]\left[\begin{array}[]{c}\alpha_{x}\\\
\alpha_{y}\\\ {p}_{{z,{0}}}\end{array}\right]$
$\displaystyle=\left[\begin{array}[]{c}-Mg-K_{k}{q}_{{z,{k}}}\\\
K_{k}{q}_{{y,{k}}}{q}_{{z,{k}}}\\\
-K_{k}{q}_{{x,{k}}}{q}_{{z,{k}}}\end{array}\right]$

Let the previous state to be
$\mathbf{s}_{0}=[\alpha_{x}^{0},\alpha_{y}^{0},{p}_{{z,{0}}}^{0}]$ and current
solution $\mathbf{s}_{1}=[\alpha_{x}^{1},\alpha_{y}^{1},{p}_{{z,{0}}}^{1}]$,
the next state is obtained by $\mathbf{s}_{2}=\mathbf{s}_{0}+t(s_{1}-s_{0})$,
where $t$ is the minimum solution among all $t_{j}$ for which $t_{j}\in[0,1]$.
Each $t_{j}$ is calculated from:

${F}_{{z,{j}}}(\mathbf{s}_{2})=K_{j}(-((1-t_{j})\alpha_{x}^{0}+t\alpha_{x}^{1}){q}_{{x,{j}}}+((1-t_{j})\alpha_{y}^{0}+t_{j}\alpha_{y}^{1}){q}_{{y,{j}}}+{q}_{{z,{j}}}+(1-t_{j}){p}_{{z,{0}}}^{0}+t_{j}{p}_{{z,{0}}}^{1}))=0$

Therefore,

$t_{j}=\frac{\alpha_{x}^{0}{q}_{{x,{j}}}-\alpha_{y}^{0}{q}_{{y,{j}}}-{q}_{{z,{j}}}-{p}_{{z,{0}}}^{0}}{-(\alpha_{x}^{1}-\alpha_{x}^{0}){q}_{{x,{j}}}+(\alpha_{y}^{1}-\alpha_{y}^{0}){q}_{{y,{j}}}+({p}_{{z,{0}}}^{1}-{p}_{{z,{0}}}^{0})},\qquad
t=\min\\{t_{1},\cdots t_{N}\\}$

Figure 9: Searching iterations for a 3D robot. The iterations labeled 0-3 in
the search space (left) correspond to the labels below the robot (right). The
colored surfaces are the contact switching surfaces for legs that has active
contact switches in this iteration (the color of a surface corresponds to the
color of a leg on the robot). It starts with the initial condition, where
z-axis force balance is satisfied. Then it solves force and moment balance for
pitch, roll and height, where the grey dashed line is the direction of the
solution given current contact indices. It goes along that direction, and
stops at the first intersection with a color line, where a leg switches its
contact state. The iteration continues, until a solution is found within its
current polygon, i.e. contact indices do not change. The green circle
highlights the legs switching contact at that iteration.

####  6.1.2 Special case: one leg in contact

When there is only one leg in contact with the ground, denoting it with index
1, we rotate the body frame about the contacting point along the normal
direction connecting the contacting leg and origin. The direction of rotation
is the direction of the torque resulted by robot gravity at the contacting
point. The direction towards next state is
$s_{1}-s_{0}=[-{p}_{{y,{1}}}^{2}-{p}_{{x,{1}}}^{2},{p}_{{y,{1}}},-{p}_{{x,{1}}}]$.
The magnitude of state change $t$ is the minimum amount that results in
another contact leg, which can be solved the same as previous section.

####  6.1.3 Special case: two legs in contact

When there are only two legs in contact with the ground, denoting it with
index 1 and 2, we rotate the body frame about the line connecting these two
contact points. The direction of rotation is the direction of the torque
resulted by robot gravity about the line. Therefore, let
$a={p}_{{y,{1}}}-{p}_{{y,{2}}},b={p}_{{x,{2}}}-{p}_{{x,{1}}},c={p}_{{y,{1}}}{p}_{{x,{2}}}-{p}_{{y,{2}}}{p}_{{x,{1}}}$,
the direction of the change in pitch and roll becomes
$n=[-ca/\sqrt{c^{2}a^{2}+c^{2}b^{2}},cb/\sqrt{c^{2}a^{2}+c^{2}b^{2}}]$.
Together with the change in height, the direction towards next state is:

$s_{1}-s_{0}=[(-n_{0}^{2}-n_{1}^{2})|c|/\sqrt{(a^{2}+b^{2})},n_{0},n_{1}]$

The magnitude of state change $t$ is again the minimum amount that results in
another contact leg, which can be solved the same as previous section.

###  6.2 Planar motion

The planar force exerted by a leg:

$\displaystyle{F}_{{xy,{k}}}$
$\displaystyle=\mathrm{H}_{k}\left({\mathrm{R}_{\theta}}\left[\dot{\theta}{\mathrm{R}_{\theta}}^{-1}\mathrm{S}{\mathrm{R}_{\theta}}{q}_{{xy,{k}}}+{\dot{q}}_{{xy,{k}}}\right]+{\dot{p}}_{{xy,{0}}}\right)$
$\displaystyle=\left(\mathrm{H}_{k}\mathrm{S}{\mathrm{R}_{\theta}}{q}_{{xy,{k}}}\right)\dot{\theta}+\mathrm{H}_{k}\,{\dot{p}}_{{xy,{0}}}+\left(\mathrm{H}_{k}{\mathrm{R}_{\theta}}{\dot{q}}_{{xy,{k}}}\right)$
$\displaystyle=\left({\mathrm{R}_{\theta}}\mathrm{H}_{q,k}\mathrm{S}{q}_{{xy,{k}}}\right)\dot{\theta}+{\mathrm{R}_{\theta}}\mathrm{H}_{q,k}{\mathrm{R}_{\theta}}^{-1}\,{\dot{p}}_{{xy,{0}}}+\left({\mathrm{R}_{\theta}}\mathrm{H}_{q,k}{\dot{q}}_{{xy,{k}}}\right)$

The $z-$moment exerted by a leg:

$\displaystyle{M}_{{z,{k}}}$
$\displaystyle:=({p}_{{xy,{k}}}-{p}_{{xy,{0}}})^{\mathsf{T}}\mathrm{S}{F}_{{xy,{k}}}$
$\displaystyle=\left(({p}_{{xy,{k}}}-{p}_{{xy,{0}}})^{\mathsf{T}}\mathrm{S}\mathrm{H}_{k}\mathrm{S}{\mathrm{R}_{\theta}}{q}_{{xy,{k}}}\right)\dot{\theta}+\left(({p}_{{xy,{k}}}-{p}_{{xy,{0}}})^{\mathsf{T}}\mathrm{S}\mathrm{H}_{k}\right){\dot{p}}_{{xy,{0}}}+\left(({p}_{{xy,{k}}}-{p}_{{xy,{0}}})^{\mathsf{T}}\mathrm{S}\mathrm{H}_{k}{\mathrm{R}_{\theta}}{\dot{q}}_{{xy,{k}}}\right)$
$\displaystyle=\left(({p}_{{xy,{k}}}-{p}_{{xy,{0}}})^{\mathsf{T}}\mathrm{S}{\mathrm{R}_{\theta}}\mathrm{H}_{q,k}\mathrm{S}{q}_{{xy,{k}}}\right)\dot{\theta}+\left(({p}_{{xy,{k}}}-{p}_{{xy,{0}}})^{\mathsf{T}}\mathrm{S}{\mathrm{R}_{\theta}}\mathrm{H}_{q,k}{\mathrm{R}_{\theta}}^{-1}\right){\dot{p}}_{{xy,{0}}}$
$\displaystyle\quad+\left(({p}_{{xy,{k}}}-{p}_{{xy,{0}}})^{\mathsf{T}}\mathrm{S}{\mathrm{R}_{\theta}}\mathrm{H}_{q,k}{\dot{q}}_{{xy,{k}}}\right)$
$\displaystyle=\left({q}_{{xy,{j}}}^{\mathsf{T}}\mathrm{S}\mathrm{H}_{q,k}\mathrm{S}{q}_{{xy,{k}}}\right)\dot{\theta}+\left({q}_{{xy,{j}}}^{\mathsf{T}}\mathrm{S}\mathrm{H}_{q,k}{\mathrm{R}_{\theta}}^{-1}\right){\dot{p}}_{{xy,{0}}}+\left({q}_{{xy,{j}}}^{\mathsf{T}}\mathrm{S}\mathrm{H}_{q,k}{\dot{q}}_{{xy,{k}}}\right)$

The linear system to solve for planar force and moment balance is:

$\displaystyle\left[\begin{array}[]{ccc}\sum_{k=1}^{N_{k}}{F}_{{xy,{k}}}\\\
\sum_{k=1}^{N_{k}}{M}_{{z,{k}}}\end{array}\right]$ $\displaystyle=0$
$\displaystyle\sum_{k=1}^{N_{k}}\left[\begin{array}[]{cc}{\mathrm{R}_{\theta}}\mathrm{H}_{q,k}{\mathrm{R}_{\theta}}^{-1}&{\mathrm{R}_{\theta}}\mathrm{H}_{q,k}\mathrm{S}{q}_{{xy,{k}}}\\\
{q}_{{xy,{j}}}^{\mathsf{T}}\mathrm{S}\mathrm{H}_{q,k}{\mathrm{R}_{\theta}}^{-1}&{q}_{{xy,{j}}}^{\mathsf{T}}\mathrm{S}\mathrm{H}_{q,k}\mathrm{S}{q}_{{xy,{k}}}\end{array}\right]\left[\begin{array}[]{c}{\dot{p}}_{{xy,{0}}}\\\
\dot{\theta}\end{array}\right]$
$\displaystyle=\sum_{k=1}^{N_{k}}\left[\begin{array}[]{c}{\mathrm{R}_{\theta}}\mathrm{H}_{q,k}{\dot{q}}_{{xy,{k}}}\\\
{q}_{{xy,{j}}}^{\mathsf{T}}\mathrm{S}\mathrm{H}_{q,k}{\dot{q}}_{{xy,{k}}}\end{array}\right]$
$\displaystyle\left[\begin{array}[]{cc}{\mathrm{R}_{\theta}}&0\\\
0&1\end{array}\right]\sum_{k=1}^{N_{k}}\left[\begin{array}[]{cc}\mathrm{H}_{q,k}&\mathrm{H}_{q,k}\mathrm{S}{q}_{{xy,{k}}}\\\
{q}_{{xy,{j}}}^{\mathsf{T}}\mathrm{S}\mathrm{H}_{q,k}&{q}_{{xy,{j}}}^{\mathsf{T}}\mathrm{S}\mathrm{H}_{q,k}\mathrm{S}{q}_{{xy,{k}}}\end{array}\right]\left[\begin{array}[]{c}{\mathrm{R}_{\theta}}^{-1}{\dot{p}}_{{xy,{0}}}\\\
\dot{\theta}\end{array}\right]$
$\displaystyle=\left[\begin{array}[]{cc}{\mathrm{R}_{\theta}}&0\\\
0&1\end{array}\right]\sum_{k=1}^{N_{k}}\left[\begin{array}[]{c}\mathrm{H}_{q,k}{\dot{q}}_{{xy,{k}}}\\\
{q}_{{xy,{j}}}^{\mathsf{T}}\mathrm{S}\mathrm{H}_{q,k}{\dot{q}}_{{xy,{k}}}\end{array}\right]$
$\displaystyle\left[\begin{array}[]{c}{\mathrm{R}_{\theta}}^{-1}{\dot{p}}_{{xy,{0}}}\\\
\dot{\theta}\end{array}\right]=\left(\sum_{k=1}^{N_{k}}\left[\begin{array}[]{cc}\mathrm{H}_{q,k}&\mathrm{H}_{q,k}\mathrm{S}{q}_{{xy,{k}}}\\\
{q}_{{xy,{j}}}^{\mathsf{T}}\mathrm{S}\mathrm{H}_{q,k}&{q}_{{xy,{j}}}^{\mathsf{T}}\mathrm{S}\mathrm{H}_{q,k}\mathrm{S}{q}_{{xy,{k}}}\end{array}\right]\right)^{-1}\sum_{k=1}^{N_{k}}$
$\displaystyle\left[\begin{array}[]{c}\mathrm{H}_{q,k}\\\
{q}_{{xy,{j}}}^{\mathsf{T}}\mathrm{S}\mathrm{H}_{q,k}\end{array}\right]{\dot{q}}_{{xy,{k}}}$