# Stability Guarantees for Continuous RL Control

Bing Song1, Jean-Jacques Slotine2 and Quang-Cuong Pham3 1Bing Song is with HP-
NTU Digital Manufacturing Corporate Lab, Singapore, Singapore 639798,
<EMAIL_ADDRESS>Slotine is with Faculty of Mechanical
Engineering, Massachusetts Institute of Technology, Cambridge MA 02139-4307,
<EMAIL_ADDRESS>Pham is with HP-NTU Digital Manufacturing Corporate
Lab, Singapore, Singapore 639798<EMAIL_ADDRESS>
###### Abstract

Lack of stability guarantees strongly limits the use of reinforcement learning
(RL) in safety critical robotic applications. Here we propose a control system
architecture for continuous RL control and derive corresponding stability
theorems via contraction analysis, yielding constraints on the network weights
to ensure stability. The control architecture can be implemented in general RL
algorithms and improve their stability, robustness, and sample efficiency. We
demonstrate the importance and benefits of such guarantees for RL on two
standard examples, PPO learning of a 2D problem and HIRO learning of maze
tasks.

##  I Introduction

Reinforcement learning (RL) is expected to enable robots to learn multi-level
reasoning tasks beyond the capability of classic control. However, lack of
stability guarantees has prevented RL from stable performance comparable to
classic control and hence wide applications in robotics, especially safety
critical applications.

Here we propose a control architecture (Fig. 1) for continuous RL control and
derive stability theorems for this architecture via contraction analysis [1].
Accordingly stability can be guaranteed by constraining the network weights.
This architecture can be implemented in general RL algorithms to improve
stability, robustness, and sample efficiency. For demonstration, we provide
two examples: PPO [2] learning a 2D math problem and HIRO [3] learning maze
tasks (Fig. 1).

Figure 1: Control system architecture for continuous RL control and multi-
level tasks for demonstration. The control block diagram (a) is used to derive
stability theorems via contraction analysis, which includes a coordinate
transformation $W_{y}$ and $W_{u}$ as well as an inner policy $\pi_{z}$. We
tested HIRO with stability guarantee via an ant maze (b) and a contact maze
task (c).

A dynamic model $\dot{\mathbf{y}}=\mathbf{f}(\mathbf{y},\mathbf{u},t)$
(specifically $\partial\mathbf{f}/\partial\mathbf{y}$ and
$\partial\mathbf{f}/\partial\mathbf{u}$) is required to compute the stability
constraints. When this model only represents the robot dynamics, it can be
obtained analytically or from system identification. Given known models, the
implementation method for stability guarantees is in Section III.

To learn multi-level tasks with stability guarantees for low-level control,
the model $\dot{\mathbf{y}}=\mathbf{f}(\mathbf{y},\mathbf{u},t)$ represents
not only the robot dynamics but also a high-level policy as well as other
dynamics. In this case, the model is changing as learning progresses because
of an updating high-level policy. It is infeasible to learn such a model by
system identification simultaneously. Thus we propose to leverage impedance
control techniques to make the end effector behave like a spring-mass-damper
system in the Cartesian space. Based on this spring-mass-damper model made of
controller gains, we can apply stability analysis. This can be used in low-
level control for hierarchical RL to learn multi-level tasks and improve its
sample efficiency. The control block diagram of this method is plotted in Fig.
2 and the corresponding stability theorem is in Section IV.

To demonstrate our methods, we provide two examples: PPO learning a 2D math
problem and HIRO learning maze tasks in Section V.

In the PPO example, we seek the answer to the question: does PPO learn a
nonlinear controller that is stable in the sense of classic control for this
simple problem or a curve fitting that cannot be extended for situations
beyond training (e.g., longer running time for control)? Results show that PPO
can learn stable control for this problem but very sensitive to neural network
sizes and training time. Deteriorating performance is observed as training
progresses after the initial convergence. With stability guarantees, PPO can
maintain the stable performance regardless of training time as well as with a
large range of network sizes. PPO with stability guarantees also shows better
robustness to model errors.

In the HIRO example, we want to examine if sample efficiency can be improved
by constraining the low-level exploration inside the safe space that satisfies
stability theorems. For the two maze tasks, HIRO with low-level stability
guarantees outperforms the already-efficient original HIRO in sample
efficiency. We also observed that the performance is very sensitive to the
control architecture (Fig. 2) and we will perform an ablative study with
respect to the architecture and analyze the learning process.

Figure 2: Control block diagram for rigid body dynamics learning multi-level
tasks. For multi-level tasks, the 1st-order nonlinear model represents not
only the robot dynamics but also the high-level policy and so on. We leverage
Cartesian space impedance control as well as other task-dependent control
techniques to make the end effector behave like a spring-mass-damper system,
i.e.,
$\Lambda_{d}\ddot{\mathbf{e}}+K_{d}{\dot{\mathbf{e}}}+K_{p}{\mathbf{e}}=-\mathbf{u}$,
where $\mathbf{e}=\mathbf{x}_{d}-\mathbf{x}$ and $\mathbf{u}$ is the low-level
action. By defining a composite variable
$\mathbf{y}(t)=K_{1}\mathbf{e}(t)+K_{2}\dot{\mathbf{e}}(t)$, we can obtain a
1st-order model $\dot{\mathbf{y}}=\mathbf{f}(\mathbf{y},\mathbf{u},t)$ to
compute the stability constraints.

Existing studies in stability usually use Lyapunov theory [4]. Here we use its
differential counterpart: contraction theory [1]. The advantage of contraction
analysis over Lyapunov is that it avoids to find a proper Lyapunov function
and provides explicit sufficient conditions for stability in terms of
generalized Jacobians.

Besides Lyapunov, some works have studied stability for specific systems or
controllers [5, 6], which are less general than studies based on Lyapunov
theory or contraction theory.

Contraction theory has been used in the research area of learning for system
identification [7, 8, 9] or as a guidance [10, 11, 12] via Control Contraction
Metrics [13, 14, 15, 16] and its extensions, which are also applied in
nonlinear control for stability [17, 18]. Unlike those works searching for
stable solutions via Control Contraction Metrics, we solve for the safe set of
networks with respect to a given metric.

Instead of recursively searching in a safe set of networks, we propose to
limit the control network weights for stability, at the cost of complex in the
control architecture in Fig. 1, which (1) makes the mapping from $\mathbf{y}$
to $\mathbf{u}$ (i.e.,$\pi$) a first-order nonlinear equation by adding an
integral into the input to the control networks and (2) decouples the system
dynamics via the coordinate transformation $W_{y}$ and $W_{u}$ that allows
contraction analysis to be applied in each subspace.

To our best knowledge, this is the first time that applies contraction
analysis for continuous RL control.

##  II Stability theorem via contraction analysis

We derive the stability theorem based on the nonlinear model,

$\displaystyle\dot{\mathbf{y}}=\mathbf{f}(\mathbf{y},\mathbf{u},t)$ (1)
$\displaystyle\mathbf{u}=\pi(\mathbf{y})$

with respect to
$\begin{bmatrix}\mathbf{y}^{T}&\mathbf{u}^{T}\end{bmatrix}^{T}$, where
$\mathbf{y}\in\mathcal{R}^{m}$ is the RL state; $\mathbf{u}\in\mathcal{R}^{m}$
is the mean of the action distribution; and $\pi$ denotes the control policy
without stochastic sampling for exploration. We assume that with proper
sampling strategies, the system remains stable with sampled actions if it is
stable with $\mathbf{u}$.

Stability of nonlinear systems is generally defined relative to some nominal
motion or equilibrium. Without requiring to know what the nominal motion is,
contraction theory analyzes if a system is converging to the nominal motion
regardless of initial conditions or temporary disturbances [1]. To simply put,
if all trajectories converge to some unknown nominal motion in a region, this
region is a contraction region; a system
$\dot{\mathbf{x}}=\mathbf{f}(\mathbf{x},t)$ is contracting in a region, if its
generalized Jacobian is uniformly negative in this region. The generalized
Jacobian refers to the Jacobian with respect to some metric
$\mathbf{M}(\mathbf{x},t)$ that is symmetric and continuously differentiable.
Details are in Appendix VII-A.

With a control policy $\pi$ of a particular architecture plotted in Figure 1,
which includes

  * •

a differential coordinate transformation

$\mathbf{z}=W_{y}\mathbf{y},\ \ \ \mathbf{u}=W_{u}^{-1}\mathbf{v}$ (2)

  * •

an inner policy

$\mathbf{v}=\pi_{z}(\mathbf{s}_{z})$ (3)

  * •

the integral of $\mathbf{z}$ as the input to the inner policy which is
approximated by $\sum_{0}^{k-1}\mathbf{z}\delta_{t}$, where $\delta_{t}$ is
the sampling period for control

$\mathbf{s}_{z}=\begin{bmatrix}\mathbf{s}_{1}\\\
\mathbf{s}_{2}\end{bmatrix}=\begin{bmatrix}\mathbf{z}\\\
\sum_{0}^{k-1}\mathbf{z}\delta_{t}\end{bmatrix}$ (4)

we derive sufficient conditions for a system to be contracting in a region,
i.e., stable in that region.

Theorem 1. Given a nonlinear system Eq. (1) and a control policy of the
architecture in Fig. 1, if in a region

  * •

$W_{y}$ and $W_{u}$ are invertible that can form an uniformly positive
definite metric $\mathbf{M}=\mathbf{\Theta}^{T}\mathbf{\Theta}$ where

$\mathbf{\Theta}=\begin{bmatrix}W_{y}(\mathbf{y},\mathbf{u},t)&\\\
&W_{u}(\mathbf{y},\mathbf{u},t)\end{bmatrix}$ (5)

  * •

there exists $\epsilon>0$ such that

$(F_{1}^{T}+F_{1})\prec-[\epsilon+max(\nu^{+},0)]I$ (6)

where $\nu^{+}$ is the upper bound of eigenvalues of $F_{2}^{T}+F_{2}$,

$\displaystyle
F_{1}=\begin{bmatrix}W_{y}\frac{\partial\mathbf{f}}{\partial\mathbf{y}}W_{y}^{-1}&W_{y}\frac{\partial\mathbf{f}}{\partial\mathbf{u}}W_{u}^{-1}\\\\[8.0pt]
\frac{\partial\pi_{z}}{\partial\mathbf{s}_{1}}W_{y}\frac{\partial\mathbf{f}}{\partial\mathbf{y}}W_{y}^{-1}+\frac{\partial\pi_{z}}{{\partial\mathbf{s}_{2}}}&\frac{\partial\pi_{z}}{\partial\mathbf{s}_{1}}W_{y}\frac{\partial\mathbf{f}}{\partial\mathbf{u}}W_{u}^{-1}\end{bmatrix}$
(7) $\displaystyle F_{2}=\begin{bmatrix}\dot{W}_{y}W_{y}^{-1}&0\\\\[8.0pt]
\frac{\partial\pi_{z}}{\partial\mathbf{s}_{1}}\dot{W}_{y}W_{y}^{-1}&0\end{bmatrix}$

this system Eq. (1) is contracting in that region.

Proof. The system dynamic Eq. (1) yields the differential relation

$\delta\dot{\mathbf{y}}=\frac{\partial\mathbf{f}}{\partial\mathbf{y}}(\mathbf{y},\mathbf{u},t)\delta\mathbf{y}+\frac{\partial\mathbf{f}}{\partial\mathbf{u}}(\mathbf{y},\mathbf{u},t)\delta\mathbf{u}$
(8)

where $\delta\mathbf{y}$ and $\delta\mathbf{u}$ are differential
displacements, i.e., infinitesimal displacements at fixed time. The
transformation Eq. (2) yields the differential relation

$\displaystyle\delta\dot{\mathbf{z}}=\dot{W}_{y}\delta\mathbf{y}+W_{y}\delta\dot{\mathbf{y}}$
(9)
$\displaystyle\delta\dot{\mathbf{u}}=\dot{W}_{u}^{-1}\delta\mathbf{v}+W_{u}^{-1}\delta\dot{\mathbf{v}}$

The inner policy $\pi_{z}$ Eq. (3) yields

$\delta\dot{\mathbf{v}}=\frac{\partial\pi_{z}}{\partial\mathbf{s}_{1}}\delta\dot{\mathbf{s}}_{1}+\frac{\partial\pi_{z}}{{\partial\mathbf{s}_{2}}}\delta\mathbf{s}_{1}$
(10)

because $\delta\dot{\mathbf{s}}_{2}=\delta\mathbf{s}_{1}$. Combining Eq. (8),
Eq. (9), Eq. (10), and $\mathbf{z}=\mathbf{s}_{1}$, we obtain the equation in
the Euclidean metric,

$\begin{bmatrix}\delta\dot{\mathbf{y}}\\\
\delta\dot{\mathbf{u}}\end{bmatrix}=J\begin{bmatrix}\delta\mathbf{y}\\\
\delta\mathbf{u}\end{bmatrix}=\begin{bmatrix}\frac{\partial\mathbf{f}}{\partial\mathbf{y}}&\frac{\partial\mathbf{f}}{\partial\mathbf{u}}\\\\[8.0pt]
J_{21}&J_{22}\end{bmatrix}\begin{bmatrix}\delta\mathbf{y}\\\
\delta\mathbf{u}\end{bmatrix}$ (11)

where

$\displaystyle
J_{21}=W_{u}^{-1}\frac{\partial\pi_{z}}{\partial\mathbf{s}_{1}}(\dot{W}_{y}+W_{y}\frac{\partial\mathbf{f}}{\partial\mathbf{y}})+W_{u}^{-1}\frac{\partial\pi_{z}}{{\partial\mathbf{s}_{2}}}W_{y}$
(12) $\displaystyle
J_{22}=W_{u}^{-1}\frac{\partial\pi_{z}}{\partial\mathbf{s}_{1}}W_{y}\frac{\partial\mathbf{f}}{\partial\mathbf{u}}+\frac{d(W_{u}^{-1})}{dt}W_{u}$

and the equation in the $M=\Theta^{T}\Theta$ metric,

$\begin{bmatrix}\delta\dot{\mathbf{z}}\\\
\delta\dot{\mathbf{v}}\end{bmatrix}=F\begin{bmatrix}\delta\mathbf{z}\\\
\delta\mathbf{v}\end{bmatrix}$ (13)

where

$F=\begin{bmatrix}[\dot{W}_{y}+W_{y}\frac{\partial\mathbf{f}}{\partial\mathbf{y}}]W_{y}^{-1}&W_{y}\frac{\partial\mathbf{f}}{\partial\mathbf{u}}W_{u}^{-1}\\\\[8.0pt]
F_{21}&\frac{\partial\pi_{z}}{\partial\mathbf{s}_{1}}W_{y}\frac{\partial\mathbf{f}}{\partial\mathbf{u}}W_{u}^{-1}\end{bmatrix}$
(14)

and

$F_{21}=\frac{\partial\pi_{z}}{\partial\mathbf{s}_{1}}(\dot{W}_{y}+W_{y}\frac{\partial\mathbf{f}}{\partial\mathbf{y}})W_{y}^{-1}+\frac{\partial\pi_{z}}{{\partial\mathbf{s}_{2}}}$
(15)

This $F=(\dot{\Theta}+\Theta J)\Theta^{-1}$ is the generalized Jacobian of Eq.
(1) in the $M$ metric. According to contraction theory [1] (see Appendix
VII-A), if $F$ is uniformly negative definite in some region, the system Eq.
(1) is contracting in the region.

We rewrite $F$ into the sum of the constant $F_{1}$ and the time-varying
$F_{2}$. We denote the maximum eigenvalues of $F_{2}^{T}+F_{2}$,
$F_{1}^{T}+F_{1}$, and $F^{T}+F$ by $\nu^{+}$, $\lambda^{+}$, and $\gamma^{+}$
respectively. Because $\gamma^{+}<\nu^{+}+\lambda^{+}$, given $\epsilon>0$, if
$\lambda^{+}<-(\epsilon+\nu^{+})$, $\gamma^{+}<-\epsilon$. Hence, a sufficient
condition for uniformly negative definite $F$ is that there exists
$\epsilon>0$ such that $\lambda^{+}<-(\epsilon+max(\nu^{+},0))$ in that
region. Proof is finished.

According to Theorem 1, given a coordinate transformation and the system
dynamic model, we can constrain the Jacobian of the inner policy
($\partial\pi_{z}/\partial\mathbf{s}_{1}$ and
$\partial\pi_{z}/\partial\mathbf{s}_{2}$) for stability.

##  III Stability guarantee for known systems

According to Theorem 1, we need to solve Eq. (6) to find the constraints of
$\partial\pi_{z}/\partial\mathbf{s}_{1}$ and
$\partial\pi_{z}/\partial\mathbf{s}_{2}$ for stability. Here we specify
$W_{y}$, $W_{u}$ and $\pi_{z}$ to make the problem algebraically solvable. The
first step is to transform Eq. (13) into the form, for example,

$\displaystyle\begin{bmatrix}\delta\dot{z}_{1}\\\ \delta\dot{v}_{1}\\\
\hline\cr\delta\dot{z}_{2}\\\ \delta\dot{v}_{2}\\\
\hline\cr\delta\dot{z}_{3}\\\
\delta\dot{v}_{3}\end{bmatrix}=\begin{bmatrix}\begin{array}[]{c|c|c}A_{1}&B_{12}&B_{13}\\\
\hline\cr 0&A_{2}&B_{23}\\\ \hline\cr 0&0&A_{3}\\\
\end{array}\end{bmatrix}\begin{bmatrix}\delta z_{1}\\\ \delta v_{1}\\\
\hline\cr\delta z_{2}\\\ \delta v_{2}\\\ \hline\cr\delta z_{3}\\\ \delta
v_{3}\end{bmatrix}$ (16)

For a system of the above form, if $A_{i}$’s are uniformly negative and
$B_{ij}$’s are bounded, the system Eq. (1) is contracting [1]. By solving the
characteristic functions of $A_{i}$’s, we obtain the boundaries of network
weights for stability. This is Theorem 2 and its proof is in Appendix VII-B.

Theorem 2. With the following designs,

  * •

$W_{y}$: the normalized eigenvector matrix of
$\partial\mathbf{f}/\partial\mathbf{y}$

  * •

$W_{u}$: $W_{u}=PQ^{T}$ where $Q$ is from the QR decomposition of
$P(W_{y}\partial\mathbf{f}/\partial\mathbf{u})$ and $P$ is a permutation
matrix with all ones at the skew diagonal

  * •

$\pi_{z}$: a parallel sequence of independent neural networks denoted by
$h_{i}$ as in Figure. 3

Figure 3: The inner policy $\pi_{z}$ inside the control policy $\pi$ (Fig. 1)
is made of a parallel sequence of independent neural networks $h_{i}$, the
inputs of which are
$\begin{bmatrix}s_{1i}&s_{2i}\end{bmatrix}=\begin{bmatrix}z_{i}&\sum_{0}^{k-1}z_{i}\delta_{t}\end{bmatrix}$
where $i=1,2,\cdots,m$.

and with the assumptions that $\dot{W}_{y}$ can be ignored and the partial
derivatives of $\pi_{z}$ are bounded, for a system Eq. (1) with diagonalizable
$\partial\mathbf{f}/\partial\mathbf{y}$, if in a region there exists
$\epsilon>0$ that

$\displaystyle\frac{\partial h_{i}}{\partial
s_{1i}}b_{ii}+\lambda_{i}<-\epsilon,\ \ \ \frac{\partial h_{i}}{\partial
s_{2i}}b_{ii}<-\epsilon$ (17)

where $\lambda_{i}$ is the $i$th eigenvalue of
$\partial\mathbf{f}/\partial\mathbf{y}$, $b_{ii}$ is the $i^{th}$ diagonal
component of $W_{y}(\partial\mathbf{f}/\partial\mathbf{u})W_{u}^{-1}$, the
system Eq. (1) is contracting in that region.

Given the bounds of $\lambda_{i}$ and $b_{ii}$, Theorem 2 provides inequality
constraints with respect to $\partial h_{i}/\partial s_{ij}$ ($i=1,2$ and
$j=1,2,\cdots,m$) .

Here we propose a method of implementing Theorem 2 for a general RL algorithm
(Note that implementation of the control architecture does not adapt the
learning algorithm):

  * •

Apply control techniques to make the system stable when the RL actions are
equal to zero, i.e. $\dot{\mathbf{y}}=\mathbf{f}(\mathbf{y},0,t)$ converging
to some equilibrium as time progresses. This means $\lambda_{i}\leq 0$
otherwise the solution to the differential equation
$\delta\dot{\mathbf{y}}=\partial\mathbf{f}/\partial{\mathbf{y}}\delta\mathbf{y}$
contains a diverging component $e^{\lambda_{i}t}$ which is contradict to the
premise that the system is converging. This is generally true for robotic
systems in real world.

  * •

Use a control policy of the architecture in Fig. 1 and use the inner policy of
the architecture in Fig. 3 with multilayer perceptron networks (MLPs) (e.g.,
with tanh as activation functions) as $h_{i}$’s ($i=1,2,\cdots,m$). In this
case, $\partial h_{i}/\partial s_{ij}$ ($i=1,2$ and $j=1,2,\cdots,m$) becomes
the product of the network weights scaled by the derivative of the activation
function.

  * •

Set the weights of an $l$-layered network $h_{i}$ as follows. For the $j^{th}$
layer ($j=1,2,\cdots,l-1$), the weights are bigger than some small positive
number $\epsilon$. For the last layer, the weights are bigger than $-\epsilon
b_{ii}/\lVert b_{ii}\rVert$ if $b_{ii}<0$ or smaller than $-\epsilon
b_{ii}/\lVert b_{ii}\rVert$ if $b_{ii}>0$. In code, this can be implemented by
constraining the weights by $\epsilon$ and adding another operation that
multiplies the opposite sign of $b_{ii}$.

With this implementation and the assumptions in Theorem 2, the system is
contracting except at the points with a saturating activation function. Proof
is in Appendix VII-C. With a saturating activation function, the RL action
becomes a constant $\mathbf{u}_{const}$ and
$\dot{\mathbf{y}}=\mathbf{f}(\mathbf{y},\mathbf{u}_{const},t)$, which requires
pre-designed control techniques and a proper range of RL actions for
stability.

When there exists model errors in $b_{ii}$ and $\lambda_{i}$, as long as the
inequalities Eq. (17) hold, the system remains stable.

##  IV Stability guarantee for unknown systems: rigid body dynamics learning
multi-level tasks

For rigid body dynamics to learn multi-level tasks via hierarchical RL, it is
infeasible to learn the nonlinear model that includes the high level policy.
Thus we leverage Cartesian-space impedance control to make the end-effectors
behave like a spring-mass-damper system,

$\Lambda_{d}\ddot{\mathbf{e}}+K_{d}{\dot{\mathbf{e}}}+K_{p}{\mathbf{e}}+K_{rl}\mathbf{u}_{rl}=0$
(18)

where $\mathbf{e}=\mathbf{x}_{d}-\mathbf{x}$ denotes the tracking error of the
end effectors, $\mathbf{x}_{d}$ is the high-level RL action, $\Lambda_{d}$,
$K_{d}$, and $K_{p}$ are the controller gains that are positive definite
diagonal matrices, $K_{rl}$ is the weights for the feed-forward low-level RL
actions. We used techniques from [19, 20, 21] to build the control diagram in
Fig. 2 (see Appendix VII-D).

We then use a composite variable,

$\mathbf{y}=K_{1}\mathbf{e}+K_{2}\dot{\mathbf{e}}$ (19)

to decrease the order of Eq. (18). Plugging Eq. (19) into Eq. (18), the
dynamic model becomes

$\Lambda_{y}\dot{\mathbf{y}}+B_{y}\mathbf{y}+K_{rl}\mathbf{u}_{rl}=0$ (20)

where $K_{1}$, $K_{2}$, $\Lambda_{y}$ and $B_{y}$ are positive definite
diagonal matrices made of $\Lambda_{d}$, $K_{d}$, and $K_{p}$:

  * •

When $K_{p}=0$, we use $K_{2}=K_{d}$ and $K_{1}=0$ which results in

$\Lambda_{y}=\Lambda_{d},\ \ \ B_{y}=K_{d}$ (21)

  * •

Otherwise we use

$K_{1,i}=\frac{K_{p,i}}{\lambda_{d,i}},\ \ \ \
K_{2,i}=\frac{2K_{p,i}}{K_{d,i}\pm\sqrt{K_{d,i}^{2}-4K_{p,i}\lambda_{d,i}}}$
(22)

where $K_{1,i}$, $K_{2,i}$, $\lambda_{d,i}$, $K_{p,i}$, and $K_{d,i}$ denotes
the $i^{th}$ diagonal element of $K_{1}$, $K_{2}$, $\Lambda_{d}$, $K_{p}$, and
$K_{p}$ respectively. Note this requires
$K_{d,i}^{2}-4K_{p,i}\lambda_{d,i}\geq 0$. Matrices $K_{1}$ and $K_{2}$ are
positive definite given stable controllers with positive definite
$\Lambda_{d}$, $K_{d}$ and $K_{p}$.

The coefficient matrices in Eq. (20) becomes

$\Lambda_{y}=\Lambda_{d}K_{2}^{-1},\ \ \ B_{y}=\Lambda_{d}$ (23)

are positive definite diagonal matrices since $\Lambda_{d}$, $K_{2}$ are
positive definite diagonal matrices.

This linear model Eq. 20 can be used as the system model in Theorem 2.
However, this model Eq. 20 suffers from nonlinear errors (denoted by
$\bm{\zeta}$),

$\Lambda_{y}\dot{\mathbf{y}}+B_{y}\mathbf{y}+K_{rl}\mathbf{u}_{rl}=\bm{\zeta}(\mathbf{q},\dot{\mathbf{q}},\ddot{\mathbf{q}})$
(24)

Hence we propose to use absorbing matrices to evaluate the robustness. If
there exists time-varying diagonal matrices $\tilde{\Lambda}_{y}\succ 0$ and
$\tilde{B}_{y}\succ 0$ that satisfies

$\bm{\zeta}(\mathbf{q},\dot{\mathbf{q}},\ddot{\mathbf{q}})+(\tilde{\Lambda}_{y}-\Lambda_{y})\dot{\mathbf{y}}+(\tilde{B}_{y}-B_{y}){\mathbf{y}}=0$
(25)

where $\Lambda_{y}$ and $B_{y}$ are positive definite matrices in Eq. (23),
the system remains stable.

We summarize the above discussions into Theorem 3 and its proof is in Appendix
VII-E. The implementation is to first apply the control architecture in Fig. 2
and then follow the implementation method of Theorem 2.

Theorem 3 For rigid body dynamics to learn multi-level tasks via hierarchial
RL, with the control diagram in Fig. 2 implemented for low-level control, with
the assumption that the partial derivatives of $\pi_{z}$ are bounded, if in a
region

  * •

there exists $\tilde{\Lambda}_{y}\succ 0$ and $\tilde{B}_{y}\succ 0$ that
satisfies Eq. (25)

  * •

there exists $\epsilon>0$ that

$\displaystyle\frac{\partial h_{i}}{\partial z_{i}}K_{rl,i}>\epsilon,\ \ \
\frac{\partial h_{i}}{\partial({\sum_{0}^{k-1}z_{i}})}K_{rl,i}>\epsilon$ (26)

where $K_{rl,i}$ is the $i$th diagonal component of $K_{rl}$

the system with respect to the composite variable Eq. (24) is contracting.

By solving Eq. (19) with positive definite $K_{1}$ and $K_{2}$, the tracking
error is converging to zero with some transients when the composite variable
is converging to zero, that is, asymptotically stable.

##  V Experiments

We demonstrate our methods on two algorithms: (a) PPO learning a 2D problem
and (b) HIRO for rigid body dynamics learning maze tasks (Fig. 1 (b) and (c)).
We tested 5 different seeds for training. Environment and implementation
details are in Appendix VII-F.

###  V-A PPO vs PPO with stability guarantee (C-PPO)

In the 2D example, with known $\partial\mathbf{f}/\partial\mathbf{y}$ and
$\partial\mathbf{y}/\partial\mathbf{u}$, we apply Theorem 2 via the
implementation method in Section III for stability guarantee.

The 2D example approximates a peg with 1st-order dynamics moving in a
perpendicular plane to touch $x_{d}$ with force $f_{d}$ ($x_{d}$ and $f_{d}$
are randomly sampled),

$\displaystyle\tau_{z}\dot{z}+z=u_{z}$ (27)
$\displaystyle\tau_{x}\dot{x}+x=u_{x}$ $\displaystyle f=K_{sur}\min{(z-g(x),\
0)}$

where the surface profile $g(x)=K_{1}\sin(x)+K_{2}\cos(x)$. The reward is a
weighted sum of penalties on negative distance and chattering motions.

To answer the question, does PPO learn a nonlinear controller that is stable
in the sense of classic control or a curve fitting that can not be extended
for situations beyond training (e.g., longer running time for control)? We
design the experiments as follows:

  * •

PPO is considered to have learned some stable nonlinear controller (i.e., a
stable learned policy) if the tracking accuracy can be kept for extended
controller running time in the tests. The control policy was running 8 s in
training (8 second trajectories) but tested for 16 s. The target was sampled
in a range that can be reached within a short time much less than 8 s. We
checked the resulting trajectories and most of them have reached the goals in
2 s. This implies that 8 s is reasonably long enough to learn this math
problem.

  * •

Different neural network sizes, different definitions of RL states, and other
parameters are tested to find the best performance of PPO without stability
guarantee. The best performance is on the same level as PPO with stability
guarantee (named C-PPO). The settings with the best performance are used for
the stability and robustness test.

Stability test. After the initial convergence of PPO and C-PPO (learning
curves of episodes vs negative expected returns are plotted in Fig. 4 (c)), we
tested the learned policy every 1k episodes, from Episode 4000 until 9000. An
increasing percentage of failed trajectories (error larger than 0.4) is
observed for PPO (Fig. 4 (a)) as training progresses. In PPO, We observed
three types of failed trajectories:

  * •

Chattering motions Fig. 4 (b) and (d). Chattering motions are related to
$\partial\mathbf{f}/\partial\mathbf{u}$ of Eq. (1), the value of which changes
on a large scale for the task.

  * •

Non-zero equilibrium Fig. 4 (f). To avoid non-zero equilibrium, we propose
Theorem 4 in Appendix VII-F.

  * •

Drifting away Fig. 4 (e). Some trajectories drift away after 8 s, which
implies that the networks have learned to cancel the errors for each data
point of a trajectory but failed to learn a stable controller.

C-PPO maintains the tracking accuracy in the stability test. With stability
guarantees, we limit the network to approximate a stable state-dependent PI
controller instead of unlimited nonlinear curve fittings to cancel the errors.

Robustness test. We tested robustness of the learned policy at Episode 4000
for PPO and Episode 5500 for C-PPO, with $50\%$ error in $\tau_{i}$ ($i=x,z$),
$100\%$ in $K_{sur}$ and $50\%$ in $K_{i}$ ($i=1,2$). The histograms of
tracking errors are plotted in Fig. 5. C-PPO shows much better robustness to
model errors.

For this simple problem, PPO without stability guarantee can learn a stable
control policy with a network size $3\times 32$ after initial convergence at
Episode 4000. But the performance deteriorates as training progresses. With
stability guarantee, C-PPO maintains the stable performance. C-PPO also shows
improved robustness to model errors compared to a stable PPO policy. We also
observed that C-PPO can learn with a much larger network size (Fig. 6).

Figure 4: Stability of PPO vs PPO with stability guarantee (C-PPO). In the
stability test, we extended the controller running time (trajectory length)
from 8 s to 16 s to examine if PPO can perform stable control in the sense of
classic control or if PPO learns a complex curve fitting that cannot be
extended to longer controller running time. We used the setting with which PPO
has achieved its best performance in our hyper-parameters search without
stability guarantees. With this setting, PPO has achieved the same level of
tracking accuracy as C-PPO after initial convergence. With the stable learning
curves in (c) that show both PPO and C-PPO have learned the task reasonably
well, a seemingly well-learned PPO policy without a stability guarantee has
deteriorating performance, i.e., an increasing percentage of failed
trajectories in testing as training progresses, plotted in (a). Three types of
failures are observed in PPO: (b) and (d) plot the chattering motions; (e)
plots a trajectory drifting away; (f) plots a trajectory staying at a non-zero
equilibrium. On the other hand, C-PPO maintains the stable performance
regardless of training time. Figure 5: Robustness of PPO vs PPO with stability
guarantee (C-PPO). With $50\%$ model errors in $\tau_{i}$ ($i=x,z$), $100\%$
in $K_{sur}$ and $50\%$ in $K_{i}$ ($i=1,2$), the histograms of errors from
8000 randomly sampled trials of stable policies are plotted in (a) and (b).
The largest errors from C-PPO are under 0.1 (force) and 0.01 (position), while
PPO shows less robustness: the largest are around 2.5 and 0.25.  Figure 6:
Learning with different NN sizes: PPO vs PPO with stability guarantee (C-PPO).
With gradually increasing network sizes (the value and control networks) from
$3\times 32$ until $3\times 512$, PPO fails to learn at $3\times 256$ in (a),
while C-PPO shows some instabilities in learning at $3\times 512$ and when we
changed the value network to $3\times 32$ and kept the control networks to be
$3\times 512$, C-PPO learns steadily again in (b). Figure 7: Success rates of
HIRO vs HIRO with low-level stability guarantee (C-HIRO). The success rate is
computed from 50 random trials. In the ant maze test, C-HIRO reached a $\sim
60\%$ success rate in 10M steps while HIRO reached $\sim 10\%$. For the
following 5M steps, both stayed on those levels. In the contact maze test,
C-HIRO found the paths in 1M steps while HIRO took 2.5M steps.

###  V-B HIRO vs HIRO with stability guarantee (C-HIRO)

We implemented the control architecture in Fig. 2 for HIRO and applied Theorem
3 to learn two tasks: an ant maze and a contact maze (Fig. 1 (b) and (c)).

Since HIRO can learn the original ant-maze benchmark task with high sample
efficieny [3], we increased the difficulty of the ant-maze task by making a
less structured maze, starting with random initial positions in the maze to
reach randomly sampled goals, and adding one more joint for each leg. In the
contact maze task, the end-effector is supposed to reach a randomly sampled
goal in the green box without vision (Fig. 1 (c)). Details are in Appendix
VII-D and F.

In both example, C-HIRO shows better sample efficiency (Fig. 7). The reasons
of the improvement are three-pronged: (a) stability constraints decrease the
exploration space for low-level control; (b) the high-level action is in the
Cartesian space instead of the high-dimensional generalized coordinates; and
(c) low-level actions with stability guarantee may alleviate the non-
stationary issue for the high-level policy.

We also observed that the performance is very sensitive to the control
architecture in Fig. 2. For example, with different desired speeds, C-HIRO
reached around $30\%$ to $70\%$ success rates of the ant maze task in 10M
steps, although higher than HIRO. This may be related to the violations of the
assumptions in stability theorems or the interaction between the dynamic of
the low-level control system and the learning process. We will perform an
ablative study, apply contraction theory to analyze the learning process, and
further leverage control and path planning to learn multi-level tasks.

##  VI Concluding Remarks

By employing the control architecture in Fig. 1, contraction theory can be
applied for stability analysis of continuous RL control, yielding stability
constraints on network weights. The stability analysis requires a first-order
nonlinear model $\dot{\mathbf{y}}=\mathbf{f}(\mathbf{y},\mathbf{u},t)$
(specifically $\partial\mathbf{f}/\partial\mathbf{y}$ and
$\partial\mathbf{f}/\partial\mathbf{u}$) that describes the mapping from
$\mathbf{u}$ to $\mathbf{y}$, which may include different physical dynamics
and algorithms. Leveraging existing control techniques can create such a model
with respect to the end effector behaviors in the Cartesian space. This paper
demonstrates the importance and benefits of such stability guarantees for RL.
Current research focuses on the control system dynamics. In future, we will
expand our contraction analysis to the learning process.

##  VII Appendix

###  VII-A Contraction theory

The following is from contraction theory [1].

##### Definition of a contraction region

Given the system equations $\dot{\mathbf{x}}=\mathbf{f}(\mathbf{x},t)$, a
region of the state space is called a contraction region with respect to a
uniformly positive definite metric
$\mathbf{M}(\mathbf{x},t)=\mathbf{\Theta}^{T}\mathbf{\Theta}$, if the
generalized Jacobian
$\mathbf{F}=(\dot{\mathbf{\Theta}}+\mathbf{\Theta}\partial\mathbf{f}/\partial\mathbf{x})\mathbf{\Theta}^{-1}$
is uniformly negative definite in that region.

##### Theorem

Given the system equations $\dot{\mathbf{x}}=\mathbf{f}(\mathbf{x},t)$, any
trajectory, which starts in a ball of constant radius with respect to the
metric $\mathbf{M}(\mathbf{x},t)$, centered at a given trajectory and
contained at all times in a contraction region with respect to
$\mathbf{M}(\mathbf{x},t)$, remains in that ball and converges exponentially
to this trajectory. Furthermore global exponential convergence to the given
trajectory is guaranteed if the whole state space is a contraction region with
respect to the metric $\mathbf{M}(\mathbf{x},t)$.

###  VII-B Proof of Theorem 2

####  VII-B1 Contraction theory for hierarchical combinations

The following is from [1].

Consider the smooth differential relation of the form

$\frac{d}{dt}\begin{bmatrix}\delta\mathbf{z}_{1}\\\
\delta\mathbf{z}_{2}\end{bmatrix}=\begin{bmatrix}F_{11}&F_{12}\\\
0&F_{22}\end{bmatrix}\begin{bmatrix}\delta\mathbf{z}_{1}\\\
\delta\mathbf{z}_{2}\end{bmatrix}$ (28)

and assume bounded $F_{12}$. If $F_{11}$ and $F_{22}$ are uniformly negative
definite in some region, the system is contracting in that region [1].

In Eq. (28), the second equation does not depend on the first, so that
exponential convergence of $\delta\mathbf{z}_{2}$ to zero can be concluded for
uniformly negative definite $F_{22}$. Then $F_{12}\delta\mathbf{z}_{2}$ can be
seen as an exponentially decaying disturbance in the first equation.

Next we prove that a system with bounded disturbances is contracting in a
region if its generalized Jacobian (without considering the disturbance) is
uniformly negative definite in that region.

Given a system $\dot{\mathbf{x}}=\mathbf{f}(\mathbf{x},t)$, the differential
relation in a $M$ metric is written as

$\frac{d}{dt}\delta\mathbf{z}=F(\mathbf{x},t)\delta\mathbf{z}$ (29)

Consider the distance $R=\int_{P_{1}}^{P_{2}}\|\delta\mathbf{z}\|$ between two
trajectories $P_{1}$ and $P_{2}$, contained at all times in a contraction
region characterized by maximal eigenvalues
$\lambda_{max}(\mathbf{x},t)\leq-\beta<0$ of $\mathbf{F}$. The relative
velocity between these trajectories verifies $\dot{R}+|\lambda_{max}|R\leq 0$.
Assume now, instead, that $P_{1}$ represents a desired system trajectory and
$P_{2}$ the actual system trajectory in a disturbed flow field
$\dot{\mathbf{x}}=\mathbf{f}(\mathbf{x},t)+\mathbf{d}(\mathbf{x},t)$. Then

$\dot{R}+|\lambda_{max}|R\leq\|\Theta\mathbf{d}\|$ (30)

For bounded disturbance $\|\Theta\mathbf{d}\|$ any trajectory remains in a
boundary ball of Eq. (30) around the desired trajectory. Since initial
conditions $R(t=0)$ are exponentially forgotten, we can state that any
trajectory converges exponentially to a ball of radius $R$ in Eq. (30) with
arbitrary initial condition $R(t=0)$.

Since the first equation of Eq. (28),

$\delta\dot{\mathbf{z}}_{1}=F_{11}\delta\mathbf{z}_{1}+F_{12}\delta\mathbf{z}_{2}$
(31)

can be seen as a system of $\delta\mathbf{z}_{1}$ with exponentially decaying
disturbance $\delta\mathbf{z}_{2}$, a uniformly negative definite $F_{11}$
implies exponential convergence of $\delta z_{1}$ to an exponentially decaying
ball.

In summary, with bounded $F_{12}$ and uniformly negative $F_{ii}$ ($i=1,2.$),
the whole system globally exponentially converges to a single trajectory. By
recursion, this can be extended to systems similarly partitioned in more than
two equations.

Next, we prove Theorem 2 based on the above discussion of stability for
hierarchical combinations of equations.

####  VII-B2 The hierarchical combination

Here we show that with the proposed $W_{y}$, $W_{u}$, and $\pi_{z}$, Eq. (13)
can be reorganized into the form of a hierarchical combination Eq. (16).

Firstly, we show that
$W_{y}(\partial\mathbf{f}/\partial\mathbf{y})W_{y}^{-1}$,
$\partial\pi_{z}/\partial\mathbf{s}_{1}$, and
$\partial\pi_{z}/\partial\mathbf{s}_{2}$ are diagonal, and
$W_{y}(\partial\mathbf{f}/\partial\mathbf{u})W_{u}^{-1}$ is a upper triangular
matrix.

  * •

Diagonal $W_{y}(\partial\mathbf{f}/\partial\mathbf{y})W_{y}^{-1}$. For a
system Eq. (1) with diagonalizable $\partial\mathbf{f}/\partial\mathbf{y}$,
when $W_{y}$ is the normalized eigenvector matrix,
$W_{y}(\partial\mathbf{f}/\partial\mathbf{y})W_{y}^{-1}$ becomes the diagonal
eigenvalue matrix, denoted by $\Lambda$.

  * •

Diagonal $\partial\pi_{z}/\partial\mathbf{s}_{1}$ and
$\partial\pi_{z}/\partial\mathbf{s}_{2}$. The inner policy $\pi_{z}$ in Fig. 3
can be written as

$\pi_{z}(\mathbf{s}_{z})=\pi_{z}(\begin{bmatrix}\mathbf{s}_{1}\\\
\mathbf{s}_{2}\end{bmatrix})=\begin{bmatrix}h_{1}([s_{11},s_{21}])\\\\[8.0pt]
h_{2}([s_{12},s_{22}])\\\ \vdots\\\\[8.0pt]
h_{m}([s_{1m},s_{2m}])\end{bmatrix}$ (32)

The $ij^{th}$ components of $\partial\pi_{z}/\partial\mathbf{s}_{1}$ and
$\partial\pi_{z}/\partial\mathbf{s}_{2}$ are $\partial h_{i}/\partial s_{1j}$
and $\partial h_{i}/\partial s_{2j}$ respectively. When $i\neq j$, $\partial
h_{i}/\partial s_{1j}=0$ and $\partial h_{i}/\partial s_{2j}=0$.

  * •

Upper triangular $W_{y}(\partial\mathbf{f}/\partial\mathbf{u})W_{u}^{-1}$.
Denoting $W_{y}(\partial\mathbf{f}/\partial\mathbf{u})W_{u}^{-1}$ by $B$, here
we present two methods that makes $B$ a upper triangular matrix: QR
decomposition and LU decomposition.

    * –

QR decomposition:

$\displaystyle X=P(W_{y}\frac{\partial\mathbf{f}}{\partial\mathbf{u}})$ (33)
$\displaystyle X^{T}=QR$ $\displaystyle W_{u}=PQ^{T}$

where $P$ is a permutation matrix

$P=\begin{bmatrix}&&1\\\ &\iddots&\\\ 1&&\end{bmatrix}$ (34)

The resulting $B$ becomes

$\displaystyle B$
$\displaystyle=W_{y}\frac{\partial\mathbf{f}}{\partial\mathbf{u}}W_{u}^{-1}$
(35)
$\displaystyle=P(PW_{y}\frac{\partial\mathbf{f}}{\partial\mathbf{u}}(PQ^{T})^{-1}$
$\displaystyle=P(QR)^{T}(Q^{T})^{-1}P$ $\displaystyle=PR^{T}P$

which is an upper triangular matrix.

    * –

LU decomposition:

$\displaystyle
P(\tilde{W}_{e}\frac{\partial\mathbf{f}}{\partial\mathbf{u}})P=LU$ (36)
$\displaystyle W_{u}=PUP^{T}$

The resulting $B=PLP$.

Plugging $W_{y}(\partial\mathbf{f}/\partial\mathbf{y})W_{y}^{-1}=\Lambda$,
$W_{y}(\partial\mathbf{f}/\partial\mathbf{u})W_{u}^{-1}=B$, and
$\dot{W}_{y}=0$ into Eq. (13),

$\begin{bmatrix}\delta\dot{\mathbf{z}}\\\
\delta\dot{\mathbf{v}}\end{bmatrix}=\begin{bmatrix}\Lambda&B\\\
\frac{\partial\pi_{z}}{\partial\mathbf{s}_{1}}\Lambda+\frac{\partial\pi_{z}}{{\partial\mathbf{s}_{2}}}&\frac{\partial\pi_{z}}{\partial\mathbf{s}_{1}}B\end{bmatrix}\begin{bmatrix}\delta\mathbf{z}\\\
\delta\mathbf{v}\end{bmatrix}$ (37)

Denoting the components of $\delta\mathbf{z}$ and $\delta\mathbf{v}$ that
correspond to the $i^{th}$ eigenvalue of $\Lambda$ by $\delta z_{i}$ and
$\delta v_{i}$ respectively, the two equations with respect to
$\begin{bmatrix}\delta z_{i}&\delta v_{i}\end{bmatrix}$ do not contain any
terms of $\begin{bmatrix}\delta z_{j}&\delta v_{j}\end{bmatrix}$ where $j<i$.
That is, Eq. (37) can be reorganized into a hierarchical combination of
equations with respect to $\begin{bmatrix}\delta z_{i}&\delta
v_{i}\end{bmatrix}$. Here we show an example of two real eigenvalues, i.e.,
$i=1,2.$

$\displaystyle\begin{bmatrix}\delta\dot{z}_{1}&\delta\dot{v}_{1}&|\
\delta\dot{z}_{2}&\delta\dot{v}_{2}\end{bmatrix}^{T}=$ (38)
$\displaystyle\begin{bmatrix}\begin{array}[]{cc|cc}\lambda_{1}&b_{11}&0&b_{12}\\\\[8.0pt]
\frac{\partial h_{1}}{\partial s_{11}}\lambda_{1}+\frac{\partial
h_{1}}{\partial s_{21}}&\frac{\partial h_{1}}{\partial
s_{11}}b_{11}&0&\frac{\partial h_{1}}{\partial s_{11}}b_{12}\\\\[8.0pt]
\hline\cr 0&0&\lambda_{2}&b_{22}\\\\[8.0pt] 0&0&\frac{\partial h_{2}}{\partial
s_{12}}\lambda_{2}+\frac{\partial h_{2}}{\partial s_{22}}&\frac{\partial
h_{2}}{\partial s_{12}}b_{22}\\\\[8.0pt] \end{array}\end{bmatrix}$
$\displaystyle\begin{bmatrix}\delta z_{1}&\delta v_{1}&|\ \delta z_{2}&\delta
v_{2}\end{bmatrix}^{T}$

####  VII-B3 Stability

According to the theorem for hierarchical combinations, for a hierarchical
combination like Eq. (38), if in a region,

  * •

the diagonal blocks

$\begin{bmatrix}\lambda_{i}&b_{ii}\\\\[8.0pt] \frac{\partial h_{i}}{\partial
s_{1i}}\lambda_{i}+\frac{\partial h_{i}}{\partial s_{2i}}&\frac{\partial
h_{i}}{\partial s_{1i}}b_{ii}\\\\[8.0pt] \end{bmatrix}$ (39)

are uniformly negative definite,

  * •

the off-diagonal blocks are bounded,

the system is contracting in that region.

By solving the characteristic equations of the diagonal blocks

$s^{2}-(\lambda_{i}+\frac{\partial h_{i}}{\partial
s_{1i}}b_{ii})s+b_{ii}\frac{\partial h_{i}}{\partial s_{2i}}=0$ (40)

we obtain that if there exists $\epsilon>0$,

$\frac{\partial h_{i}}{\partial s_{1i}}b_{ii}+\lambda_{i}<-\epsilon,\ \ \
\frac{\partial h_{i}}{\partial s_{2i}}b_{ii}<-\epsilon$ (41)

the diagonal blocks are uniformly negative definite. For the second condition,
given bounded $\partial\mathbf{f}/\partial\mathbf{u}$, invertible $W_{y}$ and
$W_{u}$, and bounded $\partial h_{i}/\partial s_{1i}$, the off-diagonal blocks
are bounded. Proof is finished.

One can extend this for systems with complex eigenvalues. In that case, to
compute the range of the safe set may require numerical methods. Moreover, one
might extend it for repeated eigenvalues by using Jordan canonical form, which
may not be necessary for practice and in implementation is numerically
unstable.

###  VII-C Proof of Section III

Here we prove that with the implementation method in Section III, Theorem 2 is
satisfied.

Given $\lambda_{i}\leq 0$, the inequalities Eq. (17) in Theorem 2 can be
further simplified into

$\frac{\partial h_{i}}{\partial s_{1i}}b_{ii}<-\epsilon,\ \ \ \frac{\partial
h_{i}}{\partial s_{2i}}b_{ii}<-\epsilon$ (42)

because

$\frac{\partial h_{i}}{\partial
s_{1i}}b_{ii}<-\epsilon\leq-\epsilon-\lambda_{i}$ (43)

Thus if Eq. (42) is satisfied, the inequalities in Theorem 2 are satisfied.
Next we compute $\partial h_{i}/\partial s_{1i}$ and $\partial h_{i}/\partial
s_{2i}$.

The $h_{i}$ is a MLP with the hyperbolic tangent function as the activation
function. Assuming an $l$-layered MLP, for each layer,
$x_{j+1}=\tanh(W_{j}x_{j}+b_{j})$ where $j$ denotes the $j^{th}$ layer,
$W_{j}$ is the weights, and $b_{j}$ is the bias.

The partial derivatives $\partial h_{i}/\partial s_{1i}$ and $\partial
h_{i}/\partial s_{2i}$ are the products of weights and
$[1-\tanh^{2}{(\cdot)}]$’s. We further denote the weights of the first layer
corresponding to $s_{i1}$ and $s_{i2}$ respectively by $W_{1,1}$ and
$W_{1,2}$. Then $\partial h_{i}/\partial s_{1i}$ becomes

$\displaystyle\frac{\partial h_{i}}{\partial s_{i1}}=W_{l}M_{l-1}W_{l-1}\cdots
M_{2}W_{2}M_{1}W_{1,1}$ (44)

where $M_{j}=1-\tanh^{2}(W_{j}x_{j}+b_{j})$. When the activation function is
not saturated, components of $M_{j}$ are positive. Given positively bounded
weights by $\epsilon>0$ of the layers from $1$ to $l-1$, the components of the
product $(M_{l-1}W_{l-1}\cdots M_{2}W_{2}M_{1}W_{1,1})$ are also bounded by
some positive value, denoted by $\epsilon_{1}$.

Denoting the total number of nodes of the last layer by $k$, given

$\displaystyle\begin{cases}W_{lm}>-\epsilon b_{ii}/\lVert b_{ii}\rVert\ \ \
\text{when}\ \ b_{ii}<0\\\ W_{lm}<-\epsilon b_{ii}/\lVert b_{ii}\rVert\ \ \
\text{when}\ \ b_{ii}>0\\\ \end{cases}$ (45)
$\displaystyle(M_{l-1}W_{l-1}\cdots M_{2}W_{2}M_{1}W_{1,1})_{m}>\epsilon_{1}$

where the subscript $m$ denotes the $m^{th}$ component of the vector,

$\displaystyle\frac{\partial h_{i}}{\partial s_{i1}}b_{ii}$ (46)
$\displaystyle=$ $\displaystyle W_{l}(M_{l-1}W_{l-1}\cdots
M_{2}W_{2}M_{1}W_{1,1})b_{ii}$ $\displaystyle<$
$\displaystyle-k\epsilon_{1}\epsilon\frac{b_{ii}^{2}}{\lVert b_{ii}\rVert}$

Note that $b_{ii}=0$ means the system dynamics cannot be controlled by the RL
policy, which needs pre-designed control techniques for stability.

So there exists a small positive number $k\epsilon_{1}\epsilon
b_{ii}^{2}/\lVert b_{ii}\rVert$ that the first inequality in Eq. 42 is
satisfied. By replacing the $s_{1i}$ with $s_{2i}$, one can prove the second
inequality in Eq. (42).

###  VII-D Control architecture in Fig. 2

To make an end effector behave like a linear mass-spring-damper system, we use
Cartesian impedance control, e.g., [19], [21], and [22]. The resulting
dynamics can be represented by Eq. (18). The control techniques used in our
experiments are as follows.

####  VII-D1 Locomotion

Given rigid body dynamics with a floating base,

$H(\mathbf{q})\dot{\bm{\nu}}+C(\mathbf{q},\bm{\nu})+\mathbf{g}(\mathbf{q})=S^{T}\bm{\tau}_{c}+J_{c}^{T}(\mathbf{q})\bm{\gamma}_{c}$
(47)

where

$\mathbf{q}=\begin{pmatrix}\bm{\xi_{b}}\\\
\bm{\theta}\end{pmatrix},\bm{\nu}=\begin{bmatrix}\mathbf{v}_{b}\\\
\dot{\bm{\theta}}\end{bmatrix},S=\begin{bmatrix}\mathbf{0}_{6}&\\\
&\mathbf{I}_{n}\end{bmatrix}$ (48)

and $\bm{\xi}_{b}\in\mathrm{SE}(3)$ denotes the pose of the floating base with
respect to the inertial reference frame, $\bm{\theta}\in\mathbb{R}^{n}$
denotes the joints’ angles, $\mathbf{v}_{b}$ is the twist of the base, $H$ is
the generalized mass matrix, $C$ is the Coriolis and centrifugal terms, and
$\mathbf{g}$ is the gravitational term, $S$ is the selection matrix indicating
the actuated degrees of freedom, $\bm{\gamma}_{c}\in\mathbb{R}^{k}$ is the
external contact constraint force vector, and
$J_{c}(\mathbf{q})=\begin{bmatrix}J_{cb}&J_{cs}\end{bmatrix}$ is the contact
Jacobian with $J_{cb}\in\mathbb{R}^{k\times 6}$ and
$J_{cs}\in\mathbb{R}^{k\times n}$, we want the center of mass to behave like a
mass-spring-damper system.

In the ant maze task, we use the control method in [19] for the low-level
control, which includes (a) projected inverse dynamics, (b) impedance control
with the added RL action in the Cartesian space, and (c) quadratic programming
to maintain contact constraints. When the Cartesian location of the end
effector (denoted by $\mathbf{x}$) is subject to some external disturbance
$F_{x}$, the dynamics of the end effector can be written as

$\Lambda_{d}\ddot{\mathbf{e}}+K_{d}\dot{\mathbf{e}}+K_{p}\mathbf{e}+K_{rl}\mathbf{u}_{rl}=F_{x}$
(49)

where $\mathbf{e}=\mathbf{x}_{d}-\mathbf{x}$.

##### Projected inverse dynamics

The floating system dynamics Eq. (47) can be transformed into two equations of
two orthogonal spaces, motion space and constraint space via a projector
$P=I-J_{c}^{\dagger}J_{c}$. By defining
$\mathbf{h}=C(\mathbf{q},\bm{\nu})+\mathbf{g}(\mathbf{q})$, the equation for
the motion of the system can be written as

$PM\ddot{\mathbf{q}}+P\mathbf{h}=PB\bm{\tau}$ (50)

and the other equation for contact constraints

$(I-P)(M\ddot{\mathbf{q}}+\mathbf{h})=(I-P)B\mathbf{\tau}+J^{T}_{c}\bm{\gamma}_{c}$
(51)

Note that the torque term in Eq. (51), i.e., $(I-P)B\mathbf{\tau}$ doesn’t
contribute to the motion of the system (Eq. (50)). But the torque term in Eq.
(50) can affect the constraint equation (Eq. (51)) via $\ddot{\mathbf{q}}$.
Hence we firstly compute the torques

$\bm{\tau}_{m}=PB\bm{\tau}$ (52)

that control the motion behavior via Eq. (50) and then use quadratic
programming to compute the $\bm{\tau}$ to satisfy the contact constraints.

##### Cartesian impedance control

By defining the constraint inertia matrix $M_{c}=PM+I-P$ and denoting the
Jacobian of $\mathbf{x}$ we can rewrite Eq. (50) into

$M_{c}\ddot{\mathbf{q}}+P\mathbf{h}-\dot{P}\dot{\mathbf{q}}=PB\tau+PJ_{x}^{T}F_{x}$
(53)

We then pre-multiply Eq. (53) by $J_{x}M_{c}^{-1}$, replace
$J_{x}\ddot{\mathbf{q}}$ with $\ddot{\mathbf{x}}-\dot{J}_{x}\dot{\mathbf{q}}$,
and multiply by $\Lambda_{c}=(J_{x}M_{c}^{-1}PJ_{x}^{T})^{-1}$. The resulting
equation becomes

$\Lambda_{c}\ddot{\mathbf{x}}+\mathbf{h}_{c}=\Lambda_{c}J_{x}M_{c}^{-1}PB\bm{\tau}+F_{x}$
(54)

where
$\mathbf{h}_{c}=\Lambda_{c}J_{x}M_{c}^{-1}(P\mathbf{h}-\dot{P}\dot{\mathbf{q}})-\Lambda_{c}\dot{J}_{x}\dot{\mathbf{q}}$.

The torques Eq. (52) should satisfy the under-actuation constraint
$\bm{\tau}_{m}=B\bm{\tau}_{m}$. Thus we use
$\bm{\tau}_{m}=(PB)^{\dagger}J_{x}^{T}F$ where

$F=\Lambda_{c}\ddot{\mathbf{x}}_{d}+\Lambda_{c}\Lambda_{d}^{-1}[K_{d}\dot{\mathbf{e}}+K_{p}\mathbf{e}]+\mathbf{h}_{c}+\Lambda_{c}\Lambda_{d}^{-1}K_{rl}\mathbf{u}_{rl}$
(55)

The controller gains $\Lambda_{d}$, $K_{d}$ and $K_{p}$ define the preferred
behavior.

To satisfy contact constraints, we then use quadratic programming to minimize
$\|PB\bm{\tau}-\bm{\tau}_{m}\|_{2}^{2}$ under constraints.

We use the above controller for both the center of mass and the swing leg. As
for the center of the mass, the desired goal $\mathbf{x}_{d}$ is given by the
high-level policy, which is mapped to the equivalent velocity
$\mathbf{v}_{desired}=v_{speed}(\mathbf{x}_{d}-\mathbf{x})/\lVert\mathbf{x}_{d}-\mathbf{x}\rVert$
where $v_{speed}$ is a given desired speed. As for the swing leg, we followed
[20] to compute the desired foot contact location relative to the hip
abduction/adduction joints,

$\displaystyle\mathbf{r}_{desired}=\eta[0.5(\mathbf{v}_{hip}+\mathbf{v}_{com})-\mathbf{v}_{desired}]\sqrt{\frac{h}{g}}+0.5\mathbf{v}_{hip}\delta_{t}$
(56)

where $h$ is the current height of the hip with respect to the ground, $\eta$
is a scaling parameter, and $g$ is the gravitational acceleration. Note that
because in our experiments, the velocity of the center of mass is not constant
when the gaits are learned by a contact policy. So we didn’t include the
feedback term. Also we used a large desired velocity and a small $\eta=0.002$
because of large noises in the controller. Note that Eq. (56) may not be the
optimized contact locations for the ant dynamics but this is a reasonable
baseline to start with, given the similarity between the 3-joint 4-legged ant
and the quadrupedal robot.

##### Quadratic programming

From the constraint equation Eq. (51), we can solve for $\bm{\gamma}_{c}$,

$\bm{\gamma}_{c}=\bm{\eta}+\bm{\rho}$ (57)

where

$\displaystyle\bm{\eta}=-(J_{c}^{T})^{\dagger}(I-P)(I-MM_{c}^{-1}P)B\bm{\tau}$
(58)
$\displaystyle\bm{\rho}=(J_{c}^{T})^{\dagger}(I-P)[(I-MM_{c}^{-1}P)\mathbf{h}+MM_{c}^{-1}\dot{P}\dot{\mathbf{q}}]$

The contact constrains include (a) unilateral constraint to avoid loss of
contact $\gamma_{c,z}\geq 0$ and (b) friction cone constraint to avoid
slipping $\mu\gamma_{c,z}\geq\sqrt{\gamma_{c,x}^{2}+\gamma_{c,y}^{2}}$ where
$\mu$ is the friction coefficient at the contact point. Considering the
differences in the stance leg and swing leg, the QP becomes

$\displaystyle\min_{\bm{\tau}}$
$\displaystyle\frac{1}{2}\|PB\bm{\tau}-\bm{\tau}_{m}\|_{2}^{2}$ (59) subject
to $\displaystyle\eta_{z}^{i}+\rho_{z}^{i}\geq 0$
$\displaystyle\mu(\eta_{z}^{i}+\rho_{z}^{i})\geq\sqrt{(\eta_{x}^{i}+\rho_{x}^{i})^{2}+(\eta_{y}^{i}+\rho_{y}^{i})^{2}}$
$\displaystyle\bm{\tau}_{min}\leq\bm{\tau}\leq\bm{\tau}_{max}$

where the superscript $i$ denotes the $i$th stance leg and $\bm{\tau}_{min}$
and $\bm{\tau}_{max}$ are the torque saturations. Note that the above is
quadratically constrained quadratic programming. We approximate the friction
cone with a linearized 4-edge pryamid, which gives a linearly constrained
quadratic programming which can be solved with standard QP solvers.

##### About the contact policy

In the original HIRO paper [3], the high-level policy gives the desired
position of the center of mass and desired joint angles every $c$ steps. The
low-level policy produces the torque commands to track the high-level action.
Here we adapted HIRO as follows: (1) the high-level policy includes two
independent sub-policies: one provides the desired position ($\mathbf{x}_{d}$)
that gives the direction of a desired speed for the low-level policy and the
other contact policy provides the contact time. The contact policy is to
control the time of two contact modes (stance and swing) within the $c$ steps,
denoted by $\mathbf{p}\in[-1,1]$; when $\mathbf{p}>0$, the first $\mathbf{p}c$
steps are in the stance mode and from $\mathbf{p}c+1$ to $c$ in the swing
mode; when $\mathbf{p}<0$, the first $-\mathbf{p}c$ steps are in the swing
mode and the rest in stance.

####  VII-D2 Manipulation

Given

$H(\mathbf{q})\mathbf{\ddot{q}}+C(\mathbf{q},\mathbf{\dot{q}})\mathbf{\dot{q}}+\mathbf{g}(\mathbf{q})=\bm{\tau}+J^{T}(\mathbf{q})\bm{\gamma}$
(60)

where $\mathbf{q}\in\mathbb{R}^{n}$ is the vector of joint angles, $\bm{\tau}$
represents the torques applied on joints, $\bm{\gamma}$ represents the wrench
at the end effector, $H$ is the mass matrix, $C$ is the Coriolis and
centrifugal terms, and $\mathbf{g}$ is the gravitational term, we used force-
based control in [21] (without controlling the behavior in the null space).
The contact policy and the quadratic programming to satisfy the contact
constraints in Fig. 2 are not used. The Cartesian space control is

$\mathbf{\tau}=M\bar{J}[\ddot{\mathbf{x}}_{r}-\dot{J}\dot{\mathbf{q}}]+C(\mathbf{q},\dot{\mathbf{q}})\dot{\mathbf{q}}+\mathbf{g}$
(61)

where

$\displaystyle\ddot{\mathbf{x}}_{r}=\ddot{\mathbf{x}}_{d}+K_{d}(\dot{\mathbf{x}}_{d}-\mathbf{x})+K_{p}(\mathbf{x}_{d}-\mathbf{x})$
(62) $\displaystyle\bar{J}=M^{-1}J^{T}(JM^{-1}J^{T})^{-1}$

The high level action is defined the same as in locomotion, that is, the
absolute position in the Cartesian space, which is used as the direction of a
desired speed for low-level control.

Table I: Hyperparameters in C-HIRO vs HIRO experiments Parameter | Ant maze (C-HIRO) | Ant maze (HIRO) | Contact maze(C-HIRO) | Contact maze(HIRO)
---|---|---|---|---
$\gamma$ (Discount in value function) | HL: 0.99, LL:0.1 | HL: 0.99, LL:0.99 | HL: 0.99, LL: 0.99 | HL: 0.99, LL:0.99
$\sigma$ (STD in exploration) | HL:1.0, LL:1.0 | HL:1.0, LL:1.0 | HL:1.0, LL:0.2 | HL:1.0, LL:1.0
gradient clipping | 0.1 | 0.1 | 0.1 | 0.1
target update period | 200 | 200 | 100 | 100
train steps | HL: 400, LL:400 | HL: 400, LL:400 | HL:100, LL:200 | HL:100, LL:200

###  VII-E Proof of Theorem 3.

Here we prove that if the inequalities in Theorem 3 are satisfied, Theorem 2
are satisfied. Note that $z_{i}=s_{1i}$ and
$\sum_{0}^{k-1}z_{i}\delta_{t}=s_{2i}$ from Eq. (4).

Denoting the $i^{th}$ components of the diagonal positive definite matrices
$\tilde{\Lambda}_{y}$ and $\tilde{B}_{y}$ by $\tilde{\lambda}_{y,i}$ and
$\tilde{b}_{y,i}$ respectively, plugging

$\displaystyle\lambda_{i}=-\tilde{\lambda}_{y,i}\tilde{b}_{y,i}$ (63)
$\displaystyle b_{ii}=-\tilde{\lambda}_{y,i}K_{rl,i}$

into $(\partial h_{i}/\partial s_{1i})b_{ii}+\lambda_{i}$,

$\frac{\partial h_{i}}{\partial
s_{1i}}b_{ii}+\lambda_{i}=-\tilde{\lambda}_{y,i}(\frac{\partial
h_{i}}{\partial
s_{1i}}K_{rl,i}+\tilde{b}_{y,i})<-\tilde{\lambda}_{y,i}\epsilon$ (64)

because $\tilde{b}_{y,i}>0$, $\tilde{\lambda}_{y,i}>0$, and $(\partial
h_{i}/\partial s_{1i})K_{rl,i}>\epsilon$. Hence there exists
$\epsilon_{1}=\tilde{\lambda}_{y,i}\epsilon>0$ that the first inequality in
Theorem 2 is satisfied.

When $(\partial h_{i}/\partial s_{2i})K_{rl,i}>\epsilon$,

$\displaystyle\frac{\partial h_{i}}{\partial s_{2i}}b_{ii}=$
$\displaystyle-\frac{\partial h_{i}}{\partial
s_{2i}}\tilde{\lambda}_{y,i}K_{rl,i}<$
$\displaystyle-\tilde{\lambda}_{y,i}\epsilon$ (65)

and hence the second inequality in Theorem 2 is also satisfied.

Because $W_{y}$ is identity, the assumption in Theorem 2 that $\dot{W}_{y}$
can be ignored is also satisfied. According to Theorem 2, the system is
contracting. Proof is finished.

###  VII-F Environments, Hyperparameters, and Discussions

####  VII-F1 2D Math Example

##### RL hyperparameters

We use the following hyperparameters:

  * •

entropy coefficient: 0.0

  * •

learning rate: 0.0001,

  * •

vf coefficient: 1.0

  * •

max gradient norm: 0.5

  * •

gamma:0.99

  * •

lambda: 0.95

  * •

No. of minibatches (nminibatches): 4

  * •

clip range:0.1

  * •

No. of optimized epoches (noptepochs): 10

As for the parameters not listed, we used the default parameters from the
OpenAI baselines package.

##### Environments

The coefficients of the system dynamics Eq. (27) are $\tau_{z}=0.0437$,
$\tau_{x}=0.01$, $K_{1}\in[-0.01,0.01]$, $K_{2}\in[-0.01,0.01]$, $x\in[1,3]$,
$z\in[-3,1]$, and $K_{sur}\in[1,31]$.

We tested different RL state space designs: (a)
$\begin{bmatrix}x&f\end{bmatrix}$, (b)
$\begin{bmatrix}e_{x}&e_{f}\end{bmatrix}$, (c)
$\begin{bmatrix}e_{x}&e_{f}&\sum e_{x}&\sum e_{f}\end{bmatrix}$ without NN
weight constraints, and (d) $\begin{bmatrix}e_{x}&e_{f}&\sum e_{x}&\sum
e_{f}\end{bmatrix}$ with NN weight constraints. The first two failed to reach
100% success rate (error$<0.4$). The above results are from (c) and (d).

##### Theorem 4

To avoid non-zero error equilibriums in Fig. 4 (f), we develop Theorem 4 as
follows: Given the nonlinear dynamics
$\dot{\mathbf{e}}=\mathbf{f}(\mathbf{e},\mathbf{u},t)$, we use
$\mathbf{s}_{1}=\mathbf{e}$ and $\mathbf{s}_{2}=\int_{0}^{t}\mathbf{e}dt$. The
nonlinear dynamics with a time-invariant $\pi$ can be approximated by its
locally linearized model,

$\displaystyle\dot{\mathbf{e}}=\frac{\partial\mathbf{f}}{\partial\mathbf{e}}\mathbf{e}+\frac{\partial\mathbf{f}}{\partial\mathbf{u}}\mathbf{u}+\frac{\partial\mathbf{f}}{\partial
t}t$ (66)
$\displaystyle\dot{\mathbf{u}}=\frac{\partial\pi}{\partial\mathbf{s}_{1}}\dot{\mathbf{e}}+\frac{\partial\pi}{\partial\mathbf{s}_{2}}\mathbf{e}$

Assuming $\frac{\partial\mathbf{f}}{\partial t}t$ converges to
$\mathbf{f}_{stable}$ as $t$ approaches infinity, the steady states
($\dot{\mathbf{e}}=0$ and $\dot{\mathbf{u}}=0$) are

$\begin{cases}\mathbf{e}=\mathbf{0}\\\\[8.0pt]
\frac{\partial\mathbf{f}}{\partial\mathbf{u}}\mathbf{u}+\mathbf{f}_{stable}=\mathbf{0}\end{cases}\text{and}\
\
\begin{cases}\frac{\partial\pi}{\partial\mathbf{s}_{2}}=\mathbf{0}\\\\[8.0pt]
\frac{\partial\mathbf{f}}{\partial\mathbf{e}}\mathbf{e}+\frac{\partial\mathbf{f}}{\partial\mathbf{u}}\mathbf{u}+\mathbf{f}_{stable}=0\end{cases}$
(67)

Accordingly one can limit the static equilibrium of the system to be only at
the desired value, that is, Theorem 4.

Theorem 4 If in the entire working space,

  * •

there exists $\mathbf{u}$ in the control action space that

$\frac{\partial\mathbf{f}}{\partial\mathbf{u}}\mathbf{u}+\mathbf{f}_{stable}=\mathbf{0}$
(68)

  * •

for everywhere in the action space,

$\frac{\partial\pi_{z}}{\partial\mathbf{s}_{2}}\neq\mathbf{0}$ (69)

Eq. (1) has one and only one static equilibrium at the desired value, i.e.,
$\mathbf{e}=0$.

####  VII-F2 Ant maze and contact maze

##### RL hyper-parameters

We use the same set of parameters as in the HIRO paper except those in Table.
I:

In the contracting policy, we use a $3\times 400$ MLP for each axis. Other
networks are the same as the HIRO paper.

##### Environments

The ant maze is from $x\in(-1,12)$ and $y\in(-1,14)$ and the ant is randomly
put inside the maze initially. The target is randomly sampled inside
$x\in(1,8)$ and $y\in(2,11)$, the right-side area with more blocks. The
success is defined as the L2 norm of distance is smaller than 2. The ant has
3-joint legs. The first two are at the hip within range of 60 degrees (z
direction) and 10 degrees in the perpendicular plane at the neural hip
position. The third joint is the kneel at the leg within range of 40 degrees
in the same perpendicular plane at the neural hip position. In the contact
maze, the target is sampled between $x\in(0.41,0.415)$, $y\in(-0.025,-0.02)$,
and $z\in(0.005,0.015)$ (centered at the green box) with the success defined
as the L2 norm is less than 0.005. The 7-DOF robot used is the Franka panda
robot.

## Acknowledgments

This study was supported in part under the RIE2020 Industry Alignment Fund –
Industry Collaboration Projects (IAF-ICP) Funding Initiative, as well as cash
and in-kind contribution from the industry partner, HP Inc., through the HP-
NTU Digital Manufacturing Corporate Lab.

## References

  * [1] W. Lohmiller and J.-J. E. Slotine, “On contraction analysis for non-linear systems,” _Automatica_ , vol. 34, no. 6, pp. 683–696, 1998.
  * [2] J. Schulman, F. Wolski, P. Dhariwal, A. Radford, and O. Klimov, “Proximal policy optimization algorithms,” _arXiv preprint arXiv:1707.06347_ , 2017\.
  * [3] O. Nachum, S. S. Gu, H. Lee, and S. Levine, “Data-efficient hierarchical reinforcement learning,” _Advances in neural information processing systems_ , vol. 31, 2018.
  * [4] F. Berkenkamp, M. Turchetta, A. P. Schoellig, and A. Krause, “Safe model-based reinforcement learning with stability guarantees,” _arXiv preprint arXiv:1705.08551_ , 2017.
  * [5] M. Jin and J. Lavaei, “Stability-certified reinforcement learning: A control-theoretic perspective,” _IEEE Access_ , vol. 8, pp. 229 086–229 100, 2020.
  * [6] S. A. Khader, H. Yin, P. Falco, and D. Kragic, “Stability-guaranteed reinforcement learning for contact-rich manipulation,” _IEEE Robotics and Automation Letters_ , vol. 6, no. 1, pp. 1–8, 2020.
  * [7] N. M. Boffi and J.-J. E. Slotine, “Implicit Regularization and Momentum Algorithms in Nonlinearly Parameterized Adaptive Control and Prediction,” _Neural Computation_ , vol. 33, no. 3, pp. 590–673, 03 2021. [Online]. Available: https://doi.org/10.1162/neco_a_01360
  * [8] S. Singh, S. M. Richards, V. Sindhwani, J.-J. E. Slotine, and M. Pavone, “Learning stabilizable nonlinear dynamics with contraction-based regularization,” _The International Journal of Robotics Research_ , vol. 40, no. 10-11, pp. 1123–1150, 2021.
  * [9] I. R. Manchester, M. Revay, and R. Wang, “Contraction-based methods for stable identification and robust machine learning: a tutorial,” _arXiv preprint arXiv:2110.00207_ , 2021.
  * [10] D. Sun, S. Jha, and C. Fan, “Learning certified control using contraction metric,” _arXiv preprint arXiv:2011.12569_ , 2020.
  * [11] H. Tsukamoto, S.-J. Chung, and J.-J. E. Slotine, “Contraction theory for nonlinear stability analysis and learning-based control: A tutorial overview,” _Annual Reviews in Control_ , 2021.
  * [12] H. Tsukamoto, S.-J. Chung, J.-J. Slotine, and C. Fan, “A theoretical overview of neural contraction metrics for learning-based control with guaranteed stability,” _arXiv preprint arXiv:2110.00693_ , 2021.
  * [13] I. R. Manchester and J.-J. E. Slotine, “Control contraction metrics: Convex and intrinsic criteria for nonlinear feedback design,” _IEEE Transactions on Automatic Control_ , vol. 62, no. 6, pp. 3046–3053, 2017.
  * [14] ——, “Robust control contraction metrics: A convex approach to nonlinear state-feedback $H^{\infty}$ control,” _IEEE Control Systems Letters_ , vol. 2, no. 3, pp. 333–338, 2018.
  * [15] ——, “Control contraction metrics and universal stabilizability,” _IFAC Proceedings Volumes_ , vol. 47, no. 3, pp. 8223–8228, 2014.
  * [16] H. Tsukamoto, S.-J. Chung, and J.-J. E. Slotine, “Neural stochastic contraction metrics for learning-based control and estimation,” _IEEE Control Systems Letters_ , vol. 5, no. 5, pp. 1825–1830, 2020.
  * [17] P. Zhao, A. Lakshmanan, K. Ackerman, A. Gahlawat, M. Pavone, and N. Hovakimyan, “Tube-certified trajectory tracking for nonlinear systems with robust control contraction metrics,” _arXiv preprint arXiv:2109.04453_ , 2021.
  * [18] B. T. Lopez, J.-J. E. Slotine, and J. P. How, “Dynamic tube mpc for nonlinear systems,” in _2019 American Control Conference (ACC)_ , 2019, pp. 1655–1662.
  * [19] G. Xin, W. Wolfslag, H.-C. Lin, C. Tiseo, and M. Mistry, “An optimization-based locomotion controller for quadruped robots leveraging cartesian impedance control,” _Frontiers in Robotics and AI_ , vol. 7, p. 48, 2020.
  * [20] 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 international conference on Robotics and automation_. IEEE, 2013, pp. 3287–3292.
  * [21] J. Nakanishi, R. Cory, M. Mistry, J. Peters, and S. Schaal, “Operational space control: A theoretical and empirical comparison,” _The International Journal of Robotics Research_ , vol. 27, no. 6, pp. 737–757, 2008.
  * [22] B. Siciliano and L. Villani, _Robot force control_. Springer Science & Business Media, 1999.