# Multi-Domain Walking with Reduced-Order Models of Locomotion

Min Dai, Jaemin Lee and Aaron D. Ames  The authors are with the Department of
Mechanical and Civil Engineering, California Institute of Technology,
Pasadena, CA 91125 USA. {mdai, jaemin87<EMAIL_ADDRESS>This work is
supported by NSF NRI award 1924526 and NSF CMMI award 1923239.

###### Abstract

Drawing inspiration from human multi-domain walking, this work presents a
novel reduced-order model based framework for realizing multi-domain robotic
walking. At the core of our approach is the viewpoint that human walking can
be represented by a hybrid dynamical system, with continuous phases that are
fully-actuated, under-actuated, and over-actuated and discrete changes in
actuation type occurring with changes in contact. Leveraging this perspective,
we synthesize a multi-domain linear inverted pendulum (MLIP) model of
locomotion. Utilizing the step-to-step dynamics of the MLIP model, we
successfully demonstrate multi-domain walking behaviors on the bipedal robot
Cassie—a high degree of freedom 3D bipedal robot. Thus, we show the ability to
bridge the gap between multi-domain reduced order models and full-order multi-
contact locomotion. Additionally, our results showcase the ability of the
proposed method to achieve versatile speed-tracking performance and robust
push recovery behaviors.

##  I Introduction

The agility and versatility displayed in human locomotion have long served as
an inspiration for the study of robotic bipedal locomotion. For humans,
walking involves a sequence of distinct gait phases as shown in Fig. 1. In the
context of forward walking, these phases encompass the swing foot’s heel
strike, toe strike, the transition of weight from the new stance foot’s heel
to toe, and the subsequent heel lift and ankle push-off of the stance foot
[1]. In contrast, walking robots typically rely on flat-footed gaits. This
preference often arises from mathematical convenience: a desire for a
feedback-linearizable fully-actuated system [2, 3] or a direct application of
point-foot walking methods on robots with conventional feet [4, 5].
Nevertheless, multi-domain gait presents compelling biomechanical advantages:
the heel strike effectively dampens impact forces and ankle push-off is
remarkably energy-efficient. These advantages have been substantiated in
robotic applications as in [6, 7]. Furthermore, compared to flat-foot gait,
the capacity to raise the heel permits longer strides within the same joint
constraints, resulting in faster walking speeds [8].

Researchers have explored methods for realizing multi-domain walking on
robotic platforms due to its advantages. Full-model based methods [9, 10, 11,
7] employ multi-domain trajectory optimization within the Hybrid Zero Dynamics
(HZD) framework. They entail solving a challenging nonlinear optimization
problem, which can be computationally demanding and prone to convergence
issues for obtaining a single periodic orbit. Additionally, these methods
necessitate offline trajectory generation for different periodic orbits
associated with different speeds and contact sequences. Furthermore, they are
sensitive to model discrepancies and require a heuristic foot placement
regulator [12] for stabilization—often synthesized from reduced-order models,
i.e., the “Raibert controller” [13].

Figure 1: A complete gait cycle: (top) Human multi-domain walking. (center)
MLIP walking inspired by human walking. (bottom) Cassie’s multi-domain walking
stabilized through the MLIP model.

In the realm of prosthesis legs and feet, multi-domain walking is often
utilized but with “model-free” controllers. These controllers typically
consist of combinations of low-gain impedance control [14, 15] and torque-
based control, which is determined either using predefined parameters [16] or
through a biomechanically inspired model [17]. This flexible controller design
allows for expert tuning to tailor parameters to individual users under
different scenarios, but this tuning is time-consuming. In addition, stability
is not a primary concern in this domain, as it is assumed that the human user
can stabilize themselves through stepping. Some approaches incorporate the
human into the HZD gait generation loop [18], but they encounter similar to
those previously mentioned.

In this paper, we present a novel framework that enables multi-domain walking,
including: heel-to-toe, toe-to-heel, and flat-footed behaviors. We represent
these different phases of locomotion through a hybrid dynamical system model.
Leveraging this, and as illustrated in Fig. 1, our approach begins with the
introduction of a reduced-order model, the _multi-domain linear inverted
pendulum (MLIP)_ , specifically designed to capture the intricate weight-
shifting dynamics inherent in multi-domain walking. Subsequently, we develop a
controller that stabilizes the Poincaré map associated with the hybrid
dynamics of the MILP, i.e., ensures the stability of the step-to-step
dynamics. Via the construction of outputs for the full-order system, we
synthesize a feedback controller that realizes multi-domain walking on the
bipedal robot Cassie. This results in a remarkably human-like walking gait,
adhering to the same gait cycle time distribution observed in human walking.
Leveraging the effective ankle push-off during heel-toe-toe walking, we attain
an impressive walking speed of 2.15 m/s on Cassie, which cannot be realized
using a flat-footed gait due to inherent physical joint limits. Notably, our
method surpasses existing approaches by offering versatile walking behaviors
adaptable to diverse gait parameters and commanded speeds, all while
guaranteeing stability and eliminating the need for offline optimization.
Furthermore, our framework demonstrates robustness against external
disturbances. A collection of the resulting walking behaviors is available in
the accompanying video111https://youtu.be/8u5ZiWe_qlw.

The rest of the paper is structured as follows: Section II introduces the
hybrid control problem of multi-domain walking. Sections III and IV detail the
proposed MLIP model and its integration into the full robot model,
respectively. We then present the results and evaluate the performance under
different circumstances in section V. Finally, the paper concludes with
Section VI.

##  II Hybrid Dynamics of Bipedal Robots

Bipedal walking is represented by a hybrid control system [12, 19] defined as
the tuple $\mathcal{HC}=(\Gamma,\mathcal{D},\mathcal{S},\Delta,\mathcal{FG})$.
Each component is explained as follows:

  * •

$\Gamma=(V,E)$ is a directed circle graph with a set of vertices
$V=\\{v_{\text{i}}\\}_{\text{i}\in I}$ and a set of directed edges
$E=\\{v_{\text{i}}\rightarrow v_{\text{j}}\\}_{\text{i},\text{j}\in I}$. $I$
is an indexed set of domains that we will illustrate soon.

  * •

$\mathcal{D}=\\{\mathcal{D}_{v}\\}_{v\in V}$ is a set of domains of
admissibility.

  * •

$\mathcal{S}=\\{\mathcal{S}_{e}\\}_{e\in E}$ is a set of guards.

  * •

$\Delta=\\{\Delta_{e}\\}_{e\in E}$ is a set of reset maps

  * •

$\mathcal{FG}=\\{f_{v},g_{v}\\}_{v\in V}$ is the continuous control system,
which is a set of vector fields on the state manifolds.

In any domain, the robot’s continuous dynamics can be obtained from the Euler-
Lagrangian equations:

$\displaystyle
D(\boldsymbol{q})\ddot{\boldsymbol{q}}+H(\boldsymbol{q},\dot{\boldsymbol{q}})=B\boldsymbol{\tau}+J_{\text{i}}(\boldsymbol{q})^{T}\boldsymbol{f}_{\text{i}},$
(1) $\displaystyle
J_{\text{i}}(\boldsymbol{q})\ddot{\boldsymbol{q}}+\dot{J}_{\text{i}}(\boldsymbol{q},\dot{\boldsymbol{q}})\dot{\boldsymbol{q}}=0,$
(2)

where $\boldsymbol{q}\in Q$ is a set of generalized coordinates in the
$n$-dimensional configuration space $Q$,
$D(\boldsymbol{q})\in\mathbb{R}^{n\times n}$,
$H(\boldsymbol{q},\dot{\boldsymbol{q}})\in\mathbb{R}^{n}$,
$B\in\mathbb{R}^{n\times m}$ are the inertia matrix, the collection of
centrifugal, Coriolis and gravitational forces, and the actuation matrix,
respectively. Additionally, $\boldsymbol{\tau}\in U\subseteq\mathbb{R}^{m}$
stands for input torque, $J_{\text{i}}(\boldsymbol{q})\in\mathbb{R}^{n\times
h_{\text{i}}}$ is the domain-specific Jacobian matrix related to contact
constraints, and $\boldsymbol{f}_{\text{i}}\in\mathbb{R}^{h}_{\text{i}}$
represents the corresponding constraint wrench.

Discrete impacts are assumed to be instantaneous and plastic, with solution
derivations detailed in [19]. Denoting
$\boldsymbol{x}=[\boldsymbol{q}^{T},\dot{\boldsymbol{q}}^{T}]^{T}\in\mathcal{TQ}$,
the equation of motion for the hybrid system is as follows:

$\displaystyle\begin{cases}\dot{\boldsymbol{x}}&=f_{v}(\boldsymbol{x})+g_{v}(\boldsymbol{x})\boldsymbol{\tau}\quad\boldsymbol{x}\in\mathcal{D}_{v}\setminus\mathcal{S}_{e}\\\
\boldsymbol{x}^{+}&=\Delta_{e}(\boldsymbol{x}^{-})\hskip
39.83368pt\boldsymbol{x}^{-}\in\mathcal{S}_{e}\end{cases},$ (3)

for all $v\in V$ and corresponding $e\in E$.

Inspired by human heel-to-toe walking, our hybrid system model incorporates
three domains: fully-actuated (FA), under-actuated (UA) and over-actuated
(OA). This yields $V={v_{\text{FA}},v_{\text{UA}},v_{\text{OA}}}$ and
$E={v_{\text{FA}}\rightarrow v_{\text{UA}},v_{\text{UA}}\rightarrow
v_{\text{OA}},v_{\text{OA}}\rightarrow v_{\text{FA}}}$, as illustrated in Fig.
2.

######  Definition 1.

We define the domains for a hybrid dynamical system of bipedal robots as
follows:
(FA) Fully-actuated domain: This describes the single support phase where the
stance foot has full ground contact.
(UA) Under-actuated domain: This characterizes the single support phase where
only the stance toe has ground contact.
(OA) Over-actuated domain: This corresponds to the double support phase, with
both the back leg toe and the front leg heel in contact with the ground.

######  Definition 2.

Guards are defined to represent the switching condition from one domain to the
next following defined edges. More specifically, the transition from the FA to
UA domain occurs when the stance heel lifts off the ground. The guard is given
by:

$\displaystyle\mathcal{S}_{\text{FA}\rightarrow\text{UA}}$
$\displaystyle\coloneqq\\{(\boldsymbol{x},\boldsymbol{\tau})\in\mathcal{TQ}\times
U:f^{z}_{\text{st,heel}}(\boldsymbol{x},\boldsymbol{\tau})=0\\}.$

In addition, the UA to OA domain transition takes place when the swing heel
strikes the ground, as indicated by:

$\displaystyle\mathcal{S}_{\text{UA}\rightarrow\text{OA}}$
$\displaystyle\coloneqq\\{\boldsymbol{x}\in\mathcal{TQ}:z_{\text{sw,heel}}(\boldsymbol{x})=0\\}.$

Lastly, transitioning from the OA to FA domain occurs when the swing toe
contacts the ground, and the previous stance toe lifts off, providing the
following guard condition:

$\displaystyle\mathcal{S}_{\text{OA}\rightarrow\text{FA}}\coloneqq\\{(\boldsymbol{x},\boldsymbol{\tau})\in\mathcal{TQ}\times
U:$ $\displaystyle z_{\text{sw,toe}}(\boldsymbol{x})=0,$ $\displaystyle
f^{z}_{\text{st,toe}}(\boldsymbol{x},\boldsymbol{\tau})=0\\}.$

Figure 2: Hybrid system directed graph for heel-to-toe walking

######  Remark 1.

To conserve space, we have exclusively presented the heel-to-toe walking mode.
One can equivalently define $\mathcal{HC}$ models for toe-to-heel walking and
flat-footed walking. In flat-footed walking, the UA domain corresponds to the
deactivation of the stance ankle torque.

##  III Step-to-Step Dynamics and Stabilization for Multi-Domain LIP Model

In this section, we first propose a reduced-order model, termed the multi-
domain linear inverted pendulum (MLIP) model, which is a variant of the
canonical linear inverted pendulum (LIP) model [2] that can describe multi-
domain walking. This extension involves incorporating the position of the
zero-moment point (ZMP) [20] as an additional state variable. After
characterizing its step-to-step (S2S) dynamics, we apply a linear controller
to stabilize the error dynamics, considering the discrepancy between the
actual dynamics of the robot and the reduced-order model.

Figure 3: MLIP gait cycle for heel-to-toe walking. Note that as the MLIP model
has massless legs and feet, the foot angles do not impact dynamics. The stance
and swing foot pitch angles shown are up to the user’s choice. The feet drawn
are straight lines, but curved feet with an arc length $\rho$ would result in
the same dynamics. Figure 4: MLIP gait cycle for (a) toe-to-heel and (b) flat-
footed walking with arrows indicating walking direction.

###  III-A MLIP Model and Step-to-step Dynamics

As shown in Fig. 3 and 4, the MLIP model includes a point mass and two
massless telescopic legs. It also has a constant center of mass (CoM) height
$z_{0}$ relative to the stance pivot as in the LIP model. What sets it apart
from the conventional LIP model is the inclusion of a pair of feet with a
known arc length denoted as $\rho$, which can be calculated from foot
curvature. In accordance with the bipedal locomotion domains outlined in Sec.
II, the MLIP model encompasses the UA, OA, and FA domains. The continuous
dynamics of the MLIP model in all domains are governed by the following linear
equations:

$\displaystyle\frac{d}{dt}\underbrace{\begin{bmatrix}p\\\ L\\\
p_{\text{zmp}}\end{bmatrix}}_{\boldsymbol{\xi}}=\underbrace{\begin{bmatrix}0&\frac{1}{z}&0\\\
g&0&-g\\\ 0&0&0\end{bmatrix}}_{A_{\text{ct}}}\begin{bmatrix}p\\\ L\\\
p_{\text{zmp}}\end{bmatrix}+\underbrace{\begin{bmatrix}0\\\ 0\\\
1\end{bmatrix}}_{B_{\text{ct}}}\dot{p}_{\text{zmp}},$ (4)

where $p$, $L$, and $p_{\text{zmp}}$ represent the horizontal CoM position,
mass-normalized centroidal angular momentum [21], and horizontal ZMP position,
all defined relative to the stance pivot. The stance pivot we refer to is the
stance pivot at the UA phase thus it is dependent on walking mode. When
walking in a heel-to-toe manner, the stance pivot is positioned at the toe,
while in toe-to-heel walking, it is located at the heel. For flat-footed
walking, the stance pivot can be picked anywhere from heel to toe. In this
context, we choose the pivot point to be directly under the ankle. Given that
the dynamics are linear, the end-of-domain states
$\boldsymbol{\xi}_{\text{i}}^{-}$ has closed-form solution given by:

$\displaystyle\boldsymbol{\xi}_{\text{i}}^{-}=\underbrace{e^{A_{\text{ct}}T_{\text{i}}}}_{A_{\text{i}}}\boldsymbol{\xi}_{\text{i}}^{+}+\int_{0}^{T_{\text{i}}}e^{A_{\text{ct}}(T_{\text{i}}-t)}B_{\text{ct}}\dot{p}_{\text{zmp,
i}}(t)dt,$ (5)

where superscripts $(\cdot)^{+/-}$ indicate the beginning and end of each
domain, respectively, and $T_{\text{i}}$ represents the time duration of the
i-th domain. Unlike the guards defined in Definition 2, the transitions
between domains in the MLIP model are purely time-based, as the legs are
effectively virtual. To determine $\dot{p}_{\text{zmp, i}}$ in each domain, we
draw inspiration from human walking data [22], as elaborated below.

FA:  We denote the distance that the ZMP travels during the fully-actuated
phase as $l$. In the context of heel-to-toe walking, this distance corresponds
to $\rho$. For toe-to-heel and flat-footed walking, the values are $-\rho$ and
0, respectively.

$\displaystyle\dot{p}_{\text{zmp, FA}}(t)=\frac{l}{T_{\text{FA}}}.$ (6)

UA:  During the under-actuated phase, the ZMP is considered to be fixed at the
stance pivot. Thus, $\dot{p}_{\text{zmp, UA}}(t)=0$.

OA:  During the over-actuated phase, the ZMP shifts from the previous stance
pivot to the new stance leg, thus traveling the step length denoted as $u$
within the step duration $T_{\text{OA}}$:

$\displaystyle\dot{p}_{\text{zmp, OA}}(t)=\frac{u}{T_{\text{OA}}}.$ (7)

As depicted in Fig. 3, in the context of heel-to-toe walking, the parameter
$u$ indicates the distance between the stance toe and the swing heel. However,
in the case of toe-to-heel and flat-footed walking, $u$ is adapted to denote
the distance between the stance heel and the swing toe, as well as the
separation between the stance under the ankle and the swing under the ankle.

Given $\dot{p}_{\text{zmp, i}}$ is defined to be independent of time for all
domains, we further simplified the above dynamics (5) :

$\displaystyle\int_{0}^{T_{\text{i}}}e^{A_{\text{ct}}(T_{\text{i}}-t)}B_{\text{ct}}\dot{p}_{\text{zmp,
i}}(t)dt$
$\displaystyle=\int_{0}^{T_{\text{i}}}e^{A_{\text{ct}}(T_{\text{i}}-t)}dtB_{\text{ct}}\dot{p}_{\text{zmp,
i}}$
$\displaystyle=\underbrace{\int_{0}^{T_{\text{i}}}e^{A_{\text{ct}}(T_{\text{i}}-t)}dtB_{\text{ct}}\frac{1}{T_{\text{i}}}}_{B_{\text{i}}}d_{\text{i}}.$

As a result, we can express the solution to the continuous time dynamics for
$i$-th domain as follows:

$\displaystyle\boldsymbol{\xi}_{\text{i}}^{-}=A_{\text{i}}\boldsymbol{\xi}_{\text{i}}^{+}+B_{\text{i}}d_{\text{i}}$
(MLIP-CT)

where $d_{\text{OA}}=u$, $d_{\text{FA}}=l$, and $d_{\text{UA}}=0$.

With continuous phase dynamics defined, we need to specify the impact
dynamics. For the model with massless legs, discrete state jumps resulting
from impact do not occur. Instead, the impact equation characterizes the
effects of switching between the stance and swing legs, which is defined to
happen at the transition from OA to FA domain. Thus, the impact equations are
given by:

$\displaystyle\begin{cases}\boldsymbol{\xi}_{\text{OA}}^{+}&=\boldsymbol{\xi}_{\text{UA}}^{-}\\\
\boldsymbol{\xi}_{\text{FA}}^{+}&=\boldsymbol{\xi}_{\text{OA}}^{-}+B_{\Delta}u+C_{\Delta}\\\
\boldsymbol{\xi}_{\text{UA}}^{+}&=\boldsymbol{\xi}_{\text{FA}}^{-}\end{cases},$
(MLIP-DT)

where $B_{\Delta}=\begin{bmatrix}-1&0&-1\end{bmatrix}^{T}$ and
$C_{\Delta}=\begin{bmatrix}-l&0&-l\end{bmatrix}^{T}$.

We are interested in understanding the S2S dynamics for stabilizing the
system. Consider the pre-impact state at the UA phase as $\boldsymbol{\xi}[k]$
in the $k$-th step, which is shown in red in Fig. 3 and 4, a complete step
evolution follows this sequence:
$\boldsymbol{\xi}[k]\rightarrow\boldsymbol{\xi}_{\text{OA}}^{+}\rightarrow\boldsymbol{\xi}_{\text{OA}}^{-}\rightarrow\boldsymbol{\xi}_{\text{FA}}^{+}\rightarrow\boldsymbol{\xi}_{\text{FA}}^{-}\rightarrow\boldsymbol{\xi}_{\text{UA}}^{+}\rightarrow\boldsymbol{\xi}[k+1]$.
Consequently, the S2S dynamics of the MLIP model can be described using Eq.
(MLIP-DT) and (MLIP-CT):

$\displaystyle\boldsymbol{\xi}_{k+1}=A_{\boldsymbol{\xi}}\boldsymbol{\xi}_{k}+B_{\boldsymbol{\xi}}u_{k}+C_{\boldsymbol{\xi}}$
(8)

where the detailed matrices are obtained as follows:

$\displaystyle A_{\boldsymbol{\xi}}$
$\displaystyle=A_{\text{UA}}A_{\text{FA}}A_{\text{OA}},$ (9a) $\displaystyle
B_{\boldsymbol{\xi}}$
$\displaystyle=A_{\text{UA}}A_{\text{FA}}(B_{\text{OA}}+B_{\Delta}),$ (9b)
$\displaystyle C_{\boldsymbol{\xi}}$
$\displaystyle=A_{\text{UA}}A_{\text{FA}}(B_{\text{FA}}l+C_{\Delta}).$ (9c)

Given that the Poincaré section is defined at the end of the UA phase, it
follows that $p_{\text{zmp,k}}=0$ for all $k\in\mathbb{N}$. Additionally, it
can be confirmed that $A_{\boldsymbol{\xi}}(3,3)=1$,
$B_{\boldsymbol{\xi}}(3)=0$, and $C_{\boldsymbol{\xi}}(3)=0$. Therefore, we
use $\boldsymbol{x}^{M}\in\mathbb{R}^{2}$ to represent the horizontal CoM
states, specifically $\boldsymbol{x}^{M}=[p,L]^{T}$, such that
$\boldsymbol{\xi}=[\boldsymbol{x}^{M};0]$. Using the two-dimensional states,
we can express the S2S dynamics of the MLIP as follows:

$\displaystyle\boldsymbol{x}_{k+1}^{\text{M}}=A^{\text{M}}\boldsymbol{x}_{k}^{\text{M}}+B^{\text{M}}u_{k}+C^{\text{M}},$
(MLIP-S2S)

where $A^{\text{M}}=A_{\boldsymbol{\xi}}(1:2,1:2)$,
$B^{\text{M}}=B_{\boldsymbol{\xi}}(1:2)$, and
$C^{\text{M}}=C_{\boldsymbol{\xi}}(1:2)$.

######  Remark 2.

It might appear that the MLIP model requires non-trivial FA, UA, and OA phases
to obtain well-defined S2S dynamics, given the presence of
$\frac{1}{T_{\text{i}}}$ in the definition of $B_{\text{i}}$. However, it’s
important to note that the integral from 0 to $T_{\text{i}}$ always results in
zero when $T_{\text{i}}=0$. In cases of $T_{\text{OA}}=0$, we also need to set
$B_{\Delta}=\begin{bmatrix}-1&0&0\end{bmatrix}^{T}$ to immediately transfer
the zero moment point to the new stance foot.

######  Remark 3.

The second coordinate $L$ can be replaced by linear velocity, as done in [2,
5]. However, we have chosen to use angular momentum about stance pivot for
higher data quality, given a Kalman filter for angular momentum can be easily
set up as demonstrated in [23]. Replacing the second coordinate will result in
a modified $A_{\text{ct}}$, while the linear structure of the step-to-step
dynamics remains unchanged.

######  Remark 4.

One can retrieve the Hybrid-LIP model in [5] by using the CoM linear velocity
as the second coordinate and setting $l=0$ and
$T_{\text{FA}}=T_{\text{OA}}=0$. However, there are different assumptions
regarding the motion of the ZMP in the OA phase, which is equivalent to the
double support phase in the Hybrid-LIP model. For more in-depth information,
one can refer to [5].

###  III-B Stabilization for Robot Dynamics

Having formulated the S2S dynamics for the MLIP model as a discrete-time
linear control system, we now aim to stabilize robot walking using it. In
practice, we encounter a significant challenge in obtaining the S2S dynamics
of the robot due to its inherently nonlinear nature. However, if we assume
that the robot’s control scheme ensures the availability of the next step, we
can mathematically express its S2S dynamics as the Poincaré map:

$\displaystyle\boldsymbol{x}_{k+1}^{-}=\mathcal{P}_{\boldsymbol{x}}(\boldsymbol{x}_{k}^{-},\tau(t)).$
(10)

Here, $\boldsymbol{x}_{k}^{-}\in\mathcal{TQ}$ is the robot state at the end of
the UA phase, i.e., pre-impact states. The associated Poincaré map is denoted
as $\mathcal{P}_{\boldsymbol{x}}$. Our focus lies on the pre-impact CoM
states, which can be denoted
$\boldsymbol{x}^{R}=[p^{\text{R},-},L^{\text{R},-}]^{T}\in\mathbb{R}^{2}$. We
can express the evolution of these pre-impact CoM states as

$\displaystyle\boldsymbol{x}_{k+1}^{\text{R}}=\mathcal{P}_{\boldsymbol{x}^{\text{R}}}(\boldsymbol{x}_{k}^{\text{R}},\tau(t)).$
(11)

Using the MLIP S2S dynamics as an approximation, the robot S2S dynamics can be
expressed as:

$\displaystyle\boldsymbol{x}_{k+1}^{\text{R}}=A^{\text{M}}\boldsymbol{x}_{k}^{\text{R}}+B^{\text{M}}u_{k}^{\text{R}}+C^{\text{M}}+w,$
(12)

where $u_{k}^{\text{R}}$ denotes the $k$-th step size of the robot.
$w=\mathcal{P}_{\boldsymbol{x}^{\text{R}}}(\boldsymbol{x}_{k}^{\text{R}},\tau(t))-A^{\text{M}}\boldsymbol{x}_{k}^{\text{R}}-B^{\text{M}}u_{k}^{\text{R}}-C^{\text{M}}$
represents the model discrepancy, i.e. the integrated dynamics difference
between the robot and the MLIP over a step. It is assumed that the realizable
set of walking behaviors satisfies $w\in\mathbf{W}$, where $\mathbf{W}$ is a
bounded set as mentioned in [5]. Let
$\boldsymbol{e}:=\boldsymbol{x}^{\text{R}}-\boldsymbol{x}^{\text{M}}$
represent the error state, a stabilizing controller can be designed as
follows:

$\displaystyle
u^{\text{R}}_{k}=u^{\text{M}}_{k}+K(\boldsymbol{x}^{\text{R}}_{k}-\boldsymbol{x}^{\text{M}}_{k}),$
(13)

which yields the error dynamics:

$\displaystyle\boldsymbol{e}_{k+1}=(A^{\text{M}}+B^{\text{M}}K)\boldsymbol{e}_{k}+w,$
(14)

where $K$ is the controller gain. From linear control theory, any selections
of $K$ that result in stable $A^{\text{M}}+B^{\text{M}}K$ can drive
$\boldsymbol{e}$ to converge to an invariant set $\mathbf{E}$ as in [24], i.e.
if $\boldsymbol{e}_{k}\in\mathbf{E}$, $\boldsymbol{e}_{k+1}\in\mathbf{E}$. In
this work, the controller gain $K$ was determined using the linear quadratic
regulator. It’s important to note that this process is independent of the
specific model, as long as it can be represented as a linear discrete-time
system. Similar techniques have been successfully employed in related works,
such as [6] for S2S dynamics derived from numerically linearized full robot
dynamics and in [5] for a different reduced-order model.

###  III-C MLIP Periodic Orbits

In the context of walking, the robot is often required to follow a desired
velocity. A straightforward approach is to employ closed-form reduced-order
model periodic orbits, i.e. set $\boldsymbol{x}^{\text{M}}=\boldsymbol{x}^{*}$
to be the desired state. As MLIP is a planar model, it allows for decoupled
planning of sagittal and lateral motion. We thus present the results for both
Period-1 and Period-2 orbit that are suitable for sagittal and lateral
planning, respectively. A visualization of the periodic orbits using different
parameters is shown in Fig. 5.

Period-1 Orbit: The desired step size $u^{*}$ for a Period-1 orbit is
determined by the desired walking velocity $v^{d}$ and the step duration $T$,
where $u^{*}=v^{d}T$. The corresponding desired periodic pre-impact state for
achieving $u^{*}$ is calculated by setting
$\boldsymbol{x}_{k+1}=\boldsymbol{x}_{k}=\boldsymbol{x}^{*}$ in Eq.
(MLIP-S2S):

$\displaystyle\boldsymbol{x}^{*}=(I_{2\times
2}-A^{\text{M}})^{-1}(B^{\text{M}}u^{*}+C^{\text{M}}),$ (15)

where $I_{2\times 2}\in\mathbb{R}^{2\times 2}$ is the identity matrix. In this
context, the controller can be expressed as:

$\displaystyle
u^{\text{R}}=u^{*}+K(\boldsymbol{x}^{\text{R}}-\boldsymbol{x}^{*}).$ (16)

Period-2 Orbit: Unlike P1 orbits, there is no unique solution for the P2 orbit
that achieves a desired velocity $v^{d}$. We use subscripts L/R to denote the
left or right stance leg. The step sizes must satisfy the equation
$u_{\text{L}}^{*}+u^{*}_{\text{R}}=2v^{d}T$, and the choice of one step size
determines the orbit. Solving $\boldsymbol{x}_{k+2}=\boldsymbol{x}_{k}$ yields
the desired periodic pre-impact states:

$\displaystyle\boldsymbol{x}_{\text{L/R}}^{*}=(I_{2\times
2}-(A^{\text{M}})^{2})^{-1}($ $\displaystyle
A^{\text{M}}B^{\text{M}}u_{\text{L/R}}^{*}+B^{\text{M}}u^{*}_{\text{R/L}}$
$\displaystyle+A^{\text{M}}C^{\text{M}}+C^{\text{M}}).$

Consequently, the controller can be written as
$u^{\text{R}}_{\text{L/R}}=u^{*}_{\text{L/R}}+K(\boldsymbol{x}^{\text{R}}-\boldsymbol{x}^{*}_{\text{L/R}})$.
When applying P2 orbits in the lateral plane, the state vector
$\boldsymbol{x}$ is defined as $\boldsymbol{x}:=[p_{y},-L_{x}]^{T}$ to ensure
consistency with the sign conventions used in the $A_{\text{ct}}$ matrix.

Figure 5: Phase portraits depicting various periodic orbits with $z_{0}=0.8$m:
(a) Heel-to-toe walking at speeds of 0, 0.5, 1, 2 m/s shown in blue, red,
yellow, and purple lines. (b) Comparison of heel-to-toe walking (red) and
flat-footed walking (black) at 2m/s, highlighting FA, UA, and OA phases with
solid, dashed, and dotted lines, respectively. (c) Toe-to-heel walking at -1
m/s, showcasing different total step times (T = 0.4s, 0.6s, and 0.8s) using
blue, red, and yellow lines. (d) Flat-footed P2 orbit at 0m/s with nominal
step widths of 0.3m, 0.4m, and 0.5m in blue, red, and yellow lines.

##  IV Robot Implementation

The reduced-order model generates discrete commands using step-to-step
stabilization. However, to translate these commands into real-time control
signals for the physical robot, it is imperative to construct continuous
control outputs and design corresponding feedback controllers. This process is
detailed in this section.

###  IV-A Output Definition

Figure 6: Robot Cassie schematics and output definition.

We have implemented our proposed method on robot Cassie, a 3D underactuated
robot developed by Agility Robotics [25]. As depicted in Fig. 6, each leg of
Cassie is modeled with 6 degrees of freedom (DOF), comprising five motored
joints and one passive tarsus joint. When combined with the 6 DOF for the
floating-base pelvis frame, the total DOF for the robot is 18.

To realize the MLIP-based approach on Cassie, the output design adheres to
specific requirements. The vertical CoM position $z_{\text{com}}$ relative to
the stance pivot should remain approximately constant. The vertical position
of the swing foot $z_{\text{st foot}}$ is constructed to periodically lift off
and strike the ground. The horizontal position of the swing foot
$\\{x,y\\}_{\text{sw foot}}$ relative to the stance pivot is controlled to
achieve the desired step size $\\{u_{x},u_{y}\\}$ from the MLIP-based stepping
controller Eq. (16). The pitch angles of the stance and swing foot should be
controlled to provide the desired contact location corresponding to the
desired walking mode. Additionally, the pelvis roll and pitch angles
$\\{\theta^{y},\theta^{x}\\}_{\text{pelvis}}$ and the stance and swing hip yaw
angles $\theta^{z}_{\text{st hip}}$, $\theta^{z}_{\text{sw hip}}$ should be
controlled to fully constrain the walking behaviors. With these considerations
in mind, the desired walking behavior is encoded by the virtual constraints
[19], defined as:

$\displaystyle\mathbfcal{Y}=\boldsymbol{h}^{a}-\boldsymbol{h}^{d}\in\mathbb{R}^{12},$

where $\boldsymbol{h}^{a}$ and $\boldsymbol{h}^{d}$ denote the actual and
desired output defined as follows:

$\begin{split}\boldsymbol{h}=\textrm{col}(&\\{x,y,z\\}_{\text{com}},\theta^{z}_{\text{st
hip}},\theta^{y}_{\text{st
foot}},\\{\theta^{y},\theta^{x}\\}_{\text{pelvis}},\\\ &\\{x,y,z\\}_{\text{sw
foot}},\theta^{z}_{\text{sw hip}},\theta^{y}_{\text{sw
foot}})\in\mathbb{R}^{12}.\end{split}$

In Fig. 6, we provide a visualization of these outputs, specifically for flat-
footed walking. For heel-to-toe walking, the position of CoM is defined to be
relative to the stance toe, and the position of the swing foot is represented
as the vector from the stance toe to the swing heel. Similar modifications
need to be applied to toe-to-heel walking.

Cassie employs a line foot design, which means that only the toe pitch motor
is present on the foot link. Consequently, the robot lacks motor control for
adjusting the foot roll angle. This limitation leads to inherent
underactuation in the lateral plane during single support phases. Similarly,
in the OA domain, where only the back toe and front heel make contact with the
ground for heel-to-toe walking, the ZMP can solely reside along the line
connecting these two contact points. This configuration results in coupled
control of the xCoM and yCoM components. Therefore, achieving independent
control of both xCoM and yCoM is not feasible during OA phases. These hardware
constraints lead to the following choice of selection matrices $S_{\text{i}}$
for each domain such that
$\boldsymbol{h}_{{\text{FA},\text{UA},\text{OA}}}=S_{{\text{FA},\text{UA},\text{OA}}}\boldsymbol{h}$:

$\displaystyle S_{\text{FA}}$
$\displaystyle=\textrm{diag}(\\{1,0,1\\},1,0,\\{1,1\\},\\{1,1,1\\},1,1),$
(17a) $\displaystyle S_{\text{UA}}$
$\displaystyle=\textrm{diag}(\\{0,0,1\\},1,1,\\{1,1\\},\\{1,1,1\\},1,1),$
(17b) $\displaystyle S_{\text{OA}}$
$\displaystyle=\textrm{diag}(\\{0,0,1\\},1,1,\\{1,1\\},\\{0,0,0\\},0,1).$
(17c)

The selection matrix for each domain aligns with the holonomic constraints
present in that specific domain. In the UA and OA phases, we assume patch
contact with the ground, meaning that the contact points’ positions x, y, z
and yaw angle are constrained. During the FA phase, we apply the line contact
assumption, which adds the constraint on the contact pitch angle. In the
current formulation, we allow the horizontal CoM states to evolve passively
during the OA phase. This choice helps avoid mismatched commands for xCoM and
yCoM control during that phase.

The construction of the most desired outputs for our control framework is
consistent with the approach outlined in [5]. However, there are exceptions in
the case of the horizontal COM position, stance foot pitch, and swing foot
pitch angles.

####  IV-A1 Horizontal COM position

During the FA phase, one option is to control the ZMP to transit from the heel
to the toe, as assumed in the MLIP model construction. However, this approach
requires feedback control of the ZMP position at the jerk level, which is
impractical. Instead, we opt to directly regulate the horizontal COM states at
the end of the FA phase to reach a desired state denoted as
$\boldsymbol{x}^{\text{M},*}_{\text{FA}^{-}}$. This state can be calculated
using the MLIP step-to-step dynamics, with the Poincaré section instead chosen
to be the end of the FA phase, i.e., FA-. The trajectory for this control is
defined using a Bézier polynomial that can be written as:

$\displaystyle x_{\text{com}}^{d}(s_{\text{FA}})\coloneqq
b_{\text{xcom}}(s_{\text{FA}})=A(s_{\text{FA}})\alpha_{\text{xcom}},$ (18)

where $s_{\text{FA}}\coloneqq\frac{t_{\text{FA}}}{T_{\text{FA}}}\in[0,1)$ is
the phase variable in FA domain, and $A(s_{\text{FA}})\in\mathbb{R}^{1\times
n_{b}}$ is defined using the definition of Bézier polynomial. The $k$-th
element of $A(s_{\text{FA}})$ is computed as:

$A_{k}(s_{\text{FA}})=\frac{n_{b}!}{k!(n_{b}-k)!}s_{\text{FA}}^{k}(1-s_{\text{FA}})^{n_{b}-k}.$

Additionally, $\alpha_{\text{xcom}}\in\mathbb{R}^{n_{b}}$ represents the
coefficients of a Bézier polynomial of degree $n_{b}$. These coefficients are
determined at the start of the FA phase, subject to the following linear
equality constraints:

$\displaystyle\begin{bmatrix}A(0)\\\ \dot{A}(0,T_{\text{FA}})\\\ A(1)\\\
\dot{A}(1,T_{\text{FA}})\end{bmatrix}\alpha_{\text{xcom}}=\begin{bmatrix}x_{\text{com,
FA+}}^{a}\\\ \dot{x}_{\text{com, FA+}}^{a}\\\
p^{\text{M},*}_{\text{FA}^{-}}\\\
\frac{1}{z_{0}}L^{\text{M},*}_{\text{FA}^{-}}\end{bmatrix}.$ (19)

Here, $\dot{A}(s_{\text{FA}},T_{\text{FA}})\in\mathbb{R}^{1\times n_{b}}$ is
also determined based on the definition of the Bézier polynomial, ensuring
that
$\dot{b}_{\text{xcom}}(s_{\text{FA}},T_{\text{FA}})=\dot{A}(s_{\text{FA}},T_{\text{FA}})\alpha_{\text{xcom}}$.
The actual horizontal CoM position and velocity at the beginning of the domain
is given by $x_{\text{com, FA+}}^{a}$ and $\dot{x}_{\text{com, FA+}}^{a}$.
Notably, the second coordinate in
$\boldsymbol{x}^{\text{M},*}_{\text{FA}^{-}}$ requires conversion from angular
momentum to velocity, as indicated.

####  IV-A2 Stance and swing foot pitch angle

The desired end-of-domain stance and swing foot pitch angles vary for
different contact modes in different domains. For instance, in the UA phase of
heel-to-toe walking, the final stance foot angle should be greater than zero
to enable heel-lift motion. However, these angles share a common output
structure:

$\displaystyle\\{\theta^{y}_{\text{sw/st foot}}\\}^{d}(s_{\text{i}})$
$\displaystyle\coloneqq(1-b_{\text{foot}}(s_{\text{i}}))\\{\theta^{y}_{\text{sw/st
foot}}\\}^{a}_{\text{i+}}$ (20)
$\displaystyle+b_{\text{foot}}(s_{\text{i}})\\{\theta^{y}_{\text{sw/st
foot}}\\}^{*}_{\text{i}},$

where $s_{\text{i}}\coloneqq\frac{t_{\text{i}}}{T_{\text{i}}}\in[0,1)$ is the
phase variable within each domain, $\\{\theta^{y}_{\text{sw/st
foot}}\\}^{a}_{\text{i+}}$ is the actual foot pitch angle at beginning of i-th
domain, $\\{\theta^{y}_{\text{sw/st foot}}\\}^{*}_{\text{i}}$ denotes the
desired end-of-domain pitch angle, $b_{\text{foot}}(s_{\text{i}})$ represents
a Bézier polynomial that transitions from 0 to 1.

###  IV-B Feedback Controller

Using the synthesized outputs, we employ a task-space quadratic programming
(QP) based controller [26] to ensure the tracking of desired trajectories
while respecting the constrained dynamics, physical motor torque limits, and
ground contact forces constraints. In each domain i and at each control loop,
we formulate the QP with optimization variables
$\boldsymbol{\ddot{q}},\boldsymbol{\tau},\boldsymbol{f}_{\text{i}}$ as
follows:

$\displaystyle\underset{\boldsymbol{\ddot{q}},\boldsymbol{\tau},\boldsymbol{f}_{\text{i}}}{\text{min}}$
$\displaystyle\quad||\ddot{\boldsymbol{h}}^{a}_{\text{i}}(q,\dot{q},\ddot{q})-\ddot{\boldsymbol{h}}^{d}_{\text{i}}-\ddot{\mathbfcal{Y}}^{t}_{\text{i}}||^{2}_{Q},$
(TSC-QP) s.t.
$\displaystyle\quad\text{Eq.}~{}\eqref{eq::eom},\eqref{eq::hol},$ (Dynamics)
$\displaystyle\quad
A_{\text{GRF}}\boldsymbol{f}_{\text{i}}\leq\boldsymbol{b}_{\text{GRF}},$
(Contact)
$\displaystyle\quad\boldsymbol{\tau}_{lb}\leq\boldsymbol{\tau}\leq\boldsymbol{\tau}_{ub}.$
(Torque Limit)

Here, $Q$ denotes a weight matrix, and
$\ddot{\mathbfcal{Y}}^{t}=-K_{p}\mathbfcal{Y}-K_{d}\dot{\mathbfcal{Y}}$ is the
target acceleration of the output that enables exponential tracking, where
$K_{p},K_{d}$ are the proportional and derivative gains. The affine contact
constraint on $\boldsymbol{f}_{\text{i}}$ approximates the contact friction
cone constraint. $\boldsymbol{\tau}_{lb}$ and $\boldsymbol{\tau}_{ub}$
represent the lower and upper torque bounds. Solving this QP yields the
optimal torque $\boldsymbol{\tau}$ that is applied to the robot.

##  V Results

We evaluate the proposed approach using our C++ implementation in the open-
sourced simulator [27] with the Mujoco physics engine [28] on the robot
Cassie. The output construction and corresponding low-level controller, as
described in (TSC-QP), are executed at a rate of 1kHz for real-time control. A
visual representation of the results is available in the supplementary video
provided earlier.

For our MLIP planning, we assume a constant CoM height of $z_{0}=0.8$m and a
foot length of $\rho=0.16$m, consistent with Cassie’s physical foot length.
Since lateral dynamics are always underactuated during the single support
phase, we use $T_{\text{SS}}$ to denote the step duration for the single
support phase, where $T_{\text{SS}}=T_{\text{FA}}+T_{\text{UA}}$ for both the
sagittal and lateral planes. In all tests, we use $T_{\text{SS}}=0.4$s and
$T_{\text{OA}}=0.1$s. For the MLIP model in the sagittal plane, we set
$T_{\text{FA}}=T_{\text{UA}}=\frac{T_{\text{SS}}}{2}$, resulting in a phase
distribution of 40% FA, 40% UA, and 20% OA, which aligns with actual human
walking data [29]. In the lateral plane, we have $T_{\text{FA}}=0$ and
$T_{\text{UA}}=T_{\text{SS}}$.

Figure 7: (a) Heel-to-toe walking at a speed of 2 m/s. (b) Toe-to-heel walking
at -1.5m/s. (c) Phase portrait depicting 2m/s steady-state heel-to-toe
walking. The yellow, red, and blue lines represent the FA, UA, and OA phases
for the robot, following the color code from Fig. 3. The dashed, solid, and
dotted lines correspond to the FA, UA, and OA phases for the periodic orbit in
the MLIP model. (d) Phase portrait illustrating steady-state lateral walking,
also using the same color code. Figure 8: Velocity tracking performance and
corresponding steady-state phase portrait for a range of walking speeds: 2
m/s, 1 m/s, 0.5 m/s, 0 m/s, -0.75 m/s, and -1.5 m/s, depicted by blue, yellow,
green, burgundy, purple, and red lines, respectively. The commanded velocity
is given in solid black lines.

###  V-A Versatile Walking

The versatility of the proposed approach enables us to showcase a wide range
of walking behaviors using Cassie. In all test scenarios, the robot initiates
from a stationary posture and is initially commanded to step-in-place for a
duration of five seconds. Subsequently, the reference velocity gradually
changes to the desired value.

Fig. 7(a) and (b) provide a glimpse of the robot’s steady-state walking under
speeds of 2 m/s and -1.5 m/s, respectively. Employing the multi-domain
formulation, Cassie exhibits a natural walking motion characterized by human-
like foot rolls with distinct UA, OA, and FA phases. In Fig. 7(c), we offer a
comparison of the phase portrait for the periodic orbit realized by Cassie and
the periodic orbit defined for MLIP for 2m/s heel-to-toe walking. Fig. 7(d)
depicts the phase portrait for lateral walking using a P2 orbit.

In Fig. 8, we examine the performance of our approach across a range of
reference velocities. Remarkably, at all reference velocities, the walking
stabilizes near the reference velocity, and the robot’s states remain within a
small, bounded error set relative to the nominal MLIP states as shown in the
phase portrait.

######  Remark 5.

Our primary objective in this work is not to achieve perfect velocity tracking
performance. As discussed in Sec. III-B, our proposed reduced-order model
planner is designed to stabilize the robot within a small bounded error set
around the planned reduced-order model trajectory. To improve the global
tracking performance, there are two viable strategies. One approach involves
integrating a high-level planner to adjust the desired velocity sent to the
robot. Alternatively, data-driven techniques [30] can be employed to reduce
robot-dependent model discrepancies, thereby reducing the size of the error
set.

###  V-B Push Recovery

Figure 9: (top) Cassie’s recovery from a 50N push during 1m/s walking.
(bottom) Disturbance force profile and corresponding velocity tracking results
for push recovery tests at walking speeds of 1m/s, 0.5m/s, and 0.75m/s,
represented by blue, red, and yellow lines respectively.

The proposed method uses online foot placement planning to stabilize the
robot. To evaluate the robustness of this approach, we apply unknown
disturbance forces to Cassie. Specifically, as shown in Fig. 9(b), we applied
a 50N and -50N push to the pelvis of the robot at times 15s and 20s, where
each lasts for 0.5s.

Fig. 9(a) shows the robot’s response to the 50N push during heel-to-toe
walking at a speed of 1m/s. Evidently, the robot quickly recovers from the
perturbation by taking a few larger steps. Fig. 9 (c) displays the CoM
velocity profiles for three different commanded speeds when subjected to the
same disturbance. Notably, all commanded speeds exhibit the ability to
withstand external forces and subsequently resume normal walking behavior.

###  V-C Maximum Speed

Figure 10: (a&b) Cassie’s walking behavior under maximum attainable speed for
heel-to-toe walking and flat-footed walking. (c&d) Velocity profiles and phase
portraits for the two scenarios, where blue and red lines are used for heel-
to-toe walking and flat-footed walking respectively. For flat-footed walking,
the robot oscillates around its steady-state orbits and eventually falls.

Multi-domain walking offers the potential for achieving larger footsteps
compared to flat-footed walking, ultimately leading to higher walking speeds.
Our results effectively showcase this advantage in Fig. 10, where we compare
the walking behavior for the maximum attainable speed based on the specified
gait parameters in (a) and (b). It is clear that heel-to-toe walking results
in a noticeably larger footstep. The corresponding velocity profile and phase
portrait are presented in Fig. 10(c) and (d). The proposed method using heel-
to-toe walking achieves an impressive speed of approximately 2.15 m/s, whereas
flat-footed walking only reaches a speed of 1.65 m/s. To provide context, the
average human walking speed is 1.42 m/s. We are able to realize highly dynamic
locomotion behaviors on Cassie using our proposed method.

######  Remark 6.

The previous method [5] was capable of achieving a maximum speed of
approximately 2 m/s on Cassie. However, this was only attainable with much
faster stepping motions, specifically with $T_{\text{SS}}=0.3$s and
$T_{\text{OA}}=0$s. It’s important to note that the resulting gait was
unstable without additional data-driven adaptation, as discussed in [30].

##  VI Conclusion

In conclusion, this paper introduces a novel reduced-order model based
approach to realize multi-domain walking on bipedal robots. Leveraging the S2S
dynamics of the proposed MLIP model, we have demonstrated the ability to
stabilize the robot at arbitrary walking speed, achieving a remarkable maximum
speed of 2.15 m/s. The robustness of the method is demonstrated through push
recovery tests. Importantly, this method can be readily extended to
applications in robotic assistive devices, such as exoskeletons, to enable
human-like multi-domain locomotion for the mobility impaired.

## References

  * [1] A. D. Ames, R. Vasudevan, and R. Bajcsy, “Human-data based cost of bipedal robotic walking,” in _Proceedings of the 14th international conference on Hybrid systems: computation and control_ , 2011, pp. 153–162.
  * [2] S. Kajita, F. Kanehiro, K. Kaneko, K. Yokoi, and H. Hirukawa, “The 3D linear inverted pendulum mode: a simple modeling for a biped walking pattern generation,” in _Proceedings 2001 IEEE/RSJ International Conference on Intelligent Robots and Systems_ , vol. 1. Maui, HI, USA: IEEE, 2001, pp. 239–246.
  * [3] P. V. Ambarish Goswami, Ed., _Humanoid Robotics: A Reference_. Springer Dordrecht, 2018.
  * [4] Y. Gong and J. W. Grizzle, “Zero Dynamics, Pendulum Models, and Angular Momentum in Feedback Control of Bipedal Locomotion,” _Journal of Dynamic Systems, Measurement, and Control_ , vol. 144, no. 12, p. 121006, Dec. 2022.
  * [5] X. Xiong and A. Ames, “3-D Underactuated Bipedal Walking via H-LIP Based Gait Synthesis and Stepping Stabilization,” _IEEE Transactions on Robotics_ , vol. 38, no. 4, pp. 2405–2425, Aug. 2022.
  * [6] M. Kim and S. H. Collins, “Once-Per-Step Control of Ankle Push-Off Work Improves Balance in a Three-Dimensional Simulation of Bipedal Walking,” _IEEE Transactions on Robotics_ , vol. 33, no. 2, pp. 406–418, Apr. 2017.
  * [7] J. P. Reher, A. Hereid, S. Kolathaya, C. M. Hubicki, and A. D. Ames, “Algorithmic foundations of realizing multi-contact locomotion on the humanoid robot durus,” in _Algorithmic Foundations of Robotics XII: Proceedings of the Twelfth Workshop on the Algorithmic Foundations of Robotics_. Springer, 2020, pp. 400–415.
  * [8] R. Sellaouti, O. Stasse, S. Kajita, K. Yokoi, and A. Kheddar, “Faster and Smoother Walking of Humanoid HRP-2 with Passive Toe Joints,” in _2006 IEEE/RSJ International Conference on Intelligent Robots and Systems_ , Oct. 2006, pp. 4909–4914, iSSN: 2153-0866.
  * [9] C. Chevallereau, D. Djoudi, and J. W. Grizzle, “Stable Bipedal Walking With Foot Rotation Through Direct Regulation of the Zero Moment Point,” _IEEE Transactions on Robotics_ , vol. 24, no. 2, pp. 390–401, Apr. 2008.
  * [10] H. Zhao, A. Hereid, W.-l. Ma, and A. D. Ames, “Multi-contact bipedal robotic locomotion,” _Robotica_ , vol. 35, no. 5, pp. 1072–1106, May 2017\.
  * [11] K. Chao and P. Hur, “A step towards generating human-like walking gait via trajectory optimization through contact for a bipedal robot with one-sided springs on toes,” in _2017 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_ , Sept. 2017, pp. 4848–4853, iSSN: 2153-0866.
  * [12] J. Reher, “Dynamic Bipedal Locomotion: From Hybrid Zero Dynamics to Control Lyapunov Functions via Experimentally Realizable Methods,” Ph.D. dissertation, California Institute of Technology, Feb. 2021\.
  * [13] M. H. Raibert, _Legged robots that balance_. MIT press, 1986.
  * [14] F. Sup, A. Bohara, and M. Goldfarb, “Design and Control of a Powered Transfemoral Prosthesis,” _The International Journal of Robotics Research_ , vol. 27, no. 2, pp. 263–273, Feb. 2008, publisher: SAGE Publications Ltd STM.
  * [15] N. P. Fey, A. M. Simon, A. J. Young, and L. J. Hargrove, “Controlling Knee Swing Initiation and Ankle Plantarflexion With an Active Prosthesis on Level and Inclined Surfaces at Variable Walking Speeds,” _IEEE journal of translational engineering in health and medicine_ , vol. 2, p. 2100412, 2014.
  * [16] N. Thatte, T. Shah, and H. Geyer, “Robust and Adaptive Lower Limb Prosthesis Stance Control via Extended Kalman Filter-Based Gait Phase Estimation,” _IEEE Robotics and Automation Letters_ , vol. 4, no. 4, pp. 3129–3136, Oct. 2019.
  * [17] M. F. Eilenberg, H. Geyer, and H. Herr, “Control of a Powered Ankle–Foot Prosthesis Based on a Neuromuscular Model,” _IEEE Transactions on Neural Systems and Rehabilitation Engineering_ , vol. 18, no. 2, pp. 164–173, Apr. 2010.
  * [18] R. Gehlhar, M. Tucker, A. J. Young, and A. D. Ames, “A review of current state-of-the-art control methods for lower-limb powered prostheses,” _Annual Reviews in Control_ , vol. 55, pp. 142–164, 2023.
  * [19] J. W. Grizzle, C. Chevallereau, R. W. Sinnet, and A. D. Ames, “Models, feedback control, and open problems of 3D bipedal robotic walking,” _Automatica_ , vol. 50, no. 8, pp. 1955–1988, Aug. 2014.
  * [20] M. Vukobratović and B. Borovac, “Zero-moment point — thirty five years of its life,” _International Journal of Humanoid Robotics_ , vol. 01, no. 01, pp. 157–173, Mar. 2004.
  * [21] D. E. Orin, A. Goswami, and S.-H. Lee, “Centroidal dynamics of a humanoid robot,” _Autonomous Robots_ , vol. 35, no. 2-3, pp. 161–176, Oct. 2013.
  * [22] I. Danilov, B. Gabbasov, I. Afanasyev, and E. Magid, “ZMP Trajectory from Human Body Locomotion Dynamics Evaluated by Kinect-based Motion Capture System:,” in _International Conference on Computer Vision Theory and Applications_ , 2016, pp. 160–166.
  * [23] Y. Gong and J. Grizzle, “Angular Momentum about the Contact Point for Control of Bipedal Locomotion: Validation in a LIP-based Controller,” Apr. 2021, arXiv:2008.10763 [cs, eess].
  * [24] K. Kouramas, S. Rakovic, E. Kerrigan, J. Allwright, and D. Mayne, “On the Minimal Robust Positively Invariant Set for Linear Difference Inclusions,” in _Proceedings of the 44th IEEE Conference on Decision and Control_ , Dec. 2005, pp. 2296–2301, iSSN: 0191-2216.
  * [25] A. Robotics, https://www.agilityrobotics.com/robots#cassie.
  * [26] K. Bouyarmane, K. Chappellet, J. Vaillant, and A. Kheddar, “Quadratic Programming for Multirobot and Task-Space Force Control,” _IEEE Transactions on Robotics_ , vol. 35, no. 1, pp. 64–77, Feb. 2019.
  * [27] A. Robotics, cassie-mujoco-sim https://github.com/osudrl/cassie-mujoco-sim, 2018.
  * [28] 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_ , Oct. 2012, pp. 5026–5033.
  * [29] I. Birch, W. Vernon, J. Walker, and M. Young, “Terminology and forensic gait analysis,” _Science & Justice_, vol. 55, no. 4, pp. 279–284, July 2015\.
  * [30] M. Dai, X. Xiong, J. Lee, and A. D. Ames, “Data-driven Adaptation for Robust Bipedal Locomotion with Step-to-Step Dynamics,” Aug. 2023, arXiv:2209.08458 [cs, eess].