# Adaptive Force-Based Control of Dynamic Legged Locomotion over Uneven
Terrain

Mohsen Sombolestan and Quan Nguyen M. Sombolestan and Q. Nguyen are with the
Department of Aerospace and Mechanical Engineering, University of Southern
California, Los Angeles, CA 90089, email<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Agile-legged robots have proven to be highly effective in navigating and
performing tasks in complex and challenging environments, including disaster
zones and industrial settings. However, these applications normally require
the capability of carrying heavy loads while maintaining dynamic motion.
Therefore, this paper presents a novel methodology for incorporating adaptive
control into a force-based control system. Recent advancements in the control
of quadruped robots show that force control can effectively realize dynamic
locomotion over rough terrain. By integrating adaptive control into the force-
based controller, our proposed approach can maintain the advantages of the
baseline framework while adapting to significant model uncertainties and
unknown terrain impact models. Experimental validation was successfully
conducted on the Unitree A1 robot. With our approach, the robot can carry
heavy loads (up to 50% of its weight) while performing dynamic gaits such as
fast trotting and bounding across uneven terrains.

###### Index Terms:

Adaptive control, Model predictive control (MPC), Quadruped robots, Unknown
impact model.

##  I Introduction

Legged robots have numerous potential uses, from search and rescue operations
to autonomous construction. To perform these tasks effectively, it is
important for the robot to have an accurate understanding of the environment
it will be operating in. However, due to the complexity of the robot and the
environment, the model of the robot itself might contain a significant level
of uncertainty and affect the robot’s stability, particularly when performing
agile movements. To overcome these challenges, there is a need for the
development of a control framework that can effectively compensate for these
uncertainties in real-time.

The utilization of convex model predictive control (MPC) with the single rigid
body (SRB) model in legged robots [1] has greatly enhanced the real-time
implementation of diverse walking gaits. Unlike the balance controller based
on quadratic programming [2], MPC offers the capability to perform agile
motions like jumping [3, 4] and high-speed bounding [5] for quadruped robots.
Additionally, MPC exhibits robustness in traversing rough and uneven terrains.
However, it is important to note that MPC assumes perfect knowledge of the
dynamic model.

To enhance trajectory tracking in the presence of unknown and changing
disturbances, researchers have explored the combination of MPC with adaptive
control techniques [6, 7, 8]. Additionally, parameter estimation techniques
have been employed to further improve the robustness of the control system
[9]. These approaches aim to adapt the controller and estimate system
parameters to effectively compensate for uncertainties and disturbances,
leading to improved trajectory tracking performance. It is worth noting that
all of these studies were conducted using a position-based controller model.

In this work, we tackle the legged robot locomotion issue in real-world
scenarios with a significant level of uncertainty. The uncertainty can come
from both the robot model and the environment. Since our proposed method is
based on a force controller, it retains the advantage of robustness to uneven
terrain. Thanks to MPC as our baseline controller, our framework can be
extended to different locomotion gaits and trajectories without adjusting the
controller parameters. Additionally, by incorporating the adaptive controller,
our control system can handle significant model uncertainty. As a result, our
approach enables legged robots to move across different terrains with unknown
impact models.

Figure 1: Our proposed adaptive MPC is successfully validated in an experiment
on a Unitree A1 robot while carrying an unknown load of 5 kg (almost 50% of
body weight) on rough terrain. Experimental results video:
https://youtu.be/QmwyysdTk1k.

###  I-A Related Works

####  I-A1 Offline Learning

The offline learner can either leverage a model-based control approach or
learn the control system from scratch. Using a model-based method, researchers
mainly target learning the dynamic to improve the controller performance [10].
One example of this approach is the integration of deep learning with MPC, in
which the proposed model tries to learn the cost or dynamic terms of an MPC
[11]. This hybrid method shows considerable improvement for the aerial robot
[12] when learning the dynamic model from experimental data. The major
limitation of this method is that it is restricted to the dynamic model
learned during the training phase. However, the dynamic model is prone to
frequent changes in real-world scenarios due to environmental uncertainties
and external disturbances.

To overcome the limitations of previous approaches, there has been growing
interest in utilizing reinforcement learning (RL) to train models from
scratch. The key advantage of RL models is their ability to adapt swiftly to
changes in real-world environments due to being trained in diverse
environments with varying properties. In the case of quadruped robots, an RL
model can directly predict appropriate joint torques for traversing different
types of terrain, as demonstrated by Chen et al. [13]. Additionally,
Bellegarda et al. [14] enable quadrupeds to run quickly while carrying unknown
loads by training the model to learn foot positions. However, these methods
heavily rely on domain randomization during training to generalize well to
challenging environments. Yang et al. [15] also propose an end-to-end RL
method that utilizes proprioceptive states and visual feedback to predict
environmental changes.

####  I-A2 Online Learning

To address inaccuracies in model-based controllers, researchers have explored
an alternative approach using online learning, particularly supervised
learning methods [16, 17, 18]. In this approach, the focus is on learning
disturbances online [19], and in some cases, researchers also aim to learn the
dynamics of the system itself [20]. Furthermore, this approach has been
successfully applied for online calibration of kinematic parameters in legged
robots [21]. In addition to that, in a recent study, a Lipschitz network
method has been developed to bridge the model-reality gap in real-time [22].
The online learning method shares a close relationship with adaptive control,
and numerous studies have explored the combination of these two approaches
[23]. This combination aims to leverage the advantages of both methods,
allowing for dynamic adaptation and continuous learning from real-time data to
improve control system performance. Perhaps closest to our work in terms of
online adaption is the learning method presented in [24] for legged robots.
The authors correct the model behind the controller using a supervised learner
while the robot is walking in an unknown environment. The data is collected
during the robot’s operation to learn a linear residual model which can
compensate for system errors. However, in the transition from simulation to
experiment, the acceleration estimators make noisy data required for training
the model. As a result, the method is only applied to estimate the linear
terms since the angular terms data proved to be too noisy to be helpful in the
model.

####  I-A3 Adaptive Control

The goal of adaptive control is to tune the controller’s variables online
during deployment [25]. Adaptive control has been applied for manipulation
tasks to robotic arms [26], mobile robots [27, 28, 29], and quadruped robots
[30, 31]. The conventional Model Reference Adaptive Control (MRAC)
architecture was originally designed for controlling linear systems in the
presence of parametric uncertainties [32, 33]. However, it lacks the ability
to characterize the input/output performance of the system during the
transient phase. To address this limitation and improve the transient
performance of adaptive controllers, the $L_{1}$ adaptive control offers
several advantages over traditional MRAC, such as decoupling adaptation and
robustness within a control framework [34]. In addition, by incorporating a
low-pass filter in adaptation law, the $L_{1}$ adaptive control can provide
stability [35] and transient performance [36]. Therefore, the $L_{1}$ adaptive
control technique guarantees robustness with fast adaptation [37], an
essential criterion in dynamic robotics applications. Recently, by integrating
$L_{1}$ adaptive controller and Bayesian learner, researchers leverage the
fast adaption performance of the $L_{1}$ adaptive controllers and introduce a
safe simultaneous control and learning framework [38, 39].

For legged robots, the adaptive controller has also been employed to find the
value and location of the center of mass [40]. Our work on $L_{1}$ adaptive
control for bipedal robots [41] considers a Control Lyapunov Function
(CLF)-based controller as a closed-loop nonlinear reference model for the
$L_{1}$ adaptive controller. It was validated for the robot’s walking [42] and
running [43]. However, the control framework in this prior work is based on
Hybrid Zero Dynamics [44], which uses joint position control to track the
desired trajectory from optimization for each robot joint. Moreover, in [45],
an adaptive control based on a CLF is designed for quadrupeds to interact with
unknown objects. Then, they combined the criteria derived by adaptive control
as a constraint in an MPC framework. However, adding more inequality
constraints to MPC makes the controller more complex in terms of computation.
In our approach, we compute a residual vector for compensating dynamic
uncertainty, which makes the controller more time-efficient. Additionally, by
employing our method, the robot is able to adapt to terrains with unknown
impact models.

###  I-B Contributions

A preliminary version of this research previously appeared in [46]; however,
this paper presents several novel contributions to the prior work. This work
incorporates the $L_{1}$ adaptive controller into the model predictive control
(MPC). The proposed control system leverages MPC due to its robustness to
uneven terrain, contact constraint, and generalization to different locomotion
gaits. Moreover, by integrating adaptive control into MPC, the proposed model
can compensate for significant model uncertainty. In the previous work [46],
the robot can only perform quasi-static walking; however, in this work, the
robot can perform dynamic motions thanks to MPC. Finally, the authors present
new hardware experiments to demonstrate the effectiveness of the proposed
adaptive MPC (as illustrated in Fig. 1). The main contributions of the paper
are as follows:

  * •

We introduce a novel control system that combines the $L_{1}$ adaptive control
into the force-based control system, designed to address the challenges posed
by model uncertainty in real-world applications.

  * •

Thanks to MPC, our approach offers greater versatility as it can be adapted to
a wide range of locomotion gaits and trajectories. Moreover, our method can
handle terrain uncertainty, allowing the robot to navigate rough terrains,
such as grass and gravel, as well as high-sloped terrain.

  * •

By integrating the adaptive control into MPC, it is possible for quadruped
robots to carry an unknown heavy load (up to 50% of the robot’s weight) across
challenging terrains, with the capability of executing dynamic gaits such as
fast trotting and bounding. This is a significant improvement compared to our
previous work, which only allowed the robot to perform quasi-static walking.

  * •

The combination of using MPC for both the reference model and the real model
in the adaptive controller makes the control system computationally expensive,
leading to potential delays in computation. To ensure real-time performance,
we have developed an update frequency scheme for the control system, which
allows for the optimized allocation of processing resources to each control
component.

  * •

Our proposed approach enables the control system to adapt to terrains with
unknown impact models, such as soft terrain. Traversing soft terrain is a
challenging task for quadruped robots. The A1 robot can walk on double-foam
terrain in different directions using our method. In comparison, the robot
cannot maintain its balance using the baseline controller, resulting in a
collapse.

The remainder of the paper is organized as follows. Sec. II presents the
baseline control architecture for quadruped robots and provides some knowledge
on force-based controllers. In Sec. III, we will briefly present an overview
of our control approach. Then, our proposed adaptive force-based controller
using balance controller and MPC will be elaborated in Sec. IV and Sec. V,
respectively. Furthermore, the numerical and experimental validation are shown
in Sec. VII. Finally, Sec. VIII provides concluding remarks.

##  II PRELIMINARIES

In this section, we present the background on the control architecture of
quadruped robots and describe each control component. According to [47], the
robot’s control system consists of several modules, including a high-level
controller, low-level controller, state estimation, and gait scheduler as
presented in Fig. 2.

A reference trajectory can be generated for high-level control from user input
and state estimation. The gait scheduler defines the gait timing and sequence
to switch between each leg’s swing and stance phases. The high-level part
controls the position of the swing legs and optimal ground reaction force for
stance legs based on the user commands and gait timing. As the baseline for
the stance leg controller, we will use two common approaches: 1) quadratic
program (QP) based balancing controller [2] and 2) model predictive control
(MPC) [1]. The low-level leg control converts the command generated by high-
level control into joint torques for each motor. These modules of the control
architecture will be described briefly in the following subsections. More
details can be found in [47, 1, 2].

###  II-A Gait Scheduler

The A1’s gait is defined by a finite state machine using a leg-independent
phase variable to schedule contact and swing phases for each leg [47]. The
gait scheduler utilizes independent boolean variables to define contact states
scheduled $\bm{s}_{\phi}\in\\{1=contact,0=swing\\}$ and switch each leg
between swing and stance phases. Based on the contact schedule, the controller
will execute either position control during swing or force control during
stance for each leg.

In our previous work [46], we focus on the application of load-carrying tasks,
where the load is unknown to the robot or the control system. Having more legs
on the ground during walking could also mean that the robot could produce a
larger total ground reaction force to support the heavy load. Therefore, we
used a quasi-static walking gait to maximize the number of legs on the ground
during walking (i.e., 3 stance legs and 1 swing leg throughout the gait).
However, in this paper, our framework is not limited by any specific gait.
Similar to the baseline MPC control approach [1], the approach can work for
different gaits by only changing the gait definition in the gait scheduler.

###  II-B Desired Trajectory

The desired trajectory is generated based on the robot’s velocity command. The
robot operator commands $xy$-velocity and yaw rate, then $xy$-position and yaw
are determined by integrating the corresponding velocity. $z$ position
contains a constant value of $0.3m$, and the remaining states (roll, roll
rate, pitch, pitch rate, and $z$-velocity) are always zero.

Figure 2: Baseline Control Structure. Block diagram of a control architecture
for a quadruped robot. For the stance leg control, we use two common baseline
control systems: QP-based balancing controller and MPC.

###  II-C Single Rigid Body (SRB) Model of Robot

Due to the complexity of the legged robot, a simplified rigid-body model has
been used to present the system’s dynamic. This model enables us to calculate
the ground reaction forces (GRFs) in real-time. A few assumptions have been
made to achieve simplified robot dynamics[1]:

Assumption 1: The robot has low inertia legs, so their effect is negligible.

Assumption 2: For small values of roll ($\phi$) and pitch ($\theta$), the
rotation matrix $\bm{R}$ which transforms from the body to world coordinates,
can be approximated as the rotation matrix corresponding to the yaw angle
($\psi$):

$\displaystyle\bm{R}\cong\bm{R}_{z}(\psi)=\left[\begin{array}[]{ccc}\cos(\psi)&-\sin(\psi)&0\\\
\sin(\psi)&\cos(\psi)&0\\\ 0&0&1\end{array}\right]$ (4)

Therefore, by defining the robot’s orientation as a vector of Z-Y-X Euler
angles $\bm{\Theta}=[\phi,\theta,\psi]^{T}$, the rate of change of the robot’s
orientation can be approximated as [1]:

$\displaystyle\dot{\bm{\Theta}}\cong\bm{R}_{z}(\psi)\bm{\omega}_{b}$ (5)

where $\bm{\omega}_{b}$ is the robot’s angular velocity in the world frame.

Assumption 3: For small angular velocity, the following approximation can be
made:

$\displaystyle\frac{d}{dt}(\bm{I}_{G}\mbox{\boldmath$\omega$}_{b})=\bm{I}_{G}\dot{}\mbox{\boldmath$\omega$}_{b}+\mbox{\boldmath$\omega$}_{b}\times(\bm{I}_{G}\mbox{\boldmath$\omega$}_{b})\approx\bm{I}_{G}\dot{}\mbox{\boldmath$\omega$}_{b}$
(6)

where $\bm{I}_{G}\in\mathbb{R}^{3\times 3}$ is the moment of inertia in the
world frame.

Based on the above assumptions, the state representation of the system is as
follows [1]:

$\displaystyle\left[\begin{array}[]{c}\dot{\bm{p}}_{c}\\\ \bm{\dot{\Theta}}\\\
\ddot{\bm{p}}_{c}\\\ \dot{}\mbox{\boldmath$\omega$}_{b}\end{array}\right]=$
$\displaystyle\underbrace{\left[\begin{array}[]{cccc}\bm{0}_{3}&\bm{0}_{3}&\bm{1}_{3}&\bm{0}_{3}\\\
\bm{0}_{3}&\bm{0}_{3}&\bm{0}_{3}&\bm{R}_{z}(\psi)\\\
\bm{0}_{3}&\bm{0}_{3}&\bm{0}_{3}&\bm{0}_{3}\\\
\bm{0}_{3}&\bm{0}_{3}&\bm{0}_{3}&\bm{0}_{3}\end{array}\right]}_{\bm{D}\in\mathbb{R}^{12\times
12}}\underbrace{\left[\begin{array}[]{c}{\bm{p}}_{c}\\\ \bm{\Theta}\\\
\dot{\bm{p}}_{c}\\\
\mbox{\boldmath$\omega$}_{b}\end{array}\right]}_{\bm{X}\in\mathbb{R}^{12}}+$
(19) $\displaystyle\underbrace{\left[\begin{array}[]{c}\bm{0}_{6\times 12}\\\
\bm{M}^{-1}\bm{A}\end{array}\right]}_{\bm{H}\in\mathbb{R}^{12\times
12}}\bm{F}+\left[\begin{array}[]{c}\bm{0}_{6\times 1}\\\
\bm{G}\end{array}\right]$ (24)

with

$\displaystyle\bm{M}=\left[\begin{array}[]{cc}m\bf{1}_{3}&\bm{0}_{3}\\\
\bm{0}_{3}&\bm{I}_{G}\end{array}\right]\in\mathbb{R}^{6\times 6}$ (27)
$\displaystyle\bm{A}=\left[\begin{array}[]{ccc}\bf{1}_{3}&\dots&\bf{1}_{3}\\\
[\bm{p}_{1}-\bm{p}_{c}]\times&\dots&[\bm{p}_{4}-\bm{p}_{c}]\times\end{array}\right]\in\mathbb{R}^{6\times
12}$ (30) $\displaystyle\bm{G}=\left[\begin{array}[]{c}\bm{g}\\\
\bm{0}_{3\times 1}\end{array}\right]\in\mathbb{R}^{6}$ (33)

where $m$ is the robot’s mass, $\bm{g}\in\mathbb{R}^{3}$ is the gravity
vector, $\bm{p}_{c}\in\mathbb{R}^{3}$ is the position of the center of mass
(COM), $\bm{p}_{i}\in\mathbb{R}^{3}$ ($i\in\\{1,2,3,4\\}$) are the positions
of the feet, $\ddot{\bm{p}}_{c}\in\mathbb{R}^{3}$ is body’s linear
acceleration, $\dot{}\mbox{\boldmath$\omega$}_{b}\in\mathbb{R}^{3}$ is angular
acceleration, and
$\bm{F}=[\bm{F}_{1}^{T},\bm{F}_{2}^{T},\bm{F}_{3}^{T},\bm{F}_{4}^{T}]^{T}\in\mathbb{R}^{12}$
are the ground reaction forces acting on each of the robot’s four feet. The
term $[\bm{p}_{i}-\bm{p}_{c}]\times$ is the skew-symmetric matrix representing
the cross product $(\bm{p}_{i}-\bm{p}_{c})\times\bm{F}_{i}$. Note that
$\bm{p}_{i}$ and $\bm{F}_{i}$ are presented in the world frame. Therefore, the
state representation of the system can be rewritten in the compact form:

$\displaystyle\bm{\dot{X}}=\bm{D}\bm{X}+\bm{H}\bm{F}+\left[\begin{array}[]{c}\bm{0}_{6\times
1}\\\ \bm{G}\end{array}\right]$ (36)

###  II-D Balance Controller

One of the baseline control approach for calculating GRFs for quadruped robots
is the balance controller presented in [2] based on quadratic program (QP)
solver. Based on the assumptions presented in Sec. II-C, the approximated
dynamic model between the body acceleration and GRFs is as follows:

$\displaystyle\underbrace{\left[\begin{array}[]{ccc}\bf{1}_{3}&\dots&\bf{1}_{3}\\\
[\bm{p}_{1}-\bm{p}_{c}]\times&\dots&[\bm{p}_{4}-\bm{p}_{c}]\times\end{array}\right]}_{\bm{A}\in\mathbb{R}^{6\times
12}}\bm{F}=\underbrace{\left[\begin{array}[]{c}m(\ddot{\bm{p}}_{c}+\bm{g})\\\
\bm{I}_{G}\dot{}\mbox{\boldmath$\omega$}_{b}\end{array}\right]}_{\bm{b}\in\mathbb{R}^{6}}$
(41)

and the vector $\bm{b}$ in (41) can be rewritten as:

$\displaystyle\bm{b}=\bm{M}(\left[\begin{array}[]{c}\ddot{\bm{p}}_{c}\\\
\dot{}\mbox{\boldmath$\omega$}_{b}\end{array}\right]+\bm{G}).$ (44)

Since the model (41) is linear, the controller can naturally be formulated as
the following QP problem [48], which can be solved in real-time at $1~{}kHz$:

$\displaystyle\bm{F}^{*}=\underset{\bm{F}\in\mathbb{R}^{12}}{\operatorname{argmin}}\>\>$
$\displaystyle(\bm{A}\bm{F}-\bm{b}_{d})^{T}\bm{S}(\bm{A}\bm{F}-\bm{b}_{d})$
$\displaystyle+\gamma_{1}\|\bm{F}\|^{2}+\gamma_{2}\|\bm{F}-\bm{F}_{\textrm{prev}}^{*}\|^{2}$
(45) s.t.
$\displaystyle\quad\>\>\underline{\bm{d}}\leq\bm{C}\bm{F}\leq\bar{\bm{d}}$
$\displaystyle\quad\>\>\bm{F}_{swing}^{z}=0$

where $\bm{b}_{d}$ is the desired dynamics. The idea of designing $\bm{b}_{d}$
will be elaborated in Sec. IV-A. The cost function in (II-D) includes terms
that consider three goals, including (1) driving the COM position and
orientation to the desired trajectories; (2) minimizing the force commands;
and (3) minimizing the change of the current solution $\bm{F}^{*}$ with
respect to the solution from the previous time-step, $\bm{F}^{*}_{prev}$. The
priority of each goal in the cost function is defined by the weight parameters
$\bm{S}\in\mathbb{R}^{6\times 6}$, $\gamma_{1}$, $\gamma_{2}$ respectively.

The constraints in the QP formulation enforce friction constraints, input
saturation, and contact constraints. The constraint
$\underline{\bm{d}}\leq\bm{C}\bm{F}\leq\bar{\bm{d}}$ ensures that the
optimized forces lie inside the friction pyramid and the normal forces stay
within a feasible range. More details can be found in [2]. Besides the
friction constraint, we will enforce the force constraints for the swing legs,
$\bm{F}_{swing}=\bm{0}$. The swing legs are then kept in the posing position
until it switches to the stance phase. More details on swing leg control are
provided in Sec. II-F.

###  II-E SRB-based Convex MPC

The calculation of GRFs in quadruped robots is often approached through Model
Predictive Control (MPC) [1]. This method determines the optimal sequence of
inputs over a finite-time horizon, taking into account any constraints within
the dynamic model. Every time MPC is executed in the control system, only the
first computed control input from the MPC cycle is applied. The inputs
determined over the finite time horizon are only used for the optimization
problem and are not directly applied in the control system.

To have the dynamic equation in the convenient state-space form, gravity
should be added to the state. So, the system can represent as:

$\displaystyle\bm{\dot{X}}^{c}=\bm{D}^{c}\bm{X}^{c}+\bm{H}^{c}\bm{F}$ (46)

where

$\displaystyle\bm{X}^{c}=\left[\begin{array}[]{c}\bm{p}_{c}\\\ \bm{\Theta}\\\
\dot{\bm{p}}_{c}\\\ \mbox{\boldmath$\omega$}_{b}\\\
||\bm{g}||\end{array}\right]\in\mathbb{R}^{13}$ (52)
$\displaystyle\bm{D}^{c}=\left[\begin{array}[]{ccccc}\bm{0}_{3}&\bm{0}_{3}&\bm{1}_{3}&\bm{0}_{3}&\bm{0}_{3\times
1}\\\ \bm{0}_{3}&\bm{0}_{3}&\bm{0}_{3}&\bm{R}_{z}(\psi)&\bm{0}_{3\times 1}\\\
\bm{0}_{3}&\bm{0}_{3}&\bm{0}_{3}&\bm{0}_{3}&\frac{\bm{g}}{||\bm{g}||}\\\
\bm{0}_{3}&\bm{0}_{3}&\bm{0}_{3}&\bm{0}_{3}&\bm{0}_{3\times 1}\\\
\bm{0}_{1\times 3}&\bm{0}_{1\times 3}&\bm{0}_{1\times 3}&\bm{0}_{1\times
3}&0\end{array}\right]\in\mathbb{R}^{13\times 13}$ (58)
$\displaystyle\bm{H}^{c}=\left[\begin{array}[]{c}\bm{0}_{6\times 12}\\\
\bm{M}^{-1}\bm{A}\\\ \bm{0}_{1\times
12}\end{array}\right]\in\mathbb{R}^{13\times 12}$ (62)

We consider a linear MPC problem with horizon length $k$ as follows:

$\displaystyle\min_{\bm{F}_{i}}\quad$
$\displaystyle\sum_{i=0}^{k-1}{\bm{e}_{i+1}}^{T}\bm{Q}_{i}{\bm{e}_{i+1}}+{\bm{F}_{i}}^{T}\bm{R}_{i}\bm{F}_{i}$
(63) s.t.
$\displaystyle\bm{X}^{c}_{i+1}=\bm{D}_{t,i}\bm{X}^{c}_{i}+\bm{H}_{t,i}\bm{F}_{i}$
$\displaystyle\underline{\bm{d}}\leq\bm{C}\bm{F}_{i}\leq\bar{\bm{d}}$

where $\bm{F}_{i}$ is the computed ground reaction forces at time step $i$,
$\bm{Q}_{i}$ and $\bm{R}_{i}$ are diagonal positive semi-definite matrices,
$\bm{D}_{t,i}$ and $\bm{H}_{t,i}$ are discrete time system dynamics matrices.
The ${\bm{e}_{i+1}}$ is the system state error at time step $i$ define as
$\bm{e}=[\bm{e}_{p},~{}\dot{\bm{e}}_{p}]^{T}\in\mathbb{R}^{12}$, with

$\displaystyle\bm{e}_{p}=\left[\begin{array}[]{c}\bm{p}_{c}-\bm{p}_{c,d}\\\
\log(\bm{R}_{d}\bm{R}^{T})\end{array}\right]\in\mathbb{R}^{6},\quad\dot{\bm{e}}_{p}=\left[\begin{array}[]{c}\dot{\bm{p}}_{c}-\dot{\bm{p}}_{c,d}\\\
\mbox{\boldmath$\omega$}_{b}-\mbox{\boldmath$\omega$}_{b,d}\end{array}\right]\in\mathbb{R}^{6},$
(68)

where $\bm{p}_{c,d}\in\mathbb{R}^{3}$ is the desired position of COM,
$\dot{\bm{p}}_{c,d}\in\mathbb{R}^{3}$ is the desired body’s linear velocity,
and $\mbox{\boldmath$\omega$}_{b,d}\in\mathbb{R}^{3}$ is the desired body’s
angular velocity. The desired and actual body orientations are described using
rotation matrices $\bm{R}_{d}\in\mathbb{R}^{3\times 3}$ and
$\bm{R}\in\mathbb{R}^{3\times 3}$, respectively. The orientation error is
obtained using the exponential map representation of rotations [49, 50], where
the $log(.):\mathbb{R}^{3\times 3}\to\mathbb{R}^{3}$ is a mapping from a
rotation matrix to the associated rotation vector [2]. The constraint
$\underline{\bm{d}}\leq\bm{C}\bm{F}_{i}\leq\bar{\bm{d}}$ is equivalent to the
constraint in equation (II-D) at time step $i$.

###  II-F Swing Leg Control

For the swing legs, the final footstep location for each leg is calculated
from the corresponding hip location using a linear combination of Raibert
heuristic [51], and a feedback term from the capture point formulation [52,
47]. The final footstep locations ($\bm{p}_{f,i}$) are projected on an assumed
ground plane and are calculated by:

$\displaystyle\bm{p}_{f,i}=\bm{p}_{h,i}+\frac{{T}_{\bm{c}_{\phi}}}{2}\bm{\dot{p}}_{c,d}+\sqrt{\frac{\bm{z}_{0}}{\|\bm{g}\|}}(\bm{\dot{p}}_{c}-\bm{\dot{p}}_{c,d})$
(69)

where ${T}_{\bm{c}_{\phi}}$ is the stance time scheduled, $\bm{z}_{0}$ is the
height of locomotion and $\bm{p}_{h,i}\in\mathbb{R}^{3}$ is the position of
the corresponding hip $i$. A Beizer curve calculates the desired swing
trajectory (including desired position $\bm{p}_{d,i}$ and velocity
$\bm{v}_{d,i}$) for swing legs which starts from the initial lift-off position
$\bm{p}_{0,i}$ and ends at the final touch-down location $\bm{p}_{f,i}$.

###  II-G Low-level Control

The low-level leg control can generate joint torque commands from the high-
level controller. For low-level force control, the controller transforms the
force vector to the hip frame by rotation matrix $\bm{R}$. Then, joint torques
are calculated as follows:

$\displaystyle\bm{\tau}_{stance,i}=-{\bm{J}(\bm{q}_{i})}^{T}\bm{R}^{T}\bm{F}_{i}$
(70)

where $\bm{J}(\bm{q}_{i})\in\mathbb{R}^{3\times 3}$ is the leg Jacobian matrix
and $\bm{q_{i}}$ is the joints angle of leg $i$-th.

To track the desired swing trajectory for each foot, a PD controller with a
feedforward term is used to compute joint torques [47]:

$\displaystyle\bm{\tau}_{swing,i}=\bm{J}(\bm{q}_{i})^{T}[\bm{K}_{p,p}(\bm{p}_{d,i}-\bm{p}_{i})+\bm{K}_{d,p}(\bm{v}_{d,i}-\bm{v}_{i})]$
(71)

where $\bm{p}_{d,i}$ and $\bm{v}_{d,i}$ are desired foot position and
velocity, respectively, $\bm{p}_{i}$ and $\bm{v}_{i}$ are actual foot position
and velocity in the robot’s frame, $\bm{K}_{p,p}\in\mathbb{R}^{3\times 3}$ and
$\bm{K}_{d,p}\in\mathbb{R}^{3\times 3}$ are the diagonal matrices of the
proportional and derivative gains.

##  III Overview of the Proposed Approach

This section will present an overview of our novel control architecture to
incorporate adaptive control into the force control framework. While our
approach is not limited to any specific adaptive control approach, we decide
to use $L_{1}$ adaptive control [37, 41] thanks to its advancement in
guaranteeing fast adaptation and smooth control signals. Note that our
proposed control system is designed for the stance leg control part in the
control architecture of the quadruped robot (see Fig. 2).

Our prior work [41] introduced an adaptive control based on Hybrid Zero
Dynamics (HZD) [53] for bipedal robots. HZD is a common control approach for
bipedal robots since it can handle hybrid and underactuated dynamics
associated with this kind of robot. In this paper, however, our approach
leverages the combination of the adaptive control and force control system,
which calculates ground reaction forces (GRFs) to achieve highly dynamic
locomotion for quadrupeds [47]. The use of force control in legged robot
systems has several key benefits, including increased robustness in the
presence of challenging terrains [2] and the ability to accommodate a wide
range of dynamic movements [1], such as various types of locomotion gaits. By
combining force control with adaptive control strategies that compensate for
model uncertainty, achieving an enhanced control system with these advantages
is possible.

The overview of our proposed adaptive force-based control system is presented
in Fig. 3a. By incorporating a $L_{1}$ adaptive controller, we aim to design a
combined controller. The force-based controller calculates the optimal GRFs
for following the desired trajectory. The adaptive controller calculates the
residual parameters for compensating the nonlinear model uncertainty
$\bm{\theta}$ in the system dynamic. Therefore, the goal is to adjust adaptive
control signal $\bm{u}_{a}$ as well as adaptation law to estimate the model
uncertainty ($\hat{\bm{\theta}}$) correctly and make the real model follows
the reference model. For the reference model, we employ a similar linear model
described in (36), and we will update the reference model in real-time using
an ODE solver. Moreover, the vector of uncertainties estimation
$\hat{\bm{\theta}}$ typically has high frequency due to fast estimation in the
adaptation law. Thus, we employ a low-pass filter to obtain smooth control
signals. We use the same swing leg control to appropriately synchronize the
reference and real models. This means that we also use the real model’s foot
position for the reference model.

In the following sections, we will elaborate on integrating two different
force-based control as the baseline controller into the adaptive control.
First, in Sec. IV, we will describe the proposed method using a QP-based
balancing controller, as presented in Fig. 3b. Then, in Sec. V, we will show
how to incorporate MPC into the adaptive controller in detail, as illustrated
in Fig. 3c.

(a) Main structure

(b) Adaptive balance controller

(c) Adaptive MPC

Figure 3: Proposed adaptive force-based control system diagram. a) The main
structure of the proposed adaptive force-based control system, b) Block
diagram of the proposed adaptive QP-based balancing controller, c) Block
diagram of the proposed adaptive MPC. Each dashed line indicates the update
frequency for control components

##  IV Adaptive force-based control Using the Balance Controller

In this section, we use the balance controller as the force-based controller,
previously demonstrated in [46]. In Sec. V, we will present our control
framework for integrating the $L_{1}$ adaptive control into MPC.

###  IV-A Closed-loop Dynamics

The $L_{1}$ adaptive control is basically designed for trajectory tracking;
however, the goal of the balance controller is to compute optimal GRFs. Hence,
to integrate the balance controller presented in Sec. II-D into $L_{1}$
adaptive control, we should relate the linear model described in (41) into the
closed-loop dynamics.

Let us consider the system state error ($\bm{e}$) according to equation (68)
as the state variable. Therefore, the closed-loop error dynamics in state-
space form can be represented as follow:

$\displaystyle\dot{\bm{e}}=\bm{D}_{l}\bm{e}+\bm{B}\bm{u},$ (72)

where

$\displaystyle\bm{D}_{l}=\left[\begin{array}[]{cc}\bm{0}_{6}&\bm{1}_{6}\\\
\bm{0}_{6}&\bm{0}_{6}\end{array}\right]\in\mathbb{R}^{12\times
12},\quad\bm{B}=\left[\begin{array}[]{c}\bm{0}_{6}\\\
\bm{1}_{6}\end{array}\right]\in\mathbb{R}^{12\times 6}$ (77)

and $\bm{u}\in\mathbb{R}^{6}$ is the control input function. By employing a PD
control law, we have

$\bm{u}=\begin{bmatrix}-\bm{K}_{P}&-\bm{K}_{D}\end{bmatrix}\bm{e},$ (78)

where $\bm{K}_{P}\in\mathbb{R}^{6\times 6}$ and
$\bm{K}_{D}\in\mathbb{R}^{6\times 6}$ are diagonal positive definite matrices.
According to definition of matrices $\bm{D}_{l}$ and $\bm{B}$, from equation
(72) it can be obtained that:

$\displaystyle\ddot{\bm{e}}_{p}=\left[\begin{array}[]{c}\ddot{\bm{p}}_{c}-\ddot{\bm{p}}_{c,d}\\\
\dot{}\mbox{\boldmath$\omega$}_{b}-\dot{}\mbox{\boldmath$\omega$}_{b,d}\end{array}\right]=\bm{u},$
(81)

where $\ddot{\bm{e}}_{p}$ is the derivative of $\dot{\bm{e}}_{p}$ presented in
(68), $\ddot{\bm{p}}_{c,d}$ and $\dot{}\mbox{\boldmath$\omega$}_{b,d}$ are the
desired COM linear acceleration and the desired angular acceleration,
respectively. Since the desired trajectory is obtained from the velocity
command, both desired accelerations $\ddot{\bm{p}}_{c,d}$ and
$\dot{}\mbox{\boldmath$\omega$}_{b,d}$ are zero vectors. Then from (44) and
(81), the desired dynamics can be given by:

$\displaystyle\bm{b}_{d}=\bm{M}(\bm{u}+\bm{G}),$ (82)

where $\bm{M}$ and $\bm{G}$ are defined in (27). By substituting (82) into the
QP problem (II-D), we can obtain the optimal GRFs as the input for the low-
level leg controller. The objective of the QP formulation in equation (II-D)
is to find a solution that ensures the actual dynamics $\bm{A}\bm{F}$ match
the desired dynamics $\bm{b}_{d}$. In general, the QP-based balance controller
is capable of achieving the desired control input function outlined in
equation (78), thus keeping the error $\bm{e}$ within a certain range.
However, if the desired dynamics vector $\bm{b}_{d}$ violates any of the
inequality constraints, such as force limits or friction constraints, the
controller may yield an optimal solution $\bm{F}^{*}$ that may not completely
align with the desired dynamics. With this solution, the optimal dynamic
${\bm{b}_{d}}^{*}$ and $\bm{u}^{*}$ can be written as:

$\displaystyle{\bm{b}_{d}}^{*}=\bm{A}\bm{F}^{*},$ (83)
$\displaystyle\bm{u}^{*}=\bm{M}^{-1}~{}{\bm{b}_{d}}^{*}-\bm{G}$ (84)

where in the appendix, we will show that the $\bm{u}^{*}$ remains within a
bounded range.

Note that the optimal ground reaction force $\bm{F}^{*}$ serves as the control
input for the robot and the variable $\bm{u}^{*}$ acts as an input for the
closed-loop dynamic. The closed-loop structure for the robot is depicted in
Fig. 3b (the green dashed line).

###  IV-B Effects of Uncertainty on Dynamic

If we consider uncertainty in the dynamic equation (36) and assume that the
matrices $\bm{D}$ and $\bm{H}$ are not accurate, then we need to present the
dynamic based on the nominal matrices $\bar{\bm{D}}$, $\bar{\bm{H}}$. The
model uncertainty mostly comes from inaccurate values for mass, inertia, and
foot position with respect to the center of mass. In addition to that, various
terrain (e.g., rough terrain or soft terrain) might have a different impact on
the robot, and it is unknown in a practical situation. Therefore, terrain
uncertainty should also be considered in the dynamic model. In this section,
we solely derive our control equations based on the model uncertainty. In Sec.
VI, we will elaborate on how our proposed control system can also consider
terrain uncertainty.

There is another parameter involved in the dynamic equation, namely the yaw
angle. This angle is obtained through the state estimation, and we assumed
that the state estimation has minimal uncertainty. According to the definition
of matrices $\bm{D}$ and $\bm{H}$ in (19), the inaccurate value of the dynamic
parameter mentioned above reflects on the $\bm{H}$ matrix. Therefore, the
dynamic equation in the presence of uncertainty can be represented as:

$\displaystyle\bm{\dot{X}}=\bm{D}\bm{X}+(\bar{\bm{H}}+\tilde{\bm{H}})\bm{F}+\left[\begin{array}[]{c}\bm{0}_{6\times
1}\\\ \bm{G}\end{array}\right]$ (87)

where $\tilde{\bm{H}}$ represent the uncertainty in matrix $\bm{H}$. It is
worth noting that according to the definition of $\bm{H}$ in equation (52),
the first six rows of $\bm{H}$ consist of zeros. Thus, we can rephrase the
dynamic equation (87) as follows:

$\displaystyle\bm{\dot{X}}=\bm{D}\bm{X}+\bar{\bm{H}}\bm{F}+\bm{B}\bm{G}+\bm{B}\bm{\theta}$
(88)

where $\bm{\theta}\in\mathbb{R}^{6}$ is the vector of uncertainty for six
corresponding equations and is defined as follows:

$\displaystyle\bm{\theta}\mathrel{\ensurestackMath{\stackon[1pt]{=}{\scriptstyle\Delta}}}\bm{B}^{T}\tilde{\bm{H}}$

With reference to the state representation given by equation (88), the vector
$\bm{\theta}$ can be interpreted as a time-varying disturbance affecting the
body and orientation accelerations.

The uncertainty vector $\bm{\theta}$ depends on both time $t$ and $\bm{F}$.
Since $\bm{F}$ is obtained through the QP problem (II-D), it is a function of
$\bm{b}_{d}$. Furthermore, $\bm{b}_{d}$ is a function of $\bm{u}$ according to
(82). Considering that $\bm{u}$ is determined by the PD control (78), we can
conclude that $\bm{\theta}$ is a function of both the tracking error $\bm{e}$
and time $t$. As a result, for any given time $t$, it is always possible to
find $\bm{\alpha}(t)\in\mathbb{R}^{6}$ and $\bm{\beta}(t)\in\mathbb{R}^{6}$
satisfying [34]:

$\displaystyle\bm{\theta}(\bm{e},t)=\bm{\alpha}(t)||\bm{e}||+\bm{\beta}(t).$
(89)

###  IV-C Designing Adaptive Controller for Compensating the Uncertainty

By incorporating $L_{1}$ adaptive controller, we want to design a combined
controller $\bm{u}=\bm{u}_{1}+\bm{u}_{2}$, where $\bm{u}_{1}$ is the control
input to follow the desired trajectory for the nominal model as presented in
(78) and $\bm{u}_{2}$ is to compensate the nonlinear model uncertainties
$\bm{\theta}$. Therefore, the goal is to adjust the control signal
$\bm{u}_{2}$ so that the real model can follow the reference model. For the
reference model, we employ a similar linear model described in (41) which,
instead of $\bm{M}$, the nominal matrix $\bar{\bm{M}}$ is being used. The
diagram of our proposed force-based adaptive control based on a balance
controller is presented in Fig. 3b.

The duplicate version of equation (88) for state space representation
presented in (72) by considering combined controller
$\bm{u}=\bm{u}_{1}+\bm{u}_{2}$ is as follows:

$\dot{\bm{e}}=\bm{D}_{l}\bm{e}+\bm{B}\bm{u}_{1}+\bm{B}(\bm{u}_{2}+\bm{\theta}).$
(90)

Note that the vector of uncertainty $\bm{\theta}$ in equations (88) and (90)
are not the same since the state vector of equation (88) is $\bm{X}$ while the
state vector of equation (90) is system error $\bm{e}$.

The state representation for the reference model can be expressed as follows:

$\displaystyle\dot{\hat{\bm{e}}}=\bm{D}_{l}\hat{\bm{e}}+\bm{B}\hat{\bm{u}}_{1}+\bm{B}(\bm{u}_{2}+\hat{\bm{\theta}}),$
(91)

where,

$\displaystyle\hat{\bm{\theta}}=\hat{\bm{\alpha}}||\bm{e}||+\hat{\bm{\beta}},$
(92)

and $\hat{\bm{u}}_{1}$ is defined as:

$\hat{\bm{u}}_{1}=\begin{bmatrix}-\bm{K}_{P}&-\bm{K}_{D}\end{bmatrix}\hat{\bm{e}}.$
(93)

To compensate the estimated uncertainty $\hat{\bm{\theta}}$, we can just
simply choose $\bm{u}_{2}=-\hat{\bm{\theta}}$ to obtain

$\dot{\hat{\bm{e}}}=\bm{D}_{l}\hat{\bm{e}}+\bm{B}\hat{\bm{u}}_{1}.$ (94)

However, $\hat{\bm{\theta}}$ typically has high frequency due to fast
estimation in the adaptation law. Therefore, we employ a low-pass filter to
obtain smooth control signals as follows:

$\displaystyle\bm{u}_{2}=-C(s)\hat{\bm{\theta}},$ (95)

where $C(s)$ is a second-order low-pass filter with a magnitude of 1:

$\displaystyle
C(s)=\frac{{\omega_{n}}^{2}}{s^{2}+2\zeta\omega_{n}s+{\omega_{n}}^{2}}.$ (96)

According to (82), the $\bm{b}_{d}$ for the real model in the presence of
uncertainty get the following form:

$\displaystyle\bm{b}_{d}=\bar{\bm{M}}(\bm{u}_{1}+\bm{u}_{2}+\bm{G}).$ (97)

Respectively, $\hat{\bm{b}}_{d}$ for reference model is as follows:

$\displaystyle\hat{\bm{b}}_{d}=\bar{\bm{M}}(\hat{\bm{u}}_{1}+\bm{u}_{2}+\hat{\bm{\theta}}+\bm{G}).$
(98)

The QP solver outlined in equation (II-D) allows us to obtain the optimal GRFs
for the real model. Similarly, the optimal GRFs $\hat{\bm{F}}$ for the
reference model can be obtained as follows:

$\displaystyle\hat{\bm{F}}^{*}=\underset{\hat{\bm{F}}\in\mathbb{R}^{12}}{\operatorname{argmin}}\>\>$
$\displaystyle(\hat{\bm{A}}\hat{\bm{F}}-\hat{\bm{b}}_{d})^{T}\bm{S}(\hat{\bm{A}}\hat{\bm{F}}-\hat{\bm{b}}_{d})$
$\displaystyle+\gamma_{1}\|\hat{\bm{F}}\|^{2}+\gamma_{2}\|\hat{\bm{F}}-\hat{\bm{F}}_{\textrm{prev}}^{*}\|^{2}$
(99) s.t.
$\displaystyle\quad\>\>\underline{\bm{d}}\leq\bm{C}\hat{\bm{F}}\leq\bar{\bm{d}}$
$\displaystyle\quad\>\>\hat{\bm{F}}_{swing}^{z}=0.$

Define the difference between the real model and the reference model
$\tilde{\bm{e}}=\hat{\bm{e}}-\bm{e}$, we then have,

$\displaystyle\dot{\tilde{\bm{e}}}=\bm{D}_{l}\tilde{\bm{e}}+\bm{B}\tilde{\bm{u}}_{1}+\bm{B}(\tilde{\bm{\alpha}}||\bm{e}||+\tilde{\bm{\beta}}),$
(100)

where

$\displaystyle\tilde{\bm{u}}_{1}=\hat{\bm{u}}_{1}-\bm{u}_{1},~{}\tilde{\bm{\alpha}}=\hat{\bm{\alpha}}-\bm{\alpha},~{}\tilde{\bm{\beta}}=\hat{\bm{\beta}}-\bm{\beta}.$
(101)

As a result, we will estimate $\bm{\theta}$ indirectly through $\bm{\alpha}$
and $\bm{\beta}$, or the values of $\hat{\bm{\alpha}}$ and $\hat{\bm{\beta}}$
computed by the following adaptation laws based on the projection operators
[54],

$\displaystyle\dot{\hat{\bm{\alpha}}}=\bm{\Gamma}\text{Proj}(\hat{\bm{\alpha}},\bm{y}_{\alpha}),~{}\dot{\hat{\bm{\beta}}}=\bm{\Gamma}\text{Proj}(\hat{\bm{\beta}},\bm{y}_{\beta})$
(102)

where $\bm{\Gamma}\in\mathbb{R}^{6\times 6}$ is a symmetric positive definite
matrix. The projection functions $\bm{y}_{\alpha}\in\mathbb{R}^{6}$ and
$\bm{y}_{\beta}\in\mathbb{R}^{6}$ are:

$\displaystyle\bm{y}_{\alpha}$
$\displaystyle=-{\bm{B}}^{T}\bm{P}\tilde{\bm{e}}||\bm{e}||,$
$\displaystyle\bm{y}_{\beta}$
$\displaystyle=-{\bm{B}}^{T}\bm{P}\tilde{\bm{e}},$ (103)

where $\bm{P}\in\mathbb{R}^{12\times 12}$ is a positive definite matrix that
is defined according to the stability criteria using the Lyapunov equation.
Moreover, the stability proof of the system is provided in the appendix.

##  V Adaptive Force-based Control using MPC

Model predictive control (MPC) has been widely used across various fields,
from finance to robotics. One of MPC’s main advantages is its ability to
handle complex systems with multiple inputs and outputs while considering hard
control constraints [55]. MPC has also been applied to quadruped robots,
providing stable locomotion [1]. Thanks to dynamic prediction in MPC, by using
the same control framework, it can achieve different dynamic locomotion gaits.
However, MPC’s limitations become evident when dealing with significant
uncertainty in the dynamic model. For instance, in the case of a quadruped
robot carrying an unknown heavy load, MPC fails to track the desired state
trajectory, resulting in unstable behavior and deviation from the desired
trajectory, especially with dynamic gaits like bounding. Furthermore, the
ability of a robot to traverse soft terrain, where the impact model is
unknown, can present a significant challenge. Our proposed approach can tackle
this challenge effectively, and we will discuss the details of how it handles
the terrain unknown impact model in Sec. VI.

In the previous section Sec. IV, we presented an adaptive force-based control
framework based on the balance controller. The balance controller relies on a
quadratic program (QP) solver, which is simple to put into practice and well-
suited for motions that are slow and safe, like standing and quasi-static
walking. Additionally, the balance controller is an instantaneous control
technique, meaning it does not predict the robot’s future movement. As a
result, the balance controller proves to be ineffective in fast-paced, highly
dynamic scenarios. On the other hand, MPC has shown great potential in
handling agile motions, even when it comes to underactuated gaits such as
bounding.

In this section, we will present a novel control architecture to integrate
adaptive control into the MPC framework. By this proposed framework, we can
achieve fast and robust locomotion in the presence of uncertainties. This
framework can also be extended to accommodate various dynamic gaits, such as
trotting and bounding, in legged robots. As we discussed in a previous
section, our approach is not restricted to a specific type of adaptive
control, but we have chosen to utilize $L_{1}$ adaptive control, which has
demonstrated advantages over other adaptive control techniques. The first step
in integrating $L_{1}$ adaptive control and MPC is to understand the
importance of a reference model and the challenges in synchronizing the real
model and reference model. We then present our proposed adaptive MPC, which
combines conventional MPC [1] with adaptive control. Finally, we address the
challenge of real-time computation while having two MPCs in our control
system. We will elaborate on how to adjust the frequency of each control
component in an optimized manner to allocate enough computation resources for
critical control parts and achieve real-time computation.

###  V-A Reference Model

Our method aims to design a combined controller based on MPC and $L_{1}$
adaptive control that the real model follows the reference model. In
accordance with our previous discussion in Sec. IV-C, the combined controller
incorporates a control signal $\bm{u}_{2}$ to account for model uncertainty,
as indicated in equation (90). In this section, the auxiliary control signal
for this purpose is $\bm{u}_{a}\in\mathbb{R}^{6}$, thus, the uncertain dynamic
equation (88) can be rewritten as follow:

$\displaystyle\bm{\dot{X}}=\bm{D}\bm{X}+\bar{\bm{H}}\bm{F}+\bm{B}\bm{G}+\bm{B}(\bm{u}_{a}+\bm{\theta}).$
(104)

The reference model is similar to the quasi-linear model described in (36)
which, instead of $\bm{H}$, the nominal matrix $\bar{\bm{H}}$ is being used.
The proposed adaptive MPC diagram is presented in Fig. 3c.

We consider a reference model for $L_{1}$ adaptive control that arises from
MPC. The MPC method is computationally expensive, but replacing it with other
simpler control methods, such as the balance controller while simulating the
robot’s performance using dynamic gaits such as bounding is impossible. The
reason is that in bounding gait, the robot’s two feet on either the front or
rear side touch the ground at each time step, making it challenging to
accurately control the height and pitch angle. The MPC approach balances the
error in the height and pitch angle and, based on the predicted dynamics of
the system in the future, computes the optimal ground reaction forces. As seen
in Fig. 4, the center of mass (COM) height oscillates around the desired
value. Thus, the underactuated nature of certain gaits like bounding
necessitates the use of MPC as the control system for the reference model.

When implementing MPC for a reference model, one challenge is ensuring that
the reference model is synchronized with the real model. This is particularly
important when the robot performs a gait with a periodic behavior, such as
bounding (see Fig. 4). In order to correctly compare the real model with the
reference model, both should have the same gait schedule. Additionally, the
adaptive MPC proposed for legs in the stance phase is independent of the swing
leg control. However, the foot position is crucial in calculating the moment
of ground reaction force around the center of mass. Therefore, to maintain
consistency between the real and reference models, it is important to ensure
that the real robot’s foot position is fed into the reference model as shown
in Fig. 3c.

Figure 4: Motion snapshot of the robot with bounding gaits. The quadruped’s
center of mass motion (yellow line) cannot be easily predicted with a simple
controller. This illustration can represent the importance of using MPC for
reference model

The reference model can be expressed as follows:

$\displaystyle\dot{\hat{\bm{X}}}=\bm{D}\hat{\bm{X}}+\bm{\bar{H}}\hat{\bm{F}}+\bm{B}\bm{G}+\bm{B}(\bm{u}_{a}+\hat{\bm{\theta}}),$
(105)

where

$\displaystyle\hat{\bm{\theta}}=\hat{\bm{\alpha}}||\bm{e}||+\hat{\bm{\beta}}.$
(106)

In this case, similar to Sec. IV, we use a second-order low-pass filter, same
as (96). Therefore, the auxiliary control signal would be:

$\displaystyle\bm{u}_{a}=-C(s)\hat{\bm{\theta}}.$ (107)

By defining the difference between the real model and the reference model
$\tilde{\bm{X}}=\hat{\bm{X}}-\bm{X}$, we then have:

$\displaystyle\dot{\tilde{\bm{X}}}=\bm{D}\tilde{\bm{X}}+\bar{\bm{H}}\tilde{\bm{F}}+\bm{B}(\tilde{\bm{\alpha}}||\bm{e}||+\tilde{\bm{\beta}}),$
(108)

where

$\displaystyle\tilde{\bm{F}}=\hat{\bm{F}}-\bm{F},~{}\tilde{\bm{\alpha}}=\hat{\bm{\alpha}}-\bm{\alpha},~{}\tilde{\bm{\beta}}=\hat{\bm{\beta}}-\bm{\beta}.$
(109)

Since the desired trajectory for both the real model and the reference model
is the same ($\bm{X}_{d}=\hat{\bm{X}}_{d}$), the difference between the real
model and reference model can be defined as:

$\displaystyle\tilde{\bm{X}}=(\hat{\bm{X}}-\hat{\bm{X}}_{d})-(\bm{X}-\bm{X}_{d})=\hat{\bm{e}}-\bm{e}=\tilde{\bm{e}}.$
(110)

Therefore, equation (108) is equal to the following equation:

$\displaystyle\dot{\tilde{\bm{e}}}=\bm{D}\tilde{\bm{e}}+\bar{\bm{H}}\tilde{\bm{F}}+\bm{B}(\tilde{\bm{\alpha}}||\bm{e}||+\tilde{\bm{\beta}}).$
(111)

The adaption laws and projection functions for computing the value of
$\bm{\alpha}$ and $\bm{\beta}$ are the same as equations (102) and (IV-C),
respectively. Moreover, the stability of the control system can be proven
using the same logic provided in the appendix.

###  V-B Adaptive MPC

After computing the auxiliary control signal $\bm{u}_{a}$ using the adaptive
controller presented in the previous subsection, we will integrate the
$\bm{u}_{a}$ with the conventional MPC for legged locomotion [1] and propose
our adaptive MPC framework. We treat the auxiliary control signal $\bm{u}_{a}$
as a residual vector in the system’s equation to compensate for dynamic
uncertainty. Therefore, the $\bm{u}_{a}$ should be combined into the state
vector and the equation (104) can be written as follow:

$\displaystyle\bm{\dot{\eta}}=\bm{D}^{e}\bm{\eta}+\bar{\bm{H}}^{e}\bm{F}+\bm{B}^{e}\bm{\theta}$
(112)

with the following extended matrices:

$\displaystyle\bm{\eta}=\left[\begin{array}[]{c}\bm{X}^{c}\\\
\hline\cr\bm{u}_{a}\end{array}\right]\in\mathbb{R}^{19}$ (115)
$\displaystyle\bm{D}^{e}=\left[\begin{array}[]{@{}c|c@{}}\begin{matrix}\bm{D}^{c}_{13\times
13}\end{matrix}&\begin{matrix}\bm{0}_{6\times 6}\\\ \bm{1}_{6\times 6}\\\
\bm{0}_{1\times 6}\\\ \end{matrix}\\\ \hline\cr\bm{0}_{6\times
13}&\bm{0}_{6\times 6}\end{array}\right]\in\mathbb{R}^{19\times 19}$ (118)
$\displaystyle\bar{\bm{H}}^{e}=\left[\begin{array}[]{c}\bar{\bm{H}}^{c}\\\
\hline\cr\bm{0}_{6\times 12}\end{array}\right]\in\mathbb{R}^{19\times 12}$
(121) $\displaystyle\bm{B}^{e}=\left[\begin{array}[]{c}\bm{B}\\\
\hline\cr\bm{0}_{7\times 6}\end{array}\right]\in\mathbb{R}^{19\times 6}$ (124)

where $\bar{\bm{H}}^{c}$ is the nominal value of $\bm{H}^{c}$. The definition
of $\bm{X}^{c}$, $\bm{D}^{c}$, and $\bm{H}^{c}$ can be found in (52). Although
$\bm{u}_{a}$ is considered a part of the state vector in (112), it is just a
residual vector for compensating dynamic uncertainty. Therefore, $\bm{u}_{a}$
is constant in the state space equation and over the horizons. To this end,
the components associated with $\bm{u}_{a}$ in matrices $\bm{D}^{e}$ and
$\bar{\bm{H}}^{e}$ are assigned zero, which means $\dot{\bm{u}}_{a}=0$. Note
that the value of $\bm{u}_{a}$ will be updated according to the adaptive law,
but it is constant during the prediction horizons.

The state representation in (112) is also convenient for discretization
methods such as zero-order hold [56] for MPC. Therefore, our adaptive MPC can
be designed according to (63) and based on the following discrete-time
dynamic:

$\displaystyle\bm{\eta}_{i+1}={\bm{D}^{e}}_{t,i}\bm{\eta}_{t,i}+{\bar{\bm{H}}^{e}}_{t,i}\bm{F}_{i}$
(125)

###  V-C Real-time Computation

The main challenge in executing our proposed adaptive MPC framework is
ensuring that the computation required is fast enough to be performed in real-
time for hardware experiments. If the controller is unable to perform updates
at a high frequency, it could result in the robot collapsing during dynamic
motion. The control system comprises two MPCs, each with 13 to 19 states
predicted over ten horizons. To ensure the robot’s balance and allocate
sufficient computation resources to each control component, we have devised a
scheme, as depicted in Fig. 3c, to update each control component in an
optimized manner.

The robot’s sensory data updates in real-time with a frequency of $1$ $kHz$.
Thus, the reference model should update with the same frequency to compare the
reference model states ($\hat{\bm{X}}$) and real model states ($\bm{X}$)
correctly. The yellow dashed line in Fig. 3c indicates the update frequency
for the reference model. We use the odeint package from Boost software in C++
[57] to solve the ODE problem associated with the dynamic equation for the
reference model.

One of the critical components in our proposed framework is the adaptive MPC,
which is responsible for computing the ground reaction force for the robot, as
shown in Fig. 3c). Through our experimentation, we have determined that for
robust locomotion with dynamic gaits, the optimal update frequency for the
adaptive MPC should be 300 Hz. In contrast, the reference MPC, which plays a
supporting role in the control system, is less sensitive and runs at a slower
rate of 30 Hz. In addition, there is a two-millisecond delay between the
running of the adaptive MPC and reference MPC to ensure sufficient
computational resources are allocated to each component. This means that the
two MPC frameworks do not run simultaneously in our control system.

##  VI Adaptation to Unknown Impact model

The dynamic formulation presented in Sec. IV and Sec. V considers the presence
of model uncertainty in real-world situations. It is assumed that the terrain
is hard enough to allow the robot receives the desired force as ground
reaction forces on its feet. However, this assumption may not hold if the
robot walks on soft or elastic terrain with an unknown impact model, which may
not generate the desired force needed for stable locomotion. Some previous
studies have included terrain knowledge and contact models in their balancing
controllers to address the soft terrain challenge, mainly using a spring-
damper model to characterize the soft terrain [58, 59]. Some control
frameworks for adapting to soft terrain in real-time have also been developed
using iterative learning [60] and whole-body control [61], without prior
knowledge about the terrain. This section demonstrates that the proposed
method in sections Sec. IV and Sec. V can also handle unknown impact models
from terrain, allowing the robot to maintain stability while walking on soft
terrains.

Assume the computed force $\bm{F}$ by MPC in (88) cannot be achieved perfectly
due to walking on soft terrain. Therefore, equation (88) can be rewritten as
follow:

$\displaystyle\bm{\dot{X}}=\bm{D}\bm{X}+\bar{\bm{H}}(\bm{F}_{a}+\tilde{\bm{F}}_{a})+\bm{B}\bm{G}+\bm{B}\bm{\theta}$
(126)

which $\bm{F}_{a}$ is the actual ground reaction force exerted on the robot
and $\tilde{\bm{F}}_{a}$ is the difference between the desired ground reaction
force and actual reaction force. Given that $\tilde{\bm{F}}_{a}$ depends on
the tracking error $\bm{e}$ and time, the uncertainty vector arising from the
ground reaction force can be incorporated with $\bm{\theta}$. Therefore, we
can reformulate equation (126) as follows:

$\displaystyle\bm{\dot{X}}=\bm{D}\bm{X}+\bar{\bm{H}}\bm{F}_{a}+\bm{B}\bm{G}+\bm{B}(\bm{\theta}+\bm{\theta}_{F}).$
(127)

where the uncertainty vector $\bm{\theta}_{F}$ is defined as follow:

$\displaystyle\bm{\theta}_{F}\mathrel{\ensurestackMath{\stackon[1pt]{=}{\scriptstyle\Delta}}}\bm{B}^{T}\bar{\bm{H}}\tilde{\bm{F}}_{a}$
(128)

The equation (127) is in the form of equation (88), which uses actual ground
reaction force instead of desired ground reaction force. Therefore, all
formulations for implementing adaptive controllers are also valid for a
situation with an unknown impact model.

##  VII Results

In this section, we validate our control approach in simulation and hardware
experiments on a Unitree A1 robot. All the hardware experiment’s computation
runs on a single PC (Intel i7-6500U, 2.5 GHz, 64-bit). For simulation, the
control system is implemented in ROS Noetic with the Gazebo 11 simulator,
which provides a high-fidelity simulation of the A1 robot. A video showcasing
the results accompanies this paper111https://youtu.be/QmwyysdTk1k.

We set the control parameters for MPC, the adaption law, and the low-pass
filter as presented in Table I. We use one set of parameters for all the
experiments with different locomotion gaits, indicating that our approach is
easily generalizable. The following subsections will introduce different
experiment results in terms of model and environment uncertainty (see Fig. 5).
In each experiment, the robot starts by using a balance controller to stand up
and then switches to the MPC framework for walking or running.

(a) gravel

(b) grass

(c) uneven terrain

(d) high-sloped terrain

Figure 5: Navigating different terrain using our proposed adaptive MPC while carrying an unknown heavy load. a) gravel, b) grass, c) rough terrain, d) high-sloped terrain. TABLE I: Controller Setting Parameter | Value
---|---
$\bm{Q}$ |  $diag(2.5,2.5,20,0.25,0.25,1.5,$ $0.2,0.2,0.2,0.1,0.1,0.3)$
$\bm{R}$ | $\gamma_{1}\bm{I}_{12}$
$\bm{\Gamma}$ | $diag(1,1,5,2,5,1)\times 10^{3}$
$\omega_{n}$ | 60
$\zeta$ | 0.7

###  VII-A Comparative Analysis

In order to evaluate the performance of our proposed adaptive MPC method, we
conduct a comparative experiment with the conventional MPC method presented in
[1]. The objective is to understand the advantages of integrating the adaptive
controller into MPC for quadrupedal locomotion.

####  VII-A1 Walking with significant model uncertainty

The experiment involves the robot walking and rotating in different
directions, using both adaptive and non-adaptive controllers while carrying an
unknown load. The results of the experiment show that the adaptive controller
provides robust locomotion, with excellent tracking error, even when carrying
an unknown 5 kg load. On the other hand, the non-adaptive controller results
in a considerable error in the COM height and eventually collapses under the
weight of just a 3 kg load. The comparative results for the adaptive and non-
adaptive controllers are shown in Fig. 6.

(a) non-adaptive controller

(b) adaptive controller

(c)

Figure 6: Comparing performance of adaptive and non-adaptive controllers. a)
Snapshots of the A1 robot with the non-adaptive controller while carrying an
unknown 3 kg load and it collapses, b) Snapshots of the A1 robot walking
robustly with the adaptive controller while carrying an unknown 5 kg load, c)
Comparative plots of the COM height for adaptive and non-adaptive controllers.

####  VII-A2 Walking on soft terrain

To evaluate the capability of our proposed control method in handling unknown
impact models, we conducted an experiment where the robot was made for walking
on a double foam, which symbolizes a soft terrain. The performance of both the
adaptive and non-adaptive controllers was evaluated and compared. The results
are depicted in Fig. 7, which represents the robot’s roll angle. The figure
clearly illustrates that the adaptive controller was able to maintain the
robot’s balance on the soft terrain, while the non-adaptive controller was
unable to do so, leading to the collapse of the robot.

(a) non-adaptive controller

(b) adaptive controller

(c)

Figure 7: Comparing performance of adaptive and non-adaptive controllers on
soft terrain. The A1 robot tries to walk on double soft foam using a) non-
adaptive and b) adaptive controllers. c) Shows the plot of the robot’s roll
angle.

###  VII-B Running with Multiple Gaits

To demonstrate the superiority of our proposed approach for dynamic gaits, we
conducted experiments with the robot running while carrying an unknown load.
These experiments were carried out for both the trotting and bounding gaits,
with an unknown load of 5 kg and 3 kg, respectively. The results of these
experiments are shown in 8. It can be seen from the figure that the tracking
of the center of mass height during the bounding gait is more unstable
compared to the trotting gait, which is due to the inherent underactuated
nature of the bounding gait.

(a) trotting

(b) bounding

(c)

Figure 8: Running experiment. The A1 robot runs with the velocity of 1 m/s
using our proposed method. a) trotting gait with an unknown 5 kg load, b)
bounding gait with an unknown 3 kg load, c) Plots of COM height.

###  VII-C Time-varying Load

To demonstrate the effectiveness of our proposed adaptive force control in
adapting to model uncertainty, we conducted simulations where the robot
carries a time-varying load of up to 92% of its weight during walking. As
shown in Fig. 9, our approach can enable the robot to adapt to time-varying
uncertainty. In the simulation, the robot starts with an unknown 5 kg load.
While increasing the robot’s velocity, the robot is subjected to a varying
external force in the z-direction that rises to 60 N, resulting in an
additional unknown 11 kg load. These results indicate that our proposed
approach effectively handles high levels of model uncertainty.

(a)

(b)

(c)

(d)  (e)

Figure 9: Simulation results for the robot carrying a time-varying load. a)
The robot starts with an unknown 5 kg load, then gradually, an unknown time-
varying force will be exerted on the robot as shown in (b) and (c) while the
robot’s velocity increases. d) Plot of COM height, e) Robot velocity tracking
in the x-direction.

###  VII-D Terrain Uncertainty

To demonstrate the capability of our proposed method to handle terrain
uncertainty, we tested the robot navigating various terrain while carrying an
unknown 5 kg load. To this end, we tried walking experiments on multiple rough
terrains as well as high-sloped terrain, and we got impressive results.

####  VII-D1 Rough terrain

We tested the robot navigating various rough terrains such as grass and
gravel. The robot walks and rotates in multiple directions while carrying an
unknown 5 kg load. Some snapshots of the robot walking on diverse rough
terrain are presented in Fig. 5. Our approach is based on a force controller
and retains the robustness features of the baseline framework, allowing the
robot to handle the rough terrain effectively.

####  VII-D2 Sloped terrain

To enable the robot to climb the sloped terrain perfectly without vision, we
need to adjust its orientation to make its body parallel to the walking
surface. This is done by using the footstep location to estimate the slope of
the ground. For each $i$-th leg, we can measure the foot position
$\bm{p}_{i}=(p_{x,i},p_{y,i},p_{z,i})$ and build the vector of feet x-position
($\bm{p}_{x}$), y-position ($\bm{p}_{y}$), and z-position ($\bm{p}_{z}$).
Thus, we can model the walking surface as a plane:

$\displaystyle z(x,y)=a_{0}+a_{1}x+a_{2}y$ (137)

and the coefficients ($a_{0}$, $a_{1}$, and $a_{2}$) will be obtained through
the solution of the least square problem using $\bm{p}_{x}$, $\bm{p}_{x}$, and
$\bm{p}_{x}$ data (see [47] for more details).

Note that the desired roll and pitch angles for the robot will be modified on
the slope according to the following:

$\displaystyle\text{roll}=\arctan(a_{2}),\quad\text{pitch}=\arctan(a_{1}).$
(138)

As a result, the reference model’s desired pitch and roll angles must be
adjusted to the non-zero values determined as described above. It’s important
to note that the reference model utilizes the actual foot position of the
robot, so there is no need to make any changes to the reference model’s
footstep planning when the robot is attempting to climb a slope.

##  VIII Conclusion

In conclusion, a novel control system has been presented that incorporates
adaptive control into force control for legged robots walking under
significant uncertainties. We have demonstrated our proposed approach’s
effectiveness using numerical and experimental validations. The experiments
show the success of the implementation of the proposed adaptive force control
on quadruped robots, allowing them to walk and run while carrying an unknown
heavy load on their trunk. The results are remarkable, with the robot being
able to carry a load of up to 5 kg (50% of its weight) while still keeping the
tracking error within a small range and maintaining stability even in all
directions. The experiment demonstrates that the proposed adaptive force
control system cannot only adapt to model uncertainty but also leverage the
benefits of force control in navigating rough terrains and soft terrain. On
the other hand, the baseline non-adaptive controller fails to track the
desired trajectory and causes the robot to collapse under uncertainty.

## Acknowledgments

The authors would like to thank Yiyu Chen at Dynamic Robotics and Control lab
(DRCL) for his help in conducting the hardware experiments.

###  -A Linear Quadratic Lyapunov Theory

According to Lyapunov theory [62], the PD control described in (78) will
asymptotically stabilize the system if

$\bm{A}_{m}=\begin{bmatrix}\bm{0}_{6}&\bm{1}_{6}\\\
-\bm{K}_{P}&-\bm{K}_{D}\end{bmatrix}\in\mathbb{R}^{12\times 12}$ (139)

is Hurwitz. This means that by choosing a control Lyapunov function candidate
as follows:

$V(\bm{e})=\bm{e}^{T}\bm{P}\bm{e},$ (140)

where $\bm{P}\in\mathbb{R}^{12\times 12}$ is the solution of the Lyapunov
equation

${\bm{A}_{m}}^{T}\bm{P}+\bm{P}\bm{A}_{m}=-\bm{Q}_{L},$ (141)

and $\bm{Q}_{L}\in\mathbb{R}^{12\times 12}$ is any symmetric positive-definite
matrix. We then have:

$\displaystyle\dot{V}(\bm{e},\bm{u})+\lambda V(\bm{e})=$
$\displaystyle~{}\bm{e}^{T}({\bm{D}_{l}}^{T}\bm{P}+\bm{P}{\bm{D}_{l}})\bm{e}$
$\displaystyle+\lambda V(\bm{e})+2\bm{e}^{T}\bm{P}{\bm{B}}\bm{u}~{}\leq 0,$
(142)

where,

$\displaystyle\lambda=\frac{\lambda_{min}(\bm{Q}_{L})}{\lambda_{max}(\bm{P})}>0.$
(143)

As a result, the state variable $\bm{e}$ and the control input $\bm{u}$ always
remain bounded:

$\displaystyle\|\bm{e}\|\leq\delta_{\eta},\quad\|\bm{u}\|\leq\delta_{u}.$
(144)

However, the control signal $\bm{u}^{*}$ (84) we construct by solving QP
problem (II-D), is not always the same as $\bm{u}$. Based on the friction
constraints present in equation (II-D), the value of $\bm{F}^{*}$ is always
bounded. Besides, according to the definition of $\bm{A}$, $\bm{M}$, and
$\bm{G}$, these matrices also have bounded values. Thus, it implies that:

$\displaystyle\|\bm{u}^{*}\|\leq\delta_{{u^{*}}}.$ (145)

Therefore, the vector of difference between $\bm{u}$ and $\bm{u}^{*}$ can be
defined as:

$\displaystyle\bm{\Delta}=\bm{u}^{*}-\bm{u}$ (146)

which is also bounded according to (145) and (144):

$\displaystyle\|\bm{\Delta}\|\leq\delta_{\Delta}.$ (147)

By substituting $\bm{u}^{*}$ in (-A), we have:

$\dot{V}(\bm{e},\bm{u}^{*})+\lambda V(\bm{e})\leq
2\bm{e}^{T}\bm{P}{\bm{B}}\bm{\Delta}\leq\epsilon_{V},$ (148)

where

$\epsilon_{V}=2\|\bm{P}\|\delta_{\eta}\delta_{\Delta}.$ (149)

###  -B Stability Analysis

Theorem: Consider the system dynamics with uncertainty described by (90), and
a reference model described by (91). Assume the use of an $L_{1}$ adaptive
controller with the optimal closed-loop control signal given by (84), the
adaptive control signal given by (95), and the adaptation laws given by (102).
Then, under the aforementioned $L_{1}$ adaptive controller, the tracking error
between the real model and reference model denoted as $\tilde{\bm{e}}$, as
well as the errors between the real and estimated uncertainty, denoted as
$\tilde{\bm{\alpha}}$ and $\tilde{\bm{\beta}}$, respectively, are bounded.

Proof: Let us consider the following control Lyapunov candidate function:

$\displaystyle\tilde{V}=\tilde{\bm{e}}^{T}\bm{P}\tilde{\bm{e}}+\tilde{\bm{\alpha}}^{T}\bm{\Gamma}^{-1}\tilde{\bm{\alpha}}+\tilde{\bm{\beta}}^{T}\bm{\Gamma}^{-1}\tilde{\bm{\beta}}.$
(150)

Therefore, its time derivative will be

$\displaystyle\dot{\tilde{V}}=\dot{\tilde{\bm{e}}}^{T}\bm{P}\tilde{\bm{e}}+\tilde{\bm{e}}^{T}\bm{P}\dot{\tilde{\bm{e}}}+\dot{\tilde{\bm{\alpha}}}^{T}\bm{\Gamma}^{-1}\tilde{\bm{\alpha}}+\tilde{\bm{\alpha}}^{T}\bm{\Gamma}^{-1}\dot{\tilde{\bm{\alpha}}}$
$\displaystyle+\dot{\tilde{\bm{\beta}}}^{T}\bm{\Gamma}^{-1}\tilde{\bm{\beta}}+\tilde{\bm{\beta}}^{T}\bm{\Gamma}^{-1}\dot{\tilde{\bm{\beta}}},$
(151)

in which we have

$\displaystyle\dot{\tilde{\bm{e}}}^{T}\bm{P}\tilde{\bm{e}}+\tilde{\bm{e}}^{T}\bm{P}\dot{\tilde{\bm{e}}}$
$\displaystyle=({\bm{D}_{l}}\tilde{\bm{e}}+{\bm{B}}{\tilde{\bm{F}}})^{T}\bm{P}\tilde{\bm{e}}+\tilde{\bm{e}}^{T}\bm{P}({\bm{D}_{l}}\tilde{\bm{e}}+{\bm{B}}{\tilde{\bm{F}}})$
$\displaystyle~{}~{}+\tilde{\bm{\alpha}}^{T}{\bm{B}}^{T}||\bm{e}||\bm{P}\tilde{\bm{e}}+\tilde{\bm{e}}^{T}\bm{P}{\bm{B}}\tilde{\bm{\alpha}}||\bm{e}||$
$\displaystyle~{}~{}+\tilde{\bm{\beta}}^{T}{\bm{B}}^{T}\bm{P}\tilde{\bm{e}}+\tilde{\bm{e}}^{T}\bm{P}{\bm{B}}\tilde{\bm{\beta}}.$
(152)

Because $\tilde{\bm{e}}=\hat{\bm{e}}-\bm{e}$ satisfies the condition imposed
by (148), it implies that:

$\displaystyle({\bm{D}_{l}}\tilde{\bm{e}}+{\bm{B}}\tilde{\bm{F}})^{T}\bm{P}\tilde{\bm{e}}~{}+$
$\displaystyle\tilde{\bm{e}}^{T}\bm{P}({\bm{D}_{l}}\tilde{\bm{e}}+{\bm{B}}\tilde{\bm{F}})\leq-\lambda\tilde{\bm{e}}^{T}\bm{P}\tilde{\bm{e}}+\epsilon_{\tilde{V}},$
(153)

where

$\epsilon_{\tilde{V}}=2\|\bm{P}\|\delta_{\tilde{e}}\delta_{\tilde{\Delta}}.$
(154)

Furthermore, with the property of the projection operator [54], we have the
following:

$\displaystyle(\hat{\bm{\alpha}}-\bm{\alpha})^{T}(\text{Proj}(\hat{\bm{\alpha}},\bm{y}_{\alpha})-\bm{y}_{\alpha})\leq
0,$
$\displaystyle(\hat{\bm{\beta}}-\bm{\beta})^{T}(\text{Proj}(\hat{\bm{\beta}},\bm{y}_{\beta})-\bm{y}_{\beta})\leq
0.$ (155)

From (102) and (-B), we can imply that

$\displaystyle\tilde{\bm{\alpha}}^{T}\bm{\Gamma}^{-1}\dot{\tilde{\bm{\alpha}}}\leq\tilde{\bm{\alpha}}^{T}\bm{y}_{\alpha}-\tilde{\bm{\alpha}}^{T}\bm{\Gamma}^{-1}\dot{\bm{\alpha}},$
$\displaystyle\tilde{\bm{\beta}}^{T}\bm{\Gamma}^{-1}\dot{\tilde{\bm{\beta}}}\leq\tilde{\bm{\beta}}^{T}\bm{y}_{\beta}-\tilde{\bm{\beta}}^{T}\bm{\Gamma}^{-1}\dot{\bm{\beta}}.$
(156)

We now replace (-B), (-B) and (-B) to (-B), which results in

$\displaystyle\dot{\tilde{V}}$
$\displaystyle\leq-\lambda\tilde{\bm{e}}^{T}\bm{P}\tilde{\bm{e}}+\epsilon_{\tilde{V}}$
$\displaystyle+\tilde{\bm{\alpha}}^{T}(\bm{y}_{\alpha}+{\bm{B}}^{T}\bm{P}\tilde{\bm{e}}||\bm{e}||)-\tilde{\bm{\alpha}}^{T}\bm{\Gamma}^{-1}\dot{\bm{\alpha}}$
$\displaystyle+(\bm{y}_{\alpha}^{T}+\tilde{\bm{e}}^{T}\bm{P}{\bm{B}}||\bm{e}||)\tilde{\bm{\alpha}}-\dot{\bm{\alpha}}^{T}\bm{\Gamma}^{-1}\bm{\alpha}$
$\displaystyle+\tilde{\bm{\beta}}^{T}(\bm{y}_{\beta}+{\bm{B}}^{T}\bm{P}\tilde{\bm{e}})-\tilde{\bm{\beta}}^{T}\bm{\Gamma}^{-1}\dot{\bm{\beta}}$
$\displaystyle+(\bm{y}_{\beta}^{T}+\tilde{\bm{e}}^{T}\bm{P}{\bm{B}})\tilde{\bm{\beta}}-\dot{\bm{\beta}}^{T}\bm{\Gamma}^{-1}\tilde{\bm{\beta}}$
(157)

So, by using the chosen projection functions (IV-C), then we conclude that:

$\displaystyle\dot{\tilde{V}}+\lambda\tilde{V}\leq\epsilon_{\tilde{V}}+\lambda\tilde{\bm{\alpha}}^{T}\bm{\Gamma}^{-1}\tilde{\bm{\alpha}}+\lambda\tilde{\bm{\beta}}^{T}\bm{\Gamma}^{-1}\tilde{\bm{\beta}}$
$\displaystyle-\tilde{\bm{\alpha}}^{T}\bm{\Gamma}^{-1}\dot{\bm{\alpha}}-\dot{\bm{\alpha}}^{T}\bm{\Gamma}^{-1}\tilde{\bm{\alpha}}$
$\displaystyle-\tilde{\bm{\beta}}^{T}\bm{\Gamma}^{-1}\dot{\bm{\beta}}-\dot{\bm{\beta}}^{T}\bm{\Gamma}^{-1}\tilde{\bm{\beta}}.$
(158)

We assume that the uncertainties $\bm{\alpha}$, $\bm{\beta}$, and their time
derivatives are bounded. Furthermore, the projection operators (102) will also
keep $\tilde{\bm{\alpha}}$ and $\tilde{\bm{\beta}}$ bounded (see [34] for a
detailed proof about these properties.) We define these bounds as follows:

$\displaystyle||\tilde{\bm{\alpha}}||\leq$
$\displaystyle\tilde{\bm{\alpha}}_{b},~{}~{}||\tilde{\bm{\beta}}||\leq\tilde{\bm{\beta}}_{b},$
$\displaystyle||\dot{\bm{\alpha}}||\leq$
$\displaystyle\dot{\bm{\alpha}}_{b},~{}~{}||\dot{\bm{\beta}}||\leq\dot{\bm{\beta}}_{b}.$
(159)

Combining this with (-B), we have,

$\displaystyle\dot{\tilde{V}}+\lambda\tilde{V}\leq\lambda\delta_{\tilde{V}},$
(160)

where

$\displaystyle\delta_{\tilde{V}}=2||\bm{\Gamma}||^{-1}(\tilde{\bm{\alpha}}_{b}^{2}+\tilde{\bm{\beta}}_{b}^{2}+\frac{1}{\lambda}\tilde{\bm{\alpha}}_{b}\dot{\bm{\alpha}}_{b}+\frac{1}{\lambda}\tilde{\bm{\beta}}_{b}\dot{\bm{\beta}}_{b})+\frac{1}{\lambda}\epsilon_{\tilde{V}}.$
(161)

Thus, if $\tilde{V}\geq\delta_{\tilde{V}}$ then $\dot{\tilde{V}}\leq 0$. As a
result, we always have $\tilde{V}\leq\delta_{\tilde{V}}$. In other words, by
choosing the adaptation gain $\bm{\Gamma}$ sufficiently large and $\bm{P}$
relatively small, we can limit the Control Lyapunov Function (150) in an
arbitrarily small neighborhood $\delta_{\tilde{V}}$ of the origin. According
to (139) and (141), achieving a small value for $\bm{P}$ depends on choosing a
proper value for $\bm{K}_{P}$, $\bm{K}_{D}$, and $\bm{Q}_{L}$. Therefore, the
value of PD gains affects the stability of the whole system. Finally, the
tracking errors between the dynamics model (90) and the reference model (91),
$\tilde{\bm{e}}$, and the error between the real and estimated uncertainty,
$\tilde{\bm{\alpha}}$, $\tilde{\bm{\beta}}$ are bounded as follows:

$\displaystyle||\tilde{\bm{e}}||\leq\sqrt{\frac{\delta_{\tilde{V}}}{||\bm{P}||}},||\tilde{\bm{\alpha}}||\leq\sqrt{||\bm{\Gamma}||\delta_{\tilde{V}}},||\tilde{\bm{\beta}}||\leq\sqrt{||\bm{\Gamma}||\delta_{\tilde{V}}}.$
(162)

## References

  * [1] J. Di Carlo, P. M. Wensing, B. Katz, G. Bledt, and S. Kim, “Dynamic Locomotion in the MIT Cheetah 3 Through Convex Model-Predictive Control,” in _IEEE Int. Conf. Intell. Robot. Syst._ IEEE, oct 2018, pp. 7440–7447.
  * [2] M. Focchi, A. del Prete, I. Havoutis, R. Featherstone, D. G. Caldwell, and C. Semini, “High-slope terrain locomotion for torque-controlled quadruped robots,” _Auton. Robots_ , vol. 41, no. 1, pp. 259–272, jan 2017.
  * [3] Q. Nguyen, M. J. Powell, B. Katz, J. D. Carlo, and S. Kim, “Optimized Jumping on the MIT Cheetah 3 Robot,” in _2019 Int. Conf. Robot. Autom._ IEEE, may 2019, pp. 7448–7454.
  * [4] H. W. Park, P. M. Wensing, and S. Kim, “Online planning for autonomous running jumps over obstacles in high-speed quadrupeds,” in _Robot. Sci. Syst._ , vol. 11, 2015, pp. 1–9.
  * [5] ——, “High-speed bounding with the MIT Cheetah 2: Control design and experiments,” _Int. J. Rob. Res._ , vol. 36, no. 2, pp. 167–192, feb 2017\.
  * [6] H. Fukushima, T.-H. Kim, and T. Sugie, “Adaptive model predictive control for a class of constrained linear systems based on the comparison model,” _Automatica_ , vol. 43, no. 2, pp. 301–308, feb 2007.
  * [7] V. Adetola, D. DeHaan, and M. Guay, “Adaptive model predictive control for constrained nonlinear systems,” _Syst. Control Lett._ , vol. 58, no. 5, pp. 320–326, may 2009.
  * [8] K. Pereida and A. P. Schoellig, “Adaptive Model Predictive Control for High-Accuracy Trajectory Tracking in Changing Conditions,” _IEEE Int. Conf. Intell. Robot. Syst._ , pp. 7831–7837, 2018.
  * [9] X. Lu, M. Cannon, and D. Koksal-Rivet, “Robust adaptive model predictive control: Performance and parameter estimation,” _Int. J. Robust Nonlinear Control_ , vol. 31, no. 18, pp. 8703–8724, dec 2021.
  * [10] A. Nagabandi, K. Konolige, S. Levine, and V. Kumar, “Deep Dynamics Models for Learning Dexterous Manipulation,” in _Proc. Conf. Robot Learn._ , ser. Proceedings of Machine Learning Research, L. P. Kaelbling, D. Kragic, and K. Sugiura, Eds., vol. 100. PMLR, nov 2020, pp. 1101–1112.
  * [11] B. Amos, I. Dario Jimenez Rodriguez, J. Sacks, B. Boots, and J. Z. Kolter, “Differentiable MPC for End-to-end Planning and Control,” _Adv. Neural Inf. Process. Syst._ , vol. 31, 2018.
  * [12] K. Y. Chee, T. Z. Jiahao, and M. A. Hsieh, “KNODE-MPC: A Knowledge-Based Data-Driven Predictive Control Framework for Aerial Robots,” _IEEE Robot. Autom. Lett._ , vol. 7, no. 2, pp. 2819–2826, apr 2022.
  * [13] S. Chen, B. Zhang, M. W. Mueller, A. Rai, and K. Sreenath, “Learning Torque Control for Quadrupedal Locomotion,” mar 2022. [Online]. Available: https://arxiv.org/pdf/2203.05194
  * [14] G. Bellegarda, Y. Chen, Z. Liu, and Q. Nguyen, “Robust High-Speed Running for Quadruped Robots via Deep Reinforcement Learning,” in _2022 IEEE/RSJ Int. Conf. Intell. Robot. Syst._ IEEE, oct 2022, pp. 10 364–10 370.
  * [15] R. Yang, M. Zhang, N. Hansen, D. H. Xu, and X. Wang, “Learning Vision-Guided Quadrupedal Locomotion End-to-End with Cross-Modal Transformers,” in _Int. Conf. Learn. Represent._ , may 2022.
  * [16] S. Shalev-Shwartz, “Online Learning and Online Convex Optimization,” _Found. Trends Mach. Learn._ , vol. 4, no. 2, pp. 107–194, 2011.
  * [17] L. Bottou and Y. Le Cun, “Large Scale Online Learning,” _Adv. Neural Inf. Process. Syst._ , vol. 16, 2003.
  * [18] A. Nagabandi, C. Finn, and S. Levine, “Deep Online Learning via Meta-Learning: Continual Adaptation for Model-Based RL,” in _7th Int. Conf. Learn. Represent. ICLR_ , 2019.
  * [19] T. Duong and N. Atanasov, “Adaptive Control of SE(3) Hamiltonian Dynamics with Learned Disturbance Features,” _IEEE Control Syst. Lett._ , vol. 6, pp. 2773–2778, 2022.
  * [20] T. Z. Jiahao, K. Y. Chee, and M. A. Hsieh, “Online Dynamics Learning for Predictive Control with an Application to Aerial Robots,” jul 2022. [Online]. Available: http://arxiv.org/abs/2207.09344
  * [21] S. Yang, H. Choset, and Z. Manchester, “Online Kinematic Calibration for Legged Robots,” _IEEE Robot. Autom. Lett._ , vol. 7, no. 3, pp. 8178–8185, jul 2022.
  * [22] S. Zhou, K. Pereida, W. Zhao, and A. P. Schoellig, “Bridging the Model-Reality Gap with Lipschitz Network Adaptation,” _IEEE Robot. Autom. Lett._ , vol. 7, no. 1, pp. 642–649, jan 2022.
  * [23] A. M. Annaswamy and A. L. Fradkov, “A historical perspective of adaptive control and learning,” _Annu. Rev. Control_ , vol. 52, pp. 18–41, jan 2021\.
  * [24] Y. Sun, W. L. Ubellacker, W.-L. Ma, X. Zhang, C. Wang, N. V. Csomay-Shanklin, M. Tomizuka, K. Sreenath, and A. D. Ames, “Online Learning of Unknown Dynamics for Model-Based Controllers in Legged Locomotion,” _IEEE Robot. Autom. Lett._ , vol. 6, no. 4, pp. 8442–8449, oct 2021.
  * [25] K. J. Åström, “Adaptive Control,” in _Math. Syst. Theory Influ. R. E. Kalman_ , A. C. Antoulas, Ed. Berlin, Heidelberg: Springer Berlin Heidelberg, 1991, pp. 437–450.
  * [26] J.-J. E. Slotine and Weiping Li, “On the Adaptive Control of Robot Manipulators,” _Int. J. Rob. Res._ , vol. 6, no. 3, pp. 49–59, sep 1987\.
  * [27] Y.-H. Liu and S. Arimoto, “Decentralized Adaptive and Nonadaptive Position/Force Controllers for Redundant Manipulators in Cooperations,” _Int. J. Rob. Res._ , vol. 17, no. 3, pp. 232–247, mar 1998.
  * [28] Z. Li, S. S. Ge, and Z. Wang, “Robust adaptive control of coordinated multiple mobile manipulators,” _Mechatronics_ , vol. 18, no. 5-6, pp. 239–250, jun 2008.
  * [29] P. Culbertson, J.-J. Slotine, and M. Schwager, “Decentralized Adaptive Control for Collaborative Manipulation of Rigid Bodies; Decentralized Adaptive Control for Collaborative Manipulation of Rigid Bodies,” _IEEE Trans. Robot._ , vol. 37, no. 6, 2021.
  * [30] M. Sombolestan and Q. Nguyen, “Hierarchical Adaptive Loco-manipulation Control for Quadruped Robots,” in _2023 IEEE Int. Conf. Robot. Autom._ IEEE, may 2023, pp. 12 156–12 162.
  * [31] ——, “Hierarchical Adaptive Control for Collaborative Manipulation of a Rigid Object by Quadrupedal Robots,” mar 2023. [Online]. Available: http://arxiv.org/abs/2303.06741
  * [32] J.-J. E. Slotine and W. Li, _Applied nonlinear control_. Prentice hall Englewood Cliffs, NJ, 1991, vol. 199.
  * [33] D. Zhang and B. Wei, “A review on model reference adaptive control of robotic manipulators,” _Annu. Rev. Control_ , vol. 43, pp. 188–198, jan 2017.
  * [34] C. Cao and N. Hovakimyan, “L1 adaptive controller for a class of systems with unknown nonlinearities: Part I,” in _2008 Am. Control Conf._ IEEE, jun 2008, pp. 4093–4098.
  * [35] ——, “Stability Margins of L1 Adaptive Controller: Part II,” in _2007 Am. Control Conf._ IEEE, jul 2007, pp. 3931–3936.
  * [36] ——, “Design and Analysis of a Novel L1 Adaptive Controller, Part II: Guaranteed Transient Performance,” in _2006 Am. Control Conf._ , vol. 2006\. IEEE, 2006, pp. 3403–3408.
  * [37] N. Hovakimyan and C. Cao, _L1 Adaptive Control Theory: Guaranteed Robustness with Fast Adaptation_. SIAM, 2010, no. 5.
  * [38] A. Gahlawat, P. Zhao, A. Patterson, N. Hovakimyan, and E. Theodorou, “L1-GP: L1 Adaptive Control with Bayesian Learning,” in _Proc. 2nd Conf. Learn. Dyn. Control_ , ser. Proceedings of Machine Learning Research, vol. 120\. PMLR, 2020, pp. 826–837.
  * [39] A. Gahlawat, A. Lakshmanan, L. Song, A. Patterson, Z. Wu, N. Hovakimyan, and E. A. Theodorou, “Contraction L1-Adaptive Control using Gaussian Processes,” in _Proc. 3rd Conf. Learn. Dyn. Control_ , ser. Proceedings of Machine Learning Research, vol. 144. PMLR, nov 2021, pp. 1027–1040.
  * [40] G. Tournois, M. Focchi, A. Del Prete, R. Orsolino, D. G. Caldwell, and C. Semini, “Online payload identification for quadruped robots,” _IEEE Int. Conf. Intell. Robot. Syst._ , vol. 2017-Septe, pp. 4889–4896, dec 2017.
  * [41] Q. Nguyen and K. Sreenath, “L1 adaptive control for bipedal robots with control Lyapunov function based quadratic programs,” in _Proc. Am. Control Conf._ , vol. 2015-July. IEEE, jul 2015, pp. 862–867.
  * [42] Q. Nguyen, A. Agrawal, W. Martin, H. Geyer, and K. Sreenath, “Dynamic bipedal locomotion over stochastic discrete terrain,” _Int. J. Rob. Res._ , vol. 37, no. 13-14, pp. 1537–1553, dec 2018.
  * [43] K. Sreenath, H.-W. Park, I. Poulakakis, and J. Grizzle, “Embedding active force control within the compliant hybrid zero dynamics to achieve stable, fast running on MABEL,” _Int. J. Rob. Res._ , vol. 32, no. 3, pp. 324–345, mar 2013.
  * [44] J. W. Grizzle, C. Chevallereau, and C. L. Shih, “HZD-based control of a five-link underactuated 3D bipedal robot,” _Proc. IEEE Conf. Decis. Control_ , pp. 5206–5213, 2008.
  * [45] M. V. Minniti, R. Grandia, F. Farshidian, and M. Hutter, “Adaptive CLF-MPC with application to quadrupedal robots,” _IEEE Robot. Autom. Lett._ , 2021\.
  * [46] M. Sombolestan, Y. Chen, and Q. Nguyen, “Adaptive Force-based Control for Legged Robots,” in _IEEE Int. Conf. Intell. Robot. Syst._ IEEE, sep 2021, pp. 7440–7447.
  * [47] G. Bledt, M. J. Powell, B. Katz, J. Di Carlo, P. M. Wensing, S. Kim, J. D. Carlo, P. M. Wensing, and S. Kim, “MIT Cheetah 3: Design and Control of a Robust, Dynamic Quadruped Robot,” in _IEEE Int. Conf. Intell. Robot. Syst._ Institute of Electrical and Electronics Engineers Inc., dec 2018, pp. 2245–2252.
  * [48] C. Gehring, S. Coros, M. Hutter, M. Bloesch, M. A. Hoepflinger, and R. Siegwart, “Control of dynamic gaits for a quadrupedal robot,” in _2013 IEEE Int. Conf. Robot. Autom._ IEEE, may 2013, pp. 3287–3292.
  * [49] F. Bullo and R. M. Murray, “Proportional Derivative (PD) Control on the Euclidean Group,” jan 1995. [Online]. Available: https://resolver.caltech.edu/CaltechCDSTR:1995.CIT-CDS-95-010
  * [50] R. M. Murray, Z. Li, and S. S. Sastry, _A Mathematical Introduction to Robotic Manipulation_. CRC Press, dec 2017\.
  * [51] M. Raibert, _Legged Robots That Balance_. Cambridge, Mass. MIT Press, 1986.
  * [52] J. Pratt, J. Carff, S. Drakunov, and A. Goswami, “Capture Point: A Step toward Humanoid Push Recovery,” in _2006 6th IEEE-RAS Int. Conf. Humanoid Robot._ IEEE, dec 2006, pp. 200–207.
  * [53] E. R. Westervelt, J. W. Grizzle, and D. E. Koditschek, “Hybrid zero dynamics of planar biped walkers,” _IEEE Trans. Automat. Contr._ , vol. 48, no. 1, pp. 42–56, jan 2003.
  * [54] E. Lavretsky and T. E. Gibson, “Projection Operator in Adaptive Systems,” dec 2011. [Online]. Available: https://arxiv.org/abs/1112.4232
  * [55] D. Q. Mayne, “Model predictive control: Recent developments and future promise,” _Automatica_ , vol. 50, no. 12, pp. 2967–2986, dec 2014.
  * [56] M. S. Fadali and A. Visioli, _Digital Control Engineering: Analysis and Design_. Elsevier Science, 2012.
  * [57] “Boost odeint library.” [Online]. Available: https://www.boost.org/doc/libs/1_64_0/libs/numeric/odeint/doc/html/index.html
  * [58] M. Azad and M. N. Mistry, “Balance control strategy for legged robots with compliant contacts,” in _2015 IEEE Int. Conf. Robot. Autom._ , vol. 2015-June, no. June. IEEE, may 2015, pp. 4391–4396.
  * [59] V. Vasilopoulos, I. S. Paraskevas, and E. G. Papadopoulos, “Monopod hopping on compliant terrains,” _Rob. Auton. Syst._ , vol. 102, pp. 13–26, apr 2018\.
  * [60] A. H. Chang, C. M. Hubicki, J. J. Aguilar, D. I. Goldman, A. D. Ames, and P. A. Vela, “Learning to jump in granular media: Unifying optimal control synthesis with Gaussian process-based regression,” in _2017 IEEE Int. Conf. Robot. Autom._ IEEE, may 2017, pp. 2154–2160.
  * [61] S. Fahmi, M. Focchi, A. Radulescu, G. Fink, V. Barasuol, and C. Semini, “STANCE: Locomotion Adaptation Over Soft Terrain,” _IEEE Trans. Robot._ , vol. 36, no. 2, pp. 443–457, apr 2020.
  * [62] A. D. Ames, K. Galloway, K. Sreenath, and J. W. Grizzle, “Rapidly Exponentially Stabilizing Control Lyapunov Functions and Hybrid Zero Dynamics,” _IEEE Trans. Automat. Contr._ , vol. 59, no. 4, pp. 876–891, apr 2014.

|  Mohsen Sombolestan received his B.Sc. degree in mechanical engineering in
2017 from Sharif University of Technology, Tehran, Iran, and his M.Sc. degree
in mechanical engineering in 2020 from Isfahan University of Technology,
Isfahan, Iran. He is working toward a Ph.D. in mechanical engineering from
University of Southern California, Los Angeles, CA, USA. His research
interests include control system design in robotic applications, especially
legged robots, focusing on adaptive control and reinforcement learning.
---|---
|  Quan Nguyen is an assistant professor of Aerospace and Mechanical
Engineering at the University of Southern California (USC). Before joining
USC, he was a Postdoctoral Associate in the Biomimetic Robotics Lab at the
Massachusetts Institute of Technology (MIT). He received his Ph.D. from
Carnegie Mellon University (CMU) in 2017 with the Best Dissertation Award. His
research interests span different control and optimization approaches for
highly dynamic robotics, including nonlinear control, trajectory optimization,
real-time optimization-based control, and robust and adaptive control. His
work on the MIT Cheetah 3 robot leaping on a desk was featured widely in many
major media channels, including CNN, BBC, NBC, ABC, etc. Nguyen won the Best
Presentation of the Session at the 2016 American Control Conference (ACC) and
the Best System Paper Finalist at the 2017 Robotics: Science & Systems
Conference (RSS).
---|---