# Stable Modular Control via Contraction Theory
for Reinforcement Learning

Bing Song  HP-NTU Digital
Manufacturing Corporate Lab
Singapore, Singapore 639798
<EMAIL_ADDRESS>Jean-Jacques Slotine  Nonlinear Systems Laboratory
Massachusetts Institute of Technology
Cambridge MA 02139-4307
<EMAIL_ADDRESS>Quang-Cuong Pham  HP-NTU Digital
Manufacturing Corporate Lab
Singapore, Singapore 639798
<EMAIL_ADDRESS>
###### Abstract

We propose a novel way to integrate control techniques with reinforcement
learning (RL) for stability, robustness, and generalization: leveraging
contraction theory to realize modularity in neural control, which ensures that
combining stable subsystems can automatically preserve the stability. We
realize such modularity via signal composition and dynamic decomposition.
Signal composition creates the latent space, within which RL applies to
maximizing rewards. Dynamic decomposition is realized by coordinate
transformation that creates an auxiliary space, within which the latent
signals are coupled in the way that their combination can preserve stability
provided each signal, that is, each subsystem, has stable self-feedbacks.
Leveraging modularity, the nonlinear stability problem is deconstructed into
algebraically solvable ones, the stability of the subsystems in the auxiliary
space, yielding linear constraints on the input gradients of control networks
that can be as simple as switching the signs of network weights. This
minimally invasive method for stability allows arguably easy integration into
the modular neural architectures in machine learning, like hierarchical RL,
and improves their performance. We demonstrate in simulation the necessity and
the effectiveness of our method: the necessity for robustness and
generalization, and the effectiveness in improving hierarchical RL for
manipulation learning.

###### Index Terms:

modularity, contraction theory, stability, reinforcement learning

##  I Introduction

Reinforcement learning (RL) has been struggling with control stability,
robustness, and generalization. How to leverage control theoretical results to
solve those issues is one fundamental question. Currently, control theoretical
results are employed in RL to create some priors[1], to design loss
functions[2], or to decrease the gap between the simulated environments and
the real-world applications by improving dynamic models[3], producing
adversarial samples[4], etc.

Here we propose a novel way to integrate control with RL: applying contraction
theory[5] to creating modular architectures in neural control. Leveraging
modularity, we develop a minimally invasive way to realize stable control with
RL.

Figure 1: Modular control framework for RL. The modularity, for example,
hierarchical combination of subsystems, is realized in the auxiliary space by
coordinate transformation ($T_{y},T_{a}$). Stability can be automatically
guaranteed provided each subsystem has stable self-feedbacks, yielding linear
constraints on the input gradients of control networks that can be as simple
as switching the signs of weights.

To guarantee stability is more challenging for RL than for nonlinear control.
RL targets at control problems likely involving large variations,
uncertainties, unknowns, and non-convexities. To handle those difficulties, RL
characterizes physical systems with distributions (rigorously, transition
functions) and optimizes neural control according to estimated values (e.g.,
value functions and gradients). However, the stability of control is
mathematically defined in nonlinear system theory for the underlying dynamics
behind those distributions and estimations.

Complexities arise from this mismatch in formulation: stabilizing control of
physical dynamics based on models vs. optimizing control for task
distributions based on data.

On one hand, it raises unique concerns. For example, in [6], we show that
local optimal policies can contain unstable control for some dynamic
parameters in the task distribution, even when there exist better policies
with higher expected returns that provide stable control for all parameters.
Overfitting to such instabilities can deteriorate generalization.

On the other hand, it creates intrinsically difficult questions to extend
control-theoretical results to RL. For example, the wide success of neural
certificates in nonlinear control[7], both empirically and theoretically, has
not yet been equally transferred to RL. The success of neural certificates
depends on that as more stable trajectories are collected to learn the
certificates, the bounds on the generalization error of the data-driven
certificates to new trajectories are shrinking[8]. The controllers are derived
or optimized based on the learned certificates to perform stable control.
Extending neural certificates to RL, the controllers need to be adapted
simultaneously to maximize the estimated value functions for possibly non-
convex problems and for possibly much larger variations in dynamic parameters.
It raises mathematically challenging questions like whether the generalization
error bounds for the certificates are still shrinking and whether the
exploration for value estimation is unnecessarily limited. Currently in safe
RL, there exists intuitive methods like weighted sum of objectives [9] and
alternative updates[10], and complicated methods like unifying the objective
functions into a multi-timescale min-max-min optimization problem[11]. All
lack the analysis on the interference between the generalization bounds of the
certificates and the estimation of value functions. Expertise in both
nonlinear control and RL is likely required for the empirical certificate loss
design as well as task-specific hyperparameter tuning.

In contrast, learning motor skills with stable control seems easy and natural
to animals. It is believed in both machine learning[12] and nonlinear
control[13] that modularity is the key. Here we introduce into RL the
modularity from nonlinear control, in particular, from contraction theory[5].

More than functional specialization, modularity in biology also ensures that
combining stable subsystems can automatically preserve the stability, which
can be mathematically characterized by contraction theory[14].

We leverage such modularity to deconstruct the nonlinear stability problem
into algebraically solvable ones, the stability of elementary subsystems,
yielding linear constraints on the input gradients of control networks
(network Jacobians).

Such modularity is realized by signal composition and dynamic decomposition.
Signal composition creates the latent space, within which RL applies to
maximizing rewards. Dynamic decomposition is realized by coordinate
transformation that creates an auxiliary space, within which the latent
signals are coupled in the way that their combination can preserve stability
as long as each signal has stable self-feedbacks, for example, the
hierarchical type in Fig. 1. The graph of the signal flows across different
spaces is illustrated in Fig. 2.

The modular architecture simplifies the nonlinear stability problem. Provided
stable elementary subsystems in the auxiliary space, the latent signals in the
latent space are stable, converging with bounded transients. Provided that the
composition of latent signals is also a contracting system, the convergence of
latent signals implies the stability of physical dynamics.

In this paper, we use control techniques to create latent signals, in
particular, composite variables that mix errors in position and velocity,
which is common in biological nervous systems[13]. This framework may be
extended to latent signals from machine learning.

To decompose dynamics in the auxiliary space, the dynamic model in the latent
space is required, specifically, the partial derivatives. For unknown
environments, one may apply system identification, learning to predict those
partial derivatives. Here we provide another solution, leveraging task space
controllers[15] to build a latent-space model composed of the controller
gains. Robot models are still required.

Figure 2: Signal flows across different spaces. Signal composition creates the
latent space, within which neural control with RL applies to maximizing
rewards. Provided that the composition of latent signals is a contracting
system, the stability of latent signals implies the stability of physical
dynamics. Dynamic decomposition in the auxiliary space further simplifies the
stability of latent signals into algebraically solvable ones via modularity.

Robustness of our method to model errors and bounded disturbances can be
characterized via contraction theory.

Our method provides a minimally invasive way for stable RL control. The
pseudocode of RL algorithms remains the same, while the graph of control
policies should be adapted into the modular control framework with coordinate
transformation and the constrained control networks, as well as adding
function compositions that define the RL states (latent signals) and
integrating RL actions with existing controllers. The stability constraints
can be as simple as switching the signs of network weights for multilayer
perceptrons (MLPs).

Following the implementation and limitations, we demonstrate the necessity and
the effectiveness of our method in simulation. The necessity for robustness
and generalization is illustrated in the complementary paper [6] and we
briefly include the results here. The effectiveness is demonstrated with
hierarchical RL (two TD3[16] agents, one for planning and the other for
control), learning two manipulation scenarios: peg push and peg maze. The
stability constraints on low-level control guarantee consistent dynamic
behaviors against perturbations and variations in initial distributions,
mitigating concerns like non-stationary distributions. Results shows that our
method outperforms (i) the state-of-the-art data-efficient hierarchical
RL[17], (ii) hierarchical RL with unconstrained modular control, and (iii)
one-level RL with stable non-neural controllers.

###  I-A Technical Contributions.

Our main contribution is the method that builds an auxiliary space to
decompose the dynamics for modular neural control, deconstructing the
nonlinear stability problem into algebraically solvable ones with linear
control techniques. This lowers the barrier to extend control-theoretical
results to RL. The minimally invasive way allows arguably easy integration
into hierarchical RL, and improves its performance by providing consistent
low-level dynamic behaviors. This framework also opens opportunities to
further integrate with machine learning, in particular, integration with
learned latent signals.

Other contributions include (i) illustrating how to transfer the stability of
latent signals to physical dynamics with the example of composite variables,
and (ii) leveraging task space control to create the latent-space dynamic
model composed of the controller gains.

###  I-B Related work.

####  I-B1 Contraction theory and safe RL

Compared with Lyapunov theory, contraction theory analyzes differential
dynamics without specifying the system equilibrium and hence robust to the
perturbations that shift the system equilibrium [18], which are common in the
real world, particularly for contact-rich manipulation learning.

Contraction theory has been introduced into learning for robust control [19],
adaptive control [20], motion planning [21], and system identification [22].

In safe RL, Krasovskii-Constrained Reinforcement Learning [23] based on
contraction analysis solves the stability constraints with a primal dual
policy gradient approach. Another work [24] learns a control contraction
metric and a tracking controller by optimizing an empirical risk function.

Our novelty lies in leveraging the modularity of contraction to deconstruct
the problem into combinations of algebraically solvable ones. To our best
knowledge, this is the first time that the modularity of contraction is
introduced into RL.

Our theoretical results on the constraints for incremental exponential
stability are in line with [25] based on integral quadratic constraints for
$\mathcal{L}_{2}$ stability, which limits the input gradients of control
networks (network Jacobians) for stable linear systems with added
nonlinearity. Here we derive the constraints for general nonlinear systems by
creating the modularity to deconstruct the nonlinear problem.

Other related work in safe RL, robustness, and generalization can be found in
[6], the complementary paper that illustrates the necessity of stability
constraints for robustness and generalization.

####  I-B2 Modularity

In control, planning, and learning, several threads leverage bioinspired
modularity to simplify intrinsically difficult or computationally complex
problems.

The well-known forms of modularity are motor primitives[26] and synergies[27].
Our framework can be possibly generalized to involve that form of modularity,
for example, by expanding the fully actuated 2D formulation of subsystems to
underactuated multi-dimensional subsystems.

Another important direction is modular motion generation, e.g., Riemannian
motion policies[28] pairing different Riemannian metrics to the subsets of the
state space defined by different mappings from states to the desired
acceleration (the polices for motion generation). Our work shares the spirit
in the way that we create new Riemannian metrics (made of the coordinate
transformation) to define the auxiliary space for dynamic deconstruction.

The modularity in machine learning[29] is studied for different purposes from
the perspective of information theory and computation, formulating the
modularity as function compositions and the like. From the perspective of
nonlinear system theory, the modularity of contraction theory solves the
stability of function compositions with their differential dynamics. We apply
those results and analyze stability with differential dynamics, the linear
form of which (network Jacobians are embedded in the coefficient matrices)
highly simplifies the analysis.

###  I-C Paper organization

  1. 1.

Dynamic decomposition.

  2. 2.

Signal composition.

  3. 3.

For unknown environments.

  4. 4.

Implementation and limitations.

  5. 5.

Example 1: necessity for robustness and generalization.

  6. 6.

Example 2: manipulation learning.

  7. 7.

Concluding remarks.

  8. 8.

Appendix

##  II Dynamic decomposition

Contraction theory analyzes the differential dynamics to guarantee
convergence[5]. Intuitively, provided that the distances between arbitrary
trajectories are converging to zero, all possible trajectories are converging
to the system limit behaviors, e.g., equilibria or limit cycles.

Provided $\mathbf{\dot{y}=\mathbf{f}(\mathbf{y},\mathbf{u})}$ in the latent
space, this section includes (i) contraction from auxiliary space to latent
space, (ii) differential dynamics in the auxiliary space, (iii) modular neural
control via dynamic decomposition, and (iv) explicit solutions to the
hierarchical-type modularity.

###  II-A Contraction from auxiliary space to latent space

Figure 3: By constraining the dynamics in the auxiliary space to converge
exponentially, the dynamics in the Euclidean space is converging with bounded
transients.

How the stability is transferred from auxiliary space ($\mathbf{z}$ and
$\mathbf{a}$) to latent space ($\mathbf{y}$ and $\mathbf{u}$) is illustrated
in Fig. 3.

Briefly speaking, provided contraction in the auxiliary space, the
transformation $T_{y}$ and $T_{a}$ maps the contraction metric $M$ in the
auxiliary space to the metric $\Theta_{M}^{T}M\Theta_{M}$ in the original
latent space, where $\Theta_{M}$ consists of $T_{y}$ and $T_{a}$. Provided
that $\Theta_{M}^{T}M\Theta_{M}$ is uniformly positive definite, the dynamics
in the original space is also converging with bounded transients. The
stability theorem and its proof is in Appendix. X-B.

###  II-B Differential dynamics in the auxiliary space

Here we derive the differential dynamics in the auxiliary space. Note that in
the following discussion, we include the integral of $\mathbf{z}$ in
$\mathbf{s}$ to force the system equilibrium at the goal (Section IV in [6]).
Proofs of dynamic convergence in this paper hold without the integral.

Given the neural-control system in the original space, by adding two interface
layers $T_{y}$ and $T_{a}$, we create the new coordinates

$\begin{bmatrix}\mathbf{z}\\\
\mathbf{a}\end{bmatrix}=\Theta_{M}\begin{bmatrix}\mathbf{y}\\\
\mathbf{u}\end{bmatrix}\ \ \text{where}\ \ \Theta_{M}=\begin{bmatrix}T_{y}&\\\
&T_{a}^{-1}\end{bmatrix}$ (1)

yielding the differential relation

$\begin{bmatrix}\delta\dot{\mathbf{z}}\\\
\delta\dot{\mathbf{a}}\end{bmatrix}=\dot{\Theta}_{M}\begin{bmatrix}\delta\mathbf{y}\\\
\delta\mathbf{u}\end{bmatrix}+\Theta_{M}\begin{bmatrix}\delta\dot{\mathbf{y}}\\\
\delta\dot{\mathbf{u}}\end{bmatrix}$ (2)

where $\delta\mathbf{y}$, $\delta\mathbf{u}$, $\delta\mathbf{z}$, and
$\delta\mathbf{a}$ are differential displacements, i.e., infinitesimal
displacements at fixed time. The control policy $\bm{\pi}$ can be written as

$\mathbf{a}=\bm{\pi}(\begin{bmatrix}\mathbf{s}_{1}\\\
\mathbf{s}_{2}\end{bmatrix})=\bm{\pi}(\begin{bmatrix}\mathbf{z}\\\
\int^{t}_{0}\mathbf{z}dt\end{bmatrix})$ (3)

which yields

$\delta\dot{\mathbf{a}}=\frac{\partial\bm{\pi}}{\partial\mathbf{s}_{1}}\delta\dot{\mathbf{z}}+\frac{\partial\bm{\pi}}{{\partial\mathbf{s}_{2}}}\delta\mathbf{z}$
(4)

given
$\delta\dot{\mathbf{s}}_{2}=\delta(\frac{d}{dt}\int^{t}_{0}\mathbf{z}dt)=\delta\mathbf{z}$
and the small change in the input gradient of $\bm{\pi}$, i.e.,
$\delta\frac{\partial\bm{\pi}}{\partial\mathbf{s}}$, can be ignored. The
dynamics $\dot{\mathbf{y}}=\mathbf{f}(\mathbf{y},\mathbf{u})$ yields the
differential relation

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

Combining Eq. (2), Eq. (4), and Eq. (5), we obtain the differential dynamics
in the auxiliary space,

$\begin{bmatrix}\delta\dot{\mathbf{z}}\\\
\delta\dot{\mathbf{a}}\end{bmatrix}=\begin{bmatrix}A&B\\\ C&D\\\
\end{bmatrix}\begin{bmatrix}\delta\mathbf{z}\\\ \delta\mathbf{a}\end{bmatrix}$
(6)

where

$\displaystyle A$
$\displaystyle=[\dot{T}_{y}+T_{y}\frac{\partial\mathbf{f}}{\partial\mathbf{y}}]T_{y}^{-1}$
(7) $\displaystyle B$
$\displaystyle=T_{y}\frac{\partial\mathbf{f}}{\partial\mathbf{u}}T_{a}$
$\displaystyle C$
$\displaystyle=\frac{\partial\bm{\pi}}{\partial\mathbf{s}_{1}}(\dot{T}_{y}+T_{y}\frac{\partial\mathbf{f}}{\partial\mathbf{y}})T_{y}^{-1}+\frac{\partial\bm{\pi}}{{\partial\mathbf{s}_{2}}}$
$\displaystyle D$
$\displaystyle=\frac{\partial\bm{\pi}}{\partial\mathbf{s}_{1}}T_{y}\frac{\partial\mathbf{f}}{\partial\mathbf{u}}T_{a}$

Note that $\mathbf{s}_{1}=\mathbf{z}$ and $\mathbf{s}_{2}=\int\mathbf{z}dt$.

###  II-C Modular neural control via dynamic decomposition

To illustrate the design of modular neural control, we use
$\mathbf{z}\in\mathbf{R}^{2}$ for simplicity, which can be generalized to
higher dimensions. Given $\mathbf{z}\in\mathbf{R}^{2}$, we rearrange Eq. (6)
by grouping $z_{i}$ and $a_{i}$, that is,

$\begin{bmatrix}\delta\dot{z}_{1}\\\ \delta\dot{a}_{1}\\\
\hline\cr\delta\dot{z}_{2}\\\
\delta\dot{a}_{2}\end{bmatrix}=\begin{bmatrix}a_{11}&b_{11}&\\!\\!\\!\framebox{}\\!\\!\\!&a_{12}&b_{12}\\\
c_{11}&d_{11}&\\!\\!\\!\framebox{}\\!\\!\\!&c_{12}&d_{12}\\\ \hline\cr
a_{21}&b_{21}&\\!\\!\\!\framebox{}\\!\\!\\!&a_{22}&b_{22}\\\
c_{21}&d_{21}&\\!\\!\\!\framebox{}\\!\\!\\!&c_{22}&d_{22}\\\
\end{bmatrix}\begin{bmatrix}\delta z_{1}\\\ \delta a_{1}\\\ \hline\cr\delta
z_{2}\\\ \delta a_{2}\end{bmatrix}$ (8)

where $a_{ij}$, $b_{ij}$, $c_{ij}$ and $d_{ij}$ denote the elements of $A$,
$B$, $C$, and $D$ respectively.

The system is then deconstructed into the combinations of two subsystems:

$\displaystyle\begin{bmatrix}\delta\dot{z}_{1}\\\
\delta\dot{a}_{1}\end{bmatrix}=\begin{bmatrix}a_{11}&b_{11}\\\
c_{11}&d_{11}\\\ \end{bmatrix}\begin{bmatrix}\delta z_{1}\\\ \delta
a_{1}\end{bmatrix}+\begin{bmatrix}a_{12}&b_{12}\\\ c_{12}&d_{12}\\\
\end{bmatrix}\begin{bmatrix}\delta z_{2}\\\ \delta a_{2}\end{bmatrix}$ (9)
$\displaystyle\begin{bmatrix}\delta\dot{z}_{2}\\\
\delta\dot{a}_{2}\end{bmatrix}=\begin{bmatrix}a_{22}&b_{22}\\\
c_{22}&d_{22}\\\ \end{bmatrix}\begin{bmatrix}\delta z_{2}\\\ \delta
a_{2}\end{bmatrix}+\begin{bmatrix}a_{21}&b_{21}\\\ c_{21}&d_{21}\\\
\end{bmatrix}\begin{bmatrix}\delta z_{1}\\\ \delta a_{1}\end{bmatrix}$

which is illustrated in Fig. 4(a). The two subsystems are coupled with weight
matrices

$\begin{bmatrix}a_{12}&b_{12}\\\ c_{12}&d_{12}\\\ \end{bmatrix}\ \text{and}\
\begin{bmatrix}a_{21}&b_{21}\\\ c_{21}&d_{21}\\\ \end{bmatrix}$ (10)

and the self-feedback loops are with weight matrices

$\begin{bmatrix}a_{11}&b_{11}\\\ c_{11}&d_{11}\\\ \end{bmatrix}\ \text{and}\
\begin{bmatrix}a_{22}&b_{22}\\\ c_{22}&d_{22}\\\ \end{bmatrix}$ (11)

Given a specific type of combinations, i.e., specified coupling matrices Eq.
(10), we can solve for $T_{y}$ and $T_{a}$ from Eq. (7) implicitly. In
Reference [14], basic types of combinations that can preserve stability are
provided and proved. Here we present two examples: hierarchical and feedback.

(a) Combination of two subsystems.

(b) Hierarchical combination

(c) Feedback combination

Figure 4: Combination examples.

Hierarchical combination: Comparing with the hierarchical form in Section 3.2
of [14] or the illustration in Fig. 4(b), it requires $T_{y}$ and $T_{a}$ to
cancel the upper right block:

$\begin{bmatrix}a_{12}&b_{12}\\\ c_{12}&d_{12}\\\ \end{bmatrix}=0$ (12)

Feedback combination: Comparing with the feedback form in Section 3.4 of [14]
or the illustration in Fig. 4(c), it requires $T_{y}$ and $T_{a}$ to make

$\begin{bmatrix}a_{12}&b_{12}\\\ c_{12}&d_{12}\\\
\end{bmatrix}=-\begin{bmatrix}a_{21}&b_{21}\\\ c_{21}&d_{21}\\\
\end{bmatrix}^{T}$ (13)

The existence of $T_{y}$ and $T_{a}$ needs further study. Intuitively, take
the hierarchical-type modularity as an example. There are $2n^{2}$ unknown
variables, that is, the scalar elements from $T_{y}$ and $T_{a}$, to satisfy
$2(n^{2}-n)$ equations from specifying the coupling matrices, which implies
there always exists a solution. Note that the above discussions are with the
assumption that $\dot{T}_{y}$ can be ignored (see Appendix. X-C1).

Provided that the modularity in control has been realized, the nonlinear
stability problem is then simplified to the stability of the self-feedback
parts of those subsystems, i.e., the self-feedback matrices Eq. (11). Solving
their characteristic equations for negative eigenvalues yields linear
constraints on network Jacobians for stability. Further simplification of the
constraints via existing controllers and network architectures are in Section
V-A, the implementation.

It is possible to solve the coordinate transformation algorithmically. In the
following, we provide an explicit solution to the hierarchical-type
modularity.

Figure 5: Convergence of hierarchical combinations, where
$\bm{\delta}_{i}=[\delta z_{i},\delta a_{i}]^{T}$ denotes the differential
displacements between two arbitrary trajectories at fixed time $t$ of the
subsystem. When $\bm{\delta}_{i}=0$, all possible trajectories of the $i$th
subsystem converge to its equilibrium. The subsystems converge recursively
following their hierarchies, given stable self-feedback-loops and bounded
couplings, as the lowest level with nonzero $\bm{\delta}$’s is always
independent.

###  II-D Explicit solution to the hierarchical-type modularity

As for the hierarchical-type modularity, the subsystems converge in the order
of their hierarchies, provided that each self-feedback loop is stable and the
couplings between subsystems are bounded (Fig. 5).

Stability Theorem (proof in Appendix X-C): For the system
$\mathbf{\dot{y}}=\mathbf{f}(\mathbf{y},\mathbf{u})$ with diagonalizable
$\partial\mathbf{f}/\partial\mathbf{y}$, given a neural control policy
$\bm{\pi}$, and two transformation layers $T_{y}$ and $T_{a}$ that (i) the
control commands $\mathbf{u}=T_{a}\bm{\pi}(T_{y}\mathbf{y})$, and (ii)
$T_{y}^{T}T_{y}$ and $(T_{a}^{-1})^{T}T_{a}^{-1}$ are uniformly positive
definite, if there exists $\alpha>0$ in a region such that $\forall\mathbf{x}$
and $\forall t>0$,

$\displaystyle\frac{\partial\pi^{i}}{\partial
s_{1i}}R_{ii}+\Lambda_{ii}<-\alpha,\ \ \ \frac{\partial\pi^{i}}{\partial
s_{2i}}R_{ii}<-\alpha$ (14)

where $R_{ii}$ and $\Lambda_{ii}$ are the $i$th diagonal components of
$R=T_{y}\frac{\partial\mathbf{f}}{\partial\mathbf{u}}T_{a}$ and
$\Lambda=T_{y}\frac{\partial\mathbf{f}}{\partial\mathbf{y}}T_{y}^{-1}$
respectively, the neural control system is contracting in the region, provided
that

  * •

$T_{y}$ is the eigenvector matrix from eigenvalue decomposition of
diagonalizable $\partial\mathbf{f}/\partial\mathbf{y}$.

  * •

$T_{a}=(PQ^{T})^{-1}$, where $Q$ is the left matrix from the QR decomposition
of $PT_{y}\partial\mathbf{f}/\partial\mathbf{u}$ and $P$ is the permutation
matrix with all ones at the skew diagonal.

  * •

each subsystem is controlled by an independent neural network,
$\bm{\pi}=[\pi^{1},\cdots,\pi^{n}]$, as illustrated in Fig. 6.

  * •

there exists $N\in\mathbb{R}^{+}$, $\forall t\geq N$, $\dot{T}_{y}=0$, and
$\forall t<N$, the dynamics are bounded.

Figure 6: Independent network $\pi^{i}$ for each dimension $i$.

##  III Signal composition

The convergence of latent signals does not necessarily imply the stability of
physical dynamics. Via composition of latent signals
$\mathbf{y}=\mathbf{g}(\mathbf{x},\mathbf{\dot{x}})$, one can enforce a
hierarchy of contracting systems again as follows:

$\displaystyle\mathbf{\dot{y}}-\mathbf{f}(\mathbf{y},\mathbf{u})=0$ (15)
$\displaystyle\mathbf{g}(\mathbf{x},\mathbf{\dot{x}})-\mathbf{y}=0$

where $\mathbf{f}$ denotes the resulting dynamics after creating the latent
signal $\mathbf{y}=\mathbf{g}(\mathbf{x},\mathbf{\dot{x}})$ for the task.

Provided that the latent signals are converging and the composition of latent
signals are stable, i.e.,
$\mathbf{\dot{y}}-\mathbf{f}(\mathbf{y},\mathbf{u})=0$ and
$\mathbf{g}(\mathbf{x},\mathbf{\dot{x}})=0$ are contracting, their combined
system Eq. (15) is contracting, that is, the stability of the physical
dynamics in the task space.

One example of such stable function composition is the composite variables
that mix the tracking errors in position and velocity,
$\mathbf{y}=K_{1}\mathbf{e}(t)+K_{2}\dot{\mathbf{e}}(t)$ where
$\mathbf{e}=\mathbf{x}_{d}-\mathbf{x}$ and $K_{1}$ and $K_{2}$ are positive
diagonal matrices. The solution to
$\mathbf{g}(\mathbf{x},\mathbf{\dot{x}})=K_{1}\mathbf{e}(t)+K_{2}\dot{\mathbf{e}}(t)=0$
is $\mathbf{e}=Cexp(-K_{1}K_{2}^{-1}t)$ where $C$ is a constant determined by
the initial condition. This solution converges exponentially.

It is possible to integrate with the latent signals in machine learning by
$\mathbf{y}=\bm{\pi}_{latent}(\mathbf{x},\mathbf{\dot{x}})$ with observable
$\mathbf{x}$ and $\mathbf{\dot{x}}$, which is not necessarily limited to the
physical states of end effectors in the task space. The stability constraints
on the input gradients of $\bm{\pi}_{latent}$ can be derived by applying
Theorem 2.1 in [30].

##  IV For unknown environments

Our method requires the dynamic model in the latent space, in particular, the
partial derivatives $\partial\mathbf{f}/\partial\mathbf{y}$ and
$\partial\mathbf{f}/\partial\mathbf{u}$ in Eq. (7). For unknown environments,
one may apply system identification to learn the dynamics,
$\dot{\mathbf{y}}=\mathbf{f}(\mathbf{y},\mathbf{u})$. Here we present another
solution that leverages task space controllers to create a model composed of
the controller gains for the applications that can be handled with information
from the end effectors.

Task space controllers in [15] can regulate the behaviors at the end effectors
to follow some reference $\mathbf{p}$, for example, acceleration-based
$\mathbf{p}=\mathbf{\ddot{x}}_{r}=\mathbf{\ddot{x}}_{d}+K_{d}(\mathbf{\dot{x}}_{d}-\mathbf{\dot{x}})+K_{p}(\mathbf{x}_{d}-\mathbf{x})$.
With the task space controller
$\bm{\tau}=G(\mathbf{q},\mathbf{\dot{q}},\mathbf{p})$, the dynamics of end
effectors can be represented as
$\mathbf{\ddot{e}}+K_{d}\mathbf{\dot{e}}+K_{p}\mathbf{e}=0$, where
$\mathbf{e}=\mathbf{x}_{d}-\mathbf{x}$.

Here we add the RL control into the reference, for example, acceleration-based
$\mathbf{p}=\mathbf{\ddot{x}}_{r}+\mathbf{u}$, leading to the dynamic model in
the task space
$\mathbf{\ddot{e}}+K_{d}\mathbf{\dot{e}}+K_{p}\mathbf{e}+\mathbf{u}=0$.

Applying the composite variable
$\mathbf{y}=K_{1}\mathbf{e}+K_{2}\mathbf{\dot{e}}$, one may create the latent-
space model following the equations in Appendix X-D, leading to the stability
constraints as simple as having positive input gradients of networks.

##  V Implementation and limitations

###  V-A Implementation

The pseudo-code of RL algorithms remains the same, while the graph for control
needs to be adapted as follows:

  1. 1.

Use the latent signals $\mathbf{y}=\mathbf{g}(\mathbf{x},\mathbf{\dot{x}})$ as
the RL state for control.

  2. 2.

Add two layers $T_{y}$ and $T_{a}$ in the control framework as in Fig. (1),
the values of which are updated time-stepwise according to the stability
theorem in Section. II-D

  3. 3.

Feed the accumulative error and the current tracking error to the control
networks in Fig. 6 that each dimension is controlled by an independent
network.

  4. 4.

Limit the network Jacobians $\partial\pi^{i}/\partial\mathbf{s}_{i}$
($i=1,2,\cdots,m$) by the linear inequality constraints in Eq. (14) where
$R_{ii}$ and $\Lambda_{ii}$ updated stepwise.

###  V-B Further simplification

When the dynamic without RL control is stable or has been stabilized by
existing controllers, i.e., $\sup(\Lambda_{ii})<0$, the linear constraints Eq.
(14) is further simplified into constraints on the signs of the network
Jacobians ($j=1,2$),

$\frac{\partial\pi^{i}}{\partial s_{ji}}R_{ii}<0$ (16)

Further simplification to the signs of network weights can be achieved via
network architectures.

Given MLPs with the hyperbolic tangent function as the activation function,
constraints on the signs of network Jacobians are further simplified to
constraints on the signs of MLP weights.

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\pi^{i}/\partial s_{1i}$ and
$\partial\pi^{i}/\partial s_{2i}$ are the products of weights and
$[1-\tanh^{2}{(\cdot)}]$. That is, let $W_{1,1}$ denote the weights of the
first layer corresponding to $s_{i1}$, the network Jacobian with respect to
$s_{i1}$ becomes

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

where $M_{j}=1-\tanh^{2}(W_{j}x_{j}+b_{j})$. Similarly, one can compute
$\partial\pi^{i}/\partial s_{2i}$. When the activation function is not
saturated $M_{j}\neq 0$, components of $M_{j}$ are positive. Hence, one can
bound the signs of $W_{j}$ to bound the signs of $\partial\pi^{i}/\partial
s_{ki}$ ($k=1,2$) according to Eq. (17) to satisfy the constraints in Eq.
(16).

###  V-C Limitations

##### Limited control frequency

When applying the hierarchical type of modularity in Section. II-D, the
control frequency is limited by the eigenvalue decomposition and QR
decomposition in solving for $T_{y}$ and $T_{a}$. One may design the
$\mathbf{y}=\mathbf{g}(\mathbf{x})$ that makes diagonal
$\partial\mathbf{f}/\partial\mathbf{y}$ for fast dynamics.

##### Loss of optimality

Optimality in control may be affected by forcing the coupling way in the
auxiliary space by forcing the function compositions of the latent signals. If
applying separate policies to control and planning, the loss of optimality in
control can be compensated in planning.

##  VI Example 1:
Necessity for robustness and generlization

Figure 7: Intuition. Local optimal policies can contain unstable control
because the stability boundaries may not align with the landscape of expected
returns.

This section briefly presents the main results in [6], which illustrates that
local optimal policies may contain unstable control for some dynamic
parameters and overfitting to such instabilities can deteriorate the
generalization of RL. Adding stability constraints can mitigate such concerns.

Intuitively, because the stability boundaries are functions of dynamic
parameters and the expected returns are functions of parameter distributions,
the boundaries do not necessarily align with the reward (expected returns)
landscape, which allows local optimal policies contain unstable control (Fig.
7).

The example is to move a 2D “peg” to touch a surface at the desired position
$x_{d}$ with the desired force $f_{d}$ by proximal policy optimization
(PPO)[31]. The unstable dynamics observed is oscillation and drift, the
fraction of which is increasing as learning progresses. Adding stability
constraints via the method here (denoted by C-PPO) can prevent such
instabilities and improve generalization (Fig. 8).

The implementation of stable modular control is in Appendix X-F. The latent
signals are the state $x$ and force $f$ and a nominal model of the dynamics is
provided.

##  VII Example 2: Manipulation learning

To evaluate our method, we created peg-maze and peg-push tasks following the
concepts of ant maze and ant push. A Franka Panda robot is used to move the
peg. The experiments are in simulation.

With the easier peg-maze task, we compare our method with the state of the
art, a hierarchical learning algorithm HIRO. HIRO plans and tracks in general
coordinates with added relabeling processes to alleviate the non-stationary
issue, outperforming other methods on ant push and ant maze, particularly in
data efficiency[17].

In our method, we extracted and modified the hierarchical learning framework
from the HIRO implementation in Tensorflow Model Garden. In particular, we
extracted the 2-level learning framework that includes two TD3 agents. The
high-level policy runs with a slower frequency to plan the desired state
increments; the low level agent runs faster to track the high-level action
within a fixed horizon. We removed the relabeling process in HIRO, used the
task space instead of the joint space for planning and tracking, and added the
composite errors and operational space controllers for unknown environments.

The resulting framework is illustrated in Fig. 9(a). The stable neural control
framework in Fig. 1 is implemented in the low level. Note that the integral
operator is applied to translation but not orientation errors, as the tasks
require strong control for accurate position but prefer flexible orientation.
Our method is denoted by TD3/C-TD3.

The ablation study is performed with the more difficult peg-push task, to
evaluate the effect from neural controllers and from stability constraints. We
compare our TD3/C-TD3 with

  * •

One-level learning with non-neural control:

TD3/Control and PPO/Control,

  * •

Hierarchical learning without stability constraints:

TD3/TD3.

###  VII-A TD3/C-TD3 vs. HIRO on peg maze.

##### Task setting.

The task is to reach a goal in the green zone through the unstructured
environment with fixed obstacles (top in Fig. 9(b)). Robot’s motion is limited
by a hood. Particularly, the vertical limit forces the peg going through the
unstructured zone.

##### Experiment

Both HIRO and TD3/C-TD3 can find the path with high success rates and with
high data efficiency ($\sim$ 2M and 1M respectively) given the initial
exploration in the direction of the maze. To evaluate the safety in
exploration, we tested policies for TD3/C-TD3 at 10k, 200k, 400k, 600k, 800k,
and 1M steps (training converged), and for HIRO at 10k, 500k, 1M, 1.5M, 2M,
and 2.3M steps (training converged); and we compare the distribution of
contact forces.

The distribution of contact forces is plotted in Fig. 9(c). In TD3/C-TD3,
$99.93\%$ steps are with less-than-1 N contacts, including $64.24\%$ steps
without contacts. In HIRO, $65.88\%$ steps are with less-than-1 N contacts,
including $6.22\%$ zero contacts. Trajectory examples are in
https://www.youtube.com/watch?v=nFLHwVfPIJw. This shows that TD3/C-TD3 can
explore through the obstacles with soft touches, i.e., small contact forces.

Besides the stable task-space neural control in the low level, TD3/C-TD3 also
benefits from an easier high-level planning problem: finding a path vs.
planning joint positions in HIRO.

###  VII-B Ablation study: peg-push task

##### Task setting

The peg-push task is to reach a hidden goal in a box (bottom in Fig. 9(b))
that requires interactions with moving objects. The robot needs to push open
the sliding lid, insert into a deep slot, and push away the green obstacle to
reach the small red spot. Appendix X-H lists the geometry and contact
settings.

##### Experiment

We firstly tuned the task space controller, the high-level action range, the
low-level horizon, and the fixed trajectory length to make sure that the robot
can stably explore the entire environment. Controller gains, RL
hyperparameters, and other task parameters are listed in Appendix X-H.

The average learning curve from 5 seeds are presented in Fig. 9(d). Only
TD3/C-TD3 is able to learn the task (expected return in $[-400,-200]$), and 2
trials have converged after 3M experiment steps. This causes the large
variance around 3M to 6M in the averaged learning curve. Other methods fail to
learn the task (returns $\approx-700$).

Close examination of the trajectories shows that TD3/control and PPO/control
can learn stable trajectories but haven’t figured out the entire path to the
goal within 12M steps. Without stability constraints, TD3/TD3 learns unstable
trajectories passing by the goal. Trajectory examples can be found at
https://www.youtube.com/watch?v=RpGzpZPifUw.

There still exists one question: whether the soft contact model highly
simplifies the task as we observed penetration (the penetration issue
discussed in [32]). Further study will focus on evaluation of the neural
controller for complex contact stiffness that varies across different surface
materials and changes along with relative motions.

##  VIII Concluding Remarks

Bioinspired modularity can deconstruct intrinsic difficulties and
computational complexities in control, planning, and learning. We realize
modular neural control for RL via contraction theory, leading to a minimally
invasive way for stability guarantees, which allows arguably easy integration
with the modular neural architectures in machine learning, in particular,
hierarchical RL, and improves its performance. This modular control framework
provides a new perspective to integrate control theoretical results with RL:
applying control techniques to signal composition and dynamic decomposition
for modularity. Towards full integration with the modular neural architectures
in machine learning, the concept applies to dynamics beyond physical robotics.

##  IX Appendix

The appendix at the end of this paper includes the following sections.

  1. A.

Incremental exponential stability.

  2. B.

Stability theorem for smooth coordinate transformation.

  3. C.

Proof of the theorem for the hierarchical combination. .

  4. D.

Integration of composite variables and task space controllers for unknown
environments.

  5. E.

Robustness.

  6. F.

Implementation in the PPO example

  7. G.

Peg maze example.

  8. H.

Peg push example.

(a) Peg touching by PPO.

(b) Stability: mitigation of oscillation and drift by constraints

(c) Generalization and robustness.

Figure 8: Necessity for robustness and generalization (more results in [6]).
(a) 2D peg-touching task. Both PPO and C-PPO can learn the task with similar
return levels and with similar data efficiency. (b) Stability. Starting from
the policy at Episode 5000, 8000 trajectories are tested every 1k episodes. In
PPO, the amount of oscillation and drift is increasing as learning progresses,
which implies overfitting, while in C-PPO, the fraction remains zero. (c)
Generalization and robustness. The policies at Episode 5000 are tested for
8000 trajectories, with new task parameters with $50\%$ differences in
$\tau_{i}$ ($i=x,z$), $100\%$ in $K_{sur}$ and $50\%$ in $K_{i}$ ($i=1,2$)
compared to the training set. The nominal model from the training set is used
in C-PPO. The histogram of tracking errors shows dramatic improvements in
generalization in C-PPO.

(a) Hierarchical RL with stable modular neural control in the low level
(TD3/C-TD3).

(b) Peg maze and peg push.

(c) Peg maze: Contact force distributions and the success rates.

(d) Peg push: Learning curves.

Figure 9: Effectiveness in manipulation learning (trajectory examples in
https://www.youtube.com/watch?v=nFLHwVfPIJw for peg maze and
https://www.youtube.com/watch?v=RpGzpZPifUw for peg push.) (a) Hierarchical RL
framework of TD3/C-TD3. (b) Peg maze and peg push tasks. Peg maze (top) is a
goal reaching task with fixed obstacles both near the goal (randomly sampled
in the green zone) and near the robot arm (the transparent hood). Peg push
(bottom) is a goal reaching task that requires opening a sliding lid box,
inserting into a deep slot, and pushing away the green obstacle to reach the
small red spot underneath. (c) Peg maze results: TD3/C-TD3 vs. HIRO (the
state-of-the-art data-efficient hierarchical RL). Both TD3/C-TD3 and HIRO can
quickly learn the task ($\sim$1M and 2M steps). The contact force distribution
are estimated from policies (success rates on top right) at 10k, 200k, 400k,
600k, 800k, and 1M steps for TD3/C-TD3, and at 10k, 500k, 1M, 1.5M, 2M, and
2.3M steps for HIRO. The distribution shows improvements in safe exploration
that in TD3/C-TD3, $99.93\%$ are with less than 1 N, and $65.88\%$ with HIRO.
(d) Peg push results: TD3/C-TD3 vs. RL with non-neural controllers
(PPO/Control, TD3/Control) vs. TD3/TD3 without stability constraints. Only
TD3/C-TD3 successfully learn the task within 12 M steps, with 2 trials learned
the task around 3 M steps causing the large variation. TD3/control and
PPO/control can learn stable trajectories but haven’t figured out the entire
path to the goal. Without stability constraints, TD3/TD3 learns unstable
trajectories passing by the goal.

## 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] R. Cheng, A. Verma, G. Orosz, S. Chaudhuri, Y. Yue, and J. Burdick, “Control regularization for reduced variance reinforcement learning,” in _International Conference on Machine Learning_. PMLR, 2019, pp. 1141–1150.
  * [2] F. Berkenkamp, M. Turchetta, A. P. Schoellig, and A. Krause, “Safe model-based reinforcement learning with stability guarantees,” _arXiv preprint arXiv:1705.08551_ , 2017.
  * [3] S. Singh, V. Sindhwani, J.-J. E. Slotine, and M. Pavone, “Learning stabilizable dynamical systems via control contraction metrics,” in _International Workshop on the Algorithmic Foundations of Robotics_. Springer, 2018, pp. 179–195.
  * [4] A. Mandlekar, Y. Zhu, A. Garg, L. Fei-Fei, and S. Savarese, “Adversarially robust policy learning: Active construction of physically-plausible perturbations,” in _2017 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_. IEEE, 2017, pp. 3932–3939.
  * [5] W. Lohmiller and J.-J. E. Slotine, “On contraction analysis for non-linear systems,” _Automatica_ , vol. 34, no. 6, pp. 683–696, 1998.
  * [6] B. Song, J.-J. Slotine, and Q.-C. Pham, “Example when local optimal policies contain unstable control,” _arXiv preprint arXiv:2209.07324_ , 2023.
  * [7] C. Dawson, S. Gao, and C. Fan, “Safe control with learned certificates: A survey of neural lyapunov, barrier, and contraction methods,” _arXiv preprint arXiv:2202.11762_ , 2022.
  * [8] N. Boffi, S. Tu, N. Matni, J.-J. Slotine, and V. Sindhwani, “Learning stability certificates from data,” in _Conference on Robot Learning_. PMLR, 2021, pp. 1341–1350.
  * [9] Z. Qin, K. Zhang, Y. Chen, J. Chen, and C. Fan, “Learning safe multi-agent control with decentralized neural barrier certificates,” _arXiv preprint arXiv:2101.05436_ , 2021.
  * [10] Y. Luo and T. Ma, “Learning barrier certificates: Towards safe reinforcement learning with zero training-time violations,” _Advances in Neural Information Processing Systems_ , vol. 34, pp. 25 621–25 632, 2021.
  * [11] H. Ma, C. Liu, S. E. Li, S. Zheng, and J. Chen, “Joint synthesis of safety certificate and safe control policy using constrained reinforcement learning,” in _Learning for Dynamics and Control Conference_. PMLR, 2022, pp. 97–109.
  * [12] S. Reed, K. Zolna, E. Parisotto, S. G. Colmenarejo, A. Novikov, G. Barth-Maron, M. Gimenez, Y. Sulsky, J. Kay, J. T. Springenberg _et al._ , “A generalist agent,” _arXiv preprint arXiv:2205.06175_ , 2022.
  * [13] J.-J. Slotine and W. Lohmiller, “Modularity, evolution, and the binding problem: a view from stability theory,” _Neural networks_ , vol. 14, no. 2, pp. 137–145, 2001.
  * [14] J.-J. E. Slotine, “Modular stability tools for distributed computation and control,” _International Journal of Adaptive Control and Signal Processing_ , vol. 17, no. 6, pp. 397–416, 2003.
  * [15] 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.
  * [16] S. Fujimoto, H. Hoof, and D. Meger, “Addressing function approximation error in actor-critic methods,” in _International conference on machine learning_. PMLR, 2018, pp. 1587–1596.
  * [17] O. Nachum, S. S. Gu, H. Lee, and S. Levine, “Data-efficient hierarchical reinforcement learning,” _Advances in neural information processing systems_ , vol. 31, 2018.
  * [18] E. M. Aylward, P. A. Parrilo, and J.-J. E. Slotine, “Stability and robustness analysis of nonlinear systems via contraction metrics and sos programming,” _Automatica_ , vol. 44, no. 8, pp. 2163–2170, 2008.
  * [19] 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.
  * [20] H. Tsukamoto, S.-J. Chung, and J.-J. Slotine, “Learning-based adaptive control via contraction theory,” in _IEEE CDC_ , 2021.
  * [21] H. Tsukamoto and S.-J. Chung, “Learning-based robust motion planning with guaranteed stability: A contraction theory approach,” _IEEE Robotics and Automation Letters_ , vol. 6, no. 4, pp. 6164–6171, 2021.
  * [22] 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.
  * [23] S. Lale, Y. Shi, G. Qu, K. Azizzadenesheli, A. Wierman, and A. Anandkumar, “Kcrl: Krasovskii-constrained reinforcement learning with guaranteed stability in nonlinear dynamical systems,” _arXiv preprint arXiv:2206.01704_ , 2022.
  * [24] D. Sun, S. Jha, and C. Fan, “Learning certified control using contraction metric,” _arXiv preprint arXiv:2011.12569_ , 2020.
  * [25] M. Jin and J. Lavaei, “Stability-certified reinforcement learning: A control-theoretic perspective,” _IEEE Access_ , vol. 8, pp. 229 086–229 100, 2020.
  * [26] K. A. Thoroughman and R. Shadmehr, “Learning of action through adaptive combination of motor primitives,” _Nature_ , vol. 407, no. 6805, pp. 742–747, 2000.
  * [27] M. Santello, M. Bianchi, M. Gabiccini, E. Ricciardi, G. Salvietti, D. Prattichizzo, M. Ernst, A. Moscatelli, H. Jörntell, A. M. Kappers _et al._ , “Hand synergies: Integration of robotics and neuroscience for understanding the control of biological and artificial hands,” _Physics of life reviews_ , vol. 17, pp. 1–23, 2016.
  * [28] N. D. Ratliff, J. Issac, D. Kappler, S. Birchfield, and D. Fox, “Riemannian motion policies,” _arXiv preprint arXiv:1801.02854_ , 2018.
  * [29] J. Pfeiffer, S. Ruder, I. Vulić, and E. M. Ponti, “Modular deep learning,” _arXiv preprint arXiv:2302.11529_ , 2023.
  * [30] 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.
  * [31] J. Schulman, F. Wolski, P. Dhariwal, A. Radford, and O. Klimov, “Proximal policy optimization algorithms,” _arXiv preprint arXiv:1707.06347_ , 2017\.
  * [32] M. Parmar, M. Halm, and M. Posa, “Fundamental challenges in deep learning for stiff contact dynamics,” in _2021 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_. IEEE, 2021, pp. 5181–5188.

##  X Appendix

###  X-A Incremental exponential stability

The stability here in this paper refers to the incremental exponential
stability as follows.

Definition[30]. Let $\xi_{1}(t)$ and $\xi_{2}(t)$ denote two trajectories of a
system $\dot{\mathbf{x}}=\mathbf{g}(\mathbf{x},t)$. This system
$\dot{\mathbf{x}}=\mathbf{g}(\mathbf{x},t)$ is incrementally exponentially
stable if there exist $C>0$ and $\beta>0$ such that

$\lVert\xi_{1}(t)-\xi_{2}(t)\rVert\leq Ce^{-\beta
t}\lVert\xi_{1}(0)-\xi_{2}(0)\rVert$ (18)

Note that stability of nonlinear systems can be defined in multiple ways and
the differences are discussed in [30].

###  X-B Stability theorem for smooth coordinate transformation

Theorem: Given a dynamic system
$\dot{\mathbf{y}}=\mathbf{f}(\mathbf{y},\mathbf{u},t)$, a neural control
policy $\bm{\pi}$, and two transformation layers $T_{y}$ and $T_{a}$ that (i)
the control commands $\mathbf{u}=T_{a}\bm{\pi}(T_{y}\mathbf{y})$, and (ii)
$T_{y}^{T}T_{y}$ and $(T_{a}^{-1})^{T}T_{a}^{-1}$ are uniformly positive
definite, this neural control system is incrementally exponentially stable in
a region, if $\forall t>0$, $\forall\mathbf{y}$ in the region, there exists
$\beta>0$ such that

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

where

$F_{1}=\begin{bmatrix}T_{y}\frac{\partial\mathbf{f}}{\partial\mathbf{y}}T_{y}^{-1}&T_{y}\frac{\partial\mathbf{f}}{\partial\mathbf{u}}T_{a}\\\\[8.0pt]
\frac{\partial\bm{\pi}}{\partial\mathbf{s}_{1}}T_{y}\frac{\partial\mathbf{f}}{\partial\mathbf{y}}T_{y}^{-1}+\frac{\partial\bm{\pi}}{{\partial\mathbf{s}_{2}}}&\frac{\partial\bm{\pi}}{\partial\mathbf{s}_{1}}T_{y}\frac{\partial\mathbf{f}}{\partial\mathbf{u}}T_{a}\end{bmatrix}$
(20)

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

$F_{2}=\begin{bmatrix}\dot{T}_{y}T_{y}^{-1}&0\\\\[8.0pt]
\frac{\partial\bm{\pi}}{\partial\mathbf{s}_{1}}\dot{T}_{y}T_{y}^{-1}&0\end{bmatrix}$
(21)

and
$\mathbf{s}=[\mathbf{s}_{1}^{T},\mathbf{s}_{2}^{T}]^{T}=[\mathbf{z}^{T},\int_{0}^{t}\mathbf{z}^{T}dt]^{T}$
with $\mathbf{z}=T_{y}\mathbf{y}$.

Proof:

##### Differential dynamics in the auxiliary space

Provided that $\mathbf{z}=T_{y}\mathbf{y}$ and
$\mathbf{a}=T_{a}^{-1}\mathbf{u}$, the differential dynamics in the auxiliary
space is Eq. (6), which can also be represented as

$\begin{bmatrix}\delta\dot{\mathbf{z}}\\\
\delta\dot{\mathbf{a}}\end{bmatrix}=F\begin{bmatrix}\delta\mathbf{z}\\\
\delta\mathbf{a}\end{bmatrix}$ (22)

where

$F=\begin{bmatrix}[\dot{T}_{y}+T_{y}\frac{\partial\mathbf{f}}{\partial\mathbf{y}}]T_{y}^{-1}&T_{y}\frac{\partial\mathbf{f}}{\partial\mathbf{u}}T_{a}\\\\[8.0pt]
F_{21}&\frac{\partial\bm{\pi}}{\partial\mathbf{s}_{1}}T_{y}\frac{\partial\mathbf{f}}{\partial\mathbf{u}}T_{a}\end{bmatrix}$
(23)

and

$F_{21}=\frac{\partial\bm{\pi}}{\partial\mathbf{s}_{1}}(\dot{T}_{y}+T_{y}\frac{\partial\mathbf{f}}{\partial\mathbf{y}})T_{y}^{-1}+\frac{\partial\bm{\pi}}{{\partial\mathbf{s}_{2}}}$
(24)

##### Sufficient condition that makes the squared distance between
trajectories exponentially converge in the auxiliary space

The time derivative of the squared distance in the auxiliary space can be
written as

$\displaystyle\frac{d}{dt}\Big{(}\begin{bmatrix}\delta\mathbf{z}\\\
\delta\mathbf{a}\end{bmatrix}^{T}\begin{bmatrix}\delta\mathbf{z}\\\
\delta\mathbf{a}\end{bmatrix}\Big{)}=2\begin{bmatrix}\delta\mathbf{z}\\\
\delta\mathbf{a}\end{bmatrix}^{T}\frac{d}{dt}\begin{bmatrix}\delta\mathbf{z}\\\
\delta\mathbf{a}\end{bmatrix}$ (25) $\displaystyle=$ $\displaystyle
2\begin{bmatrix}\delta\mathbf{z}\\\
\delta\mathbf{a}\end{bmatrix}^{T}F\begin{bmatrix}\delta\mathbf{z}\\\
\delta\mathbf{a}\end{bmatrix}\leq
2\lambda_{max}(F)\Big{(}\begin{bmatrix}\delta\mathbf{z}\\\
\delta\mathbf{a}\end{bmatrix}^{T}\begin{bmatrix}\delta\mathbf{z}\\\
\delta\mathbf{a}\end{bmatrix}\Big{)}$

where $\lambda_{max}(F)$ represents the function that takes the largest
eigenvalue of $\frac{1}{2}(F^{T}+F)$.

If $\lambda_{max}(F)$ is uniformly strictly negative, that is,
$\exists\beta>0$, $\forall\mathbf{x},\forall
t>0,\lambda_{max}(F)\leq-\beta<0$,

$\displaystyle(\delta\mathbf{z})^{T}(\delta\mathbf{z})+(\delta\mathbf{a})^{T}(\delta\mathbf{a})$
(26) $\displaystyle\leq$ $\displaystyle
e^{\int_{0}^{t}2\lambda_{max}(F)dt}[(\delta\mathbf{z}_{0})^{T}(\delta\mathbf{z}_{0})+(\delta\mathbf{a}_{0})^{T}(\delta\mathbf{a}_{0})]$
$\displaystyle=$ $\displaystyle e^{-2\beta
t}[(\delta\mathbf{z}_{0})^{T}(\delta\mathbf{z}_{0})+(\delta\mathbf{a}_{0})^{T}(\delta\mathbf{a}_{0})]$

where $\delta\mathbf{z}_{0}$ and $\delta\mathbf{a}_{0}$ denotes the initial
values at $t=0$

##### Distance between two arbitrary trajectories

The uniformly positive definite $T_{y}^{T}T_{y}$ means that
$\exists\underline{m}>0$, $\forall\mathbf{y},\forall t>0$,
$T_{y}^{T}T_{y}\succeq\underline{m}I$. Similarly, the uniformly positive
definite $(T_{a}^{-1})^{T}T_{a}^{-1}$ makes
$(\delta\mathbf{u})^{T}(T_{a}^{-1})^{T}T_{a}^{-1}(\delta\mathbf{u})>0$.

Given two arbitrary trajectories $\xi_{0}$ and $\xi_{1}$, the distance is
defined as the smallest path integral between two points at some fixed time
$t$, that is,

$\displaystyle\|\xi_{1}(t)-\xi_{0}(t)\|=\int_{\xi_{0}}^{\xi_{1}}\|\delta\mathbf{y}\|$
(27) $\displaystyle=$
$\displaystyle\int_{\xi_{0}}^{\xi_{1}}\frac{\sqrt{(\delta\mathbf{y})^{T}\underline{m}I\delta\mathbf{y}}}{\sqrt{\underline{m}}}$
$\displaystyle\leq$
$\displaystyle\int_{\xi_{0}}^{\xi_{1}}\frac{\sqrt{(\delta\mathbf{y})^{T}T_{y}^{T}T_{y}\delta\mathbf{y}}}{\sqrt{\underline{m}}}$
$\displaystyle\leq$
$\displaystyle\int_{\xi_{0}}^{\xi_{1}}\frac{\sqrt{(\delta\mathbf{y})^{T}T_{y}^{T}T_{y}\delta\mathbf{y}+(\delta\mathbf{u})^{T}(T_{a}^{-1})^{T}T_{a}^{-1}(\delta\mathbf{u})}}{\sqrt{\underline{m}}}$
$\displaystyle=$
$\displaystyle\int_{\xi_{0}}^{\xi_{1}}\sqrt{\frac{[(\delta\mathbf{z})^{T}(\delta\mathbf{z})+(\delta\mathbf{a})^{T}(\delta\mathbf{a})]}{\underline{m}}}$
$\displaystyle\leq$
$\displaystyle\int_{\xi_{0}}^{\xi_{1}}\sqrt{\frac{e^{-2\beta
t}[(\delta\mathbf{z}_{0})^{T}(\delta\mathbf{z}_{0})+(\delta\mathbf{a}_{0})^{T}(\delta\mathbf{a}_{0})]}{\underline{m}}}$
$\displaystyle=$ $\displaystyle Ce^{-\beta
t}\lVert\xi_{1}(0)-\xi_{2}(0)\rVert$

where
$C=\int_{\xi_{0}}^{\xi_{1}}\sqrt{\frac{[(\delta\mathbf{z}_{0})^{T}(\delta\mathbf{z}_{0})+(\delta\mathbf{a}_{0})^{T}(\delta\mathbf{a}_{0})}{\underline{m}(\delta\mathbf{x}_{0})^{T}(\delta\mathbf{x}_{0})}}$
that is a constant determined by the initial values at $t=0$.

Therefore, if $\lambda_{max}(F)$ is uniformly strictly negative in some
region, the feedback-control system is incrementally exponentially stable.

##### Boundary for the constant part in $F$

We rewrite $F$, Eq. (23), into the sum of the constant $F_{1}$, Eq. (20), and
the time-varying $F_{2}$, Eq. (21). Let $\nu^{+}$ denote
$\lambda_{max}(F_{2})$. If there exists $\beta>0$,
$\lambda_{max}(F_{1})<-[\beta+\max(\nu^{+},0)]$,

$\displaystyle\lambda_{max}(F)<$
$\displaystyle\lambda_{max}(F_{1})+\lambda_{max}(F_{2})$ (28) $\displaystyle<$
$\displaystyle-[\beta+\max(\nu^{+},0)]+\nu^{+}$ $\displaystyle\leq$
$\displaystyle-\beta$

Hence, if $\forall t>0$, $\forall\mathbf{y}$ in a region, $\exists\beta>0$,
$\lambda_{max}(F_{1})<-[\beta+\max(\nu^{+},0)]$, the feedback-control system
in Fig. 1 is incrementally exponentially stable in the region. This region is
a contraction region.

Proof is finished.

###  X-C Proof of the theorem for the hierarchical combination.

####  X-C1 Assumption that the term involving $F_{2}$ in Eq. (21) can be seen
as bounded temporary disturbances

The assumption that the term involving $F_{2}$ in Eq. (21) can be seen as
bounded temporary disturbances means there exists $N>0$,

  * •

$\forall t\geq N$, $\dot{T}_{y}=0$

  * •

$\forall t<N$, there exists upper bound $C_{N}$ of the distance between two
arbitrary trajectories $\lVert\xi_{1}(t)-\xi_{2}(t)\rVert<C_{N}$

####  X-C2 Differential dynamics in the auxiliary space

With the above assumption, the dynamic stability in the auxiliary space Eq.
(22) can be analyzed with the model

$\begin{bmatrix}\delta\dot{\mathbf{z}}\\\
\delta\dot{\mathbf{a}}\end{bmatrix}=F_{1}\begin{bmatrix}\delta\mathbf{z}\\\
\delta\mathbf{a}\end{bmatrix}$ (29)

because given Eq. (29) is contracting with the rate $\beta>0$ in Eq. (18), Eq.
(22) is also contracting with the rate $\beta$:

  * •

$\forall t\geq N$, Eq. (22) with $\dot{T}_{y}=0$ becomes Eq. (29),

  * •

$\forall t<N$, the distance between two arbitrary trajectories is also
bounded, $\lVert\xi_{1}(t)-\xi_{2}(t)\rVert<(C_{N}e^{\beta N})e^{-\beta t}$,
where $C_{N}$ denotes the upper bound of the distance for $t<N$.

####  X-C3 Diagonal $\Lambda$ via $T_{y}$

Given diagonalizable $\partial\mathbf{f}/\partial\mathbf{y}$ and $T_{y}$ being
the eigenvector matrix,

$\Lambda=T_{y}\frac{\partial\mathbf{f}}{\partial\mathbf{y}}T_{y}^{-1}=\begin{bmatrix}\Lambda_{11}&&&\\\
&\Lambda_{22}&&\\\ &&\ddots&\\\ &&&\Lambda_{nn}\end{bmatrix}$ (30)

where $\Lambda_{ii}$ denotes the eigenvalues.

####  X-C4 Upper triangular $R$ via $T_{a}$

The term $R$ equals to $PR_{qr}^{T}P$, which is an upper triangular matrix
that can be represented as

$R=T_{y}\frac{\partial\mathbf{f}}{\partial\mathbf{u}}T_{a}=\begin{bmatrix}R_{11}&R_{12}&\cdots&R_{1n}\\\
&R_{22}&&R_{2n}\\\ &&\ddots&\vdots\\\ &&&R_{nn}\end{bmatrix}$ (31)

because

$\displaystyle R=$ $\displaystyle
T_{y}\frac{\partial\mathbf{f}}{\partial\mathbf{u}}T_{a}$ (32) $\displaystyle=$
$\displaystyle(PP)T_{y}\frac{\partial\mathbf{f}}{\partial\mathbf{u}}(PQ^{T})^{-1}$
$\displaystyle=$ $\displaystyle
P[(T_{y}\frac{\partial\mathbf{f}}{\partial\mathbf{u}})^{T}P]^{T}(PQ^{T})^{-1}$
$\displaystyle=$ $\displaystyle P(QR_{qr})^{T}(Q^{T})^{-1}P$ $\displaystyle=$
$\displaystyle PR_{qr}^{T}Q^{T}(Q^{T})^{-1}P$ $\displaystyle=$ $\displaystyle
PR_{qr}^{T}P$

since $T_{a}=(PQ^{T})^{-1}$, $Q$ is the left matrix in the QR decomposition
$QR_{qr}=(T_{y}\frac{\partial\mathbf{f}}{\partial\mathbf{u}})^{T}P$ with

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

####  X-C5 Hierarchical combination of subsystems

Plugging
$\Lambda=T_{y}\frac{\partial\mathbf{f}}{\partial\mathbf{x}}T_{y}^{-1}$ and
$R=T_{y}\frac{\partial\mathbf{f}}{\partial\mathbf{u}}T_{a}$ into Eq. (29), the
dynamics then can be written as

$\displaystyle\begin{bmatrix}\delta\dot{\mathbf{z}}\\\
\delta\dot{\mathbf{a}}\end{bmatrix}$
$\displaystyle=\begin{bmatrix}\Lambda&R\\\\[8.0pt]
\frac{\partial\bm{\pi}}{\partial\mathbf{s}_{1}}\Lambda+\frac{\partial\bm{\pi}}{{\partial\mathbf{s}_{2}}}&\frac{\partial\bm{\pi}}{\partial\mathbf{s}_{1}}R\end{bmatrix}\begin{bmatrix}\delta\mathbf{z}\\\
\delta\mathbf{a}\end{bmatrix}$ (34)

where

$\displaystyle\Lambda$
$\displaystyle=T_{y}\frac{\partial\mathbf{f}}{\partial\mathbf{y}}T_{y}^{-1}=\begin{bmatrix}\Lambda_{11}&&&\\\
&\Lambda_{22}&&\\\ &&\ddots&\\\ &&&\Lambda_{nn}\end{bmatrix}$ (35)
$\displaystyle R$
$\displaystyle=T_{y}\frac{\partial\mathbf{f}}{\partial\mathbf{u}}T_{a}=\begin{bmatrix}R_{11}&R_{12}&\cdots&R_{1n}\\\
&R_{22}&&R_{2n}\\\ &&\ddots&\vdots\\\ &&&R_{nn}\end{bmatrix}$

With the policy that is made of a sequence of networks in parallel in Fig. 6,
the dynamics Eq. (34) can be written explicitly with $F_{1}$ being represented
by two column blocks ($F_{1}=[F_{1l}\ F_{1r}]$) as follows:

$\begin{bmatrix}\delta\dot{z}_{1}\\\ \delta\dot{z}_{2}\\\ \vdots\\\
\delta\dot{z}_{n}\\\ \delta\dot{a}_{1}\\\ \delta\dot{a}_{2}\\\ \vdots\\\
\delta\dot{a}_{n}\end{bmatrix}=\begin{bmatrix}F_{1l}&F_{1r}\end{bmatrix}\begin{bmatrix}\delta
z_{1}\\\ \delta z_{2}\\\ \vdots\\\ \delta z_{n}\\\ \delta a_{1}\\\ \delta
a_{2}\\\ \vdots\\\ \delta a_{n}\end{bmatrix}$ (36)

where

$\displaystyle F_{1l}=\begin{bmatrix}\Lambda_{11}&&\\\ &\ddots&\\\
&&\Lambda_{nn}\\\\[10.0pt] \frac{\partial\pi^{1}}{\partial
s_{11}}\Lambda_{11}+\frac{\partial\pi^{1}}{\partial s_{21}}&&\\\ &\ddots&\\\
&&\frac{\partial\pi^{n}}{\partial
s_{1n}}\Lambda_{nn}+\frac{\partial\pi^{n}}{\partial s_{2n}}\\\ \end{bmatrix}$
(37) $\displaystyle F_{1r}=\begin{bmatrix}R_{11}&R_{12}&\cdots&R_{1n}\\\
&R_{22}&&R_{2n}\\\ &&\ddots&\vdots\\\ &&&R_{nn}\\\\[10.0pt]
\frac{\partial\pi^{1}}{\partial s_{11}}R_{11}&\frac{\partial\pi^{1}}{\partial
s_{11}}R_{12}&\cdots&\frac{\partial\pi^{1}}{\partial s_{11}}R_{1n}\\\
&\frac{\partial\pi^{2}}{\partial
s_{12}}R_{22}&&\frac{\partial\pi^{2}}{\partial s_{12}}R_{2n}\\\
&&\ddots&\vdots\\\ &&&\frac{\partial\pi^{n}}{\partial s_{1n}}R_{nn}\\\
\end{bmatrix}$

By grouping $\delta z_{i}$ and $\delta a_{i}$ from the same dimension, Eq.
(34) is rearranged into the hierarchical combination, i.e.,

$\displaystyle\begin{bmatrix}\delta\dot{z}_{1}\\\ \delta\dot{a}_{1}\\\
\delta\dot{z}_{2}\\\ \delta\dot{a}_{2}\\\ \vdots\\\ \delta\dot{z}_{n}\\\
\delta\dot{a}_{n}\end{bmatrix}=$
$\displaystyle\begin{bmatrix}F_{11}&F_{12}&\cdots&F_{1n}\\\
0&F_{22}&&F_{2n}\\\ \vdots&&\ddots&\vdots\\\
0&\cdots&0&F_{nn}\end{bmatrix}\begin{bmatrix}\delta z_{1}\\\ \delta a_{1}\\\
\delta z_{2}\\\ \delta a_{2}\\\ \vdots\\\ \delta z_{n}\\\ \delta
a_{n}\end{bmatrix}$ (38)

where

$\displaystyle F_{ii}$ $\displaystyle=\begin{bmatrix}\Lambda_{ii}&R_{ii}\\\
\frac{\partial\pi^{i}}{\partial
s_{1i}}\Lambda_{ii}+\frac{\partial\pi^{i}}{\partial
s_{2i}}&\frac{\partial\pi^{i}}{\partial s_{1i}}R_{ii}\end{bmatrix}$ (39)
$\displaystyle F_{ij}$ $\displaystyle=\begin{bmatrix}0&R_{ij}\\\
0&\frac{\partial\pi^{i}}{\partial s_{1i}}R_{ij}\end{bmatrix},\ (i<j)$

####  X-C6 Linear constraints from the preservation of contraction

Given the hierarchical combination Eq. (34), if

$\begin{bmatrix}\delta\dot{z}_{i}\\\
\delta\dot{a}_{i}\end{bmatrix}=F_{ii}\begin{bmatrix}\delta z_{i}\\\ \delta
a_{i}\end{bmatrix}$ (40)

is contracting, the resulting combination is contracting. The proof of the
preservation of contraction can be found in Section 3.8.3 of [5], Theorem 2.7
in [30], and in [13],.

If $F_{ii}$ has real negative eigenvalues, denoted by $-\lambda_{i}$ where
$\lambda_{i}>0$ ($i=1,2$), the solution of Eq. (40) converges to zero as $t$
progresses, i.e.,

$\begin{bmatrix}\delta z_{i}\\\ \delta
a_{i}\end{bmatrix}=c_{1}\mathbf{v}_{1}e^{-\lambda_{1}t}+c_{2}\mathbf{v}_{2}e^{-\lambda_{2}t}$
(41)

where $c_{1}$ and $c_{2}$ are constants, $\mathbf{v}_{1}$ and $\mathbf{v}_{2}$
are the corresponding eigenvectors.

By solving the characteristic equations of $F_{ii}$,

$s^{2}-(\Lambda_{ii}+\frac{\partial\pi^{i}}{\partial
s_{1i}}R_{ii})s-R_{ii}\frac{\partial\pi^{i}}{\partial s_{2i}}=0$ (42)

we obtain that if there exists $\alpha>0$ in a region such that
$\forall\mathbf{x}$ and $\forall t>0$,

$\displaystyle\frac{\partial\pi^{i}}{\partial
s_{1i}}R_{ii}+\Lambda_{ii}<-\alpha,\ \ \ \frac{\partial\pi^{i}}{\partial
s_{2i}}R_{ii}<-\alpha$ (43)

the matrix $F_{ii}$ has real negative eigenvalues, and hence Eq. (40) is
contracting, resulting in the contraction of their hierarchical combination.

Proof is finished.

###  X-D Integration of composite variables and task space controllers for
unknown environments

####  X-D1 Reference $\mathbf{p}(\mathbf{x},\mathbf{u})$ in task space control

Given rigid body dynamics

$H(\mathbf{q})\ddot{\mathbf{q}}+C(\mathbf{q},\dot{\mathbf{q}})\dot{\mathbf{q}}+\mathbf{g}(\mathbf{q})=\bm{\tau}$
(44)

where $\mathbf{q}$ represents the generalized coordinates, $H(\mathbf{q})$ is
the mass matrix, $C(\mathbf{q},\dot{\mathbf{q}})$ represents the Coriolis
terms, $\mathbf{g}(\mathbf{q})$ is the gravity, the dynamics of the end
effector can be represented as

$\displaystyle\mathbf{x}$ $\displaystyle=\mathbf{h}(\mathbf{q})$ (45)
$\displaystyle\dot{\mathbf{x}}$ $\displaystyle=J(\mathbf{q})\dot{\mathbf{q}}$
$\displaystyle\ddot{\mathbf{x}}$
$\displaystyle=J(\mathbf{q})\ddot{\mathbf{q}}+\dot{J}(\mathbf{q})\dot{\mathbf{q}}$

Given desired $\mathbf{x}_{d}$, one example of task space controllers without
considering the null space can be represented as

$\displaystyle\bm{\tau}=$ $\displaystyle
H(\mathbf{q})J^{\dagger}\mathbf{p}(\mathbf{x},\mathbf{u})$ (46)
$\displaystyle+H(\mathbf{q})J^{\dagger}(\ddot{\mathbf{x}}_{d}-\dot{J}(\mathbf{q})\dot{\mathbf{q}})$
$\displaystyle+C(\mathbf{q},\dot{\mathbf{q}})\dot{\mathbf{q}}+\mathbf{g}(\mathbf{q})$

the dynamic of the end effector becomes

$\mathbf{\ddot{x}}_{d}-\mathbf{\ddot{x}}+\mathbf{p}(\mathbf{x},\mathbf{u})=0$
(47)

The dynamic behavior of end effectors is regulated by
$\mathbf{p}(\mathbf{x},\mathbf{u})$. For example, one may design

$\mathbf{p}(\mathbf{x},\mathbf{u})=\Lambda_{d}^{-1}K_{d}(\mathbf{\dot{x}}_{d}-\mathbf{\dot{x}})+\Lambda_{d}^{-1}K_{p}(\mathbf{x}_{d}-\mathbf{x})+\Lambda_{d}^{-1}\mathbf{u}$
(48)

yielding the dynamics

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

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

Note that there exists different types of task space control in Section 3.2.3
in [15].

####  X-D2 Composite variables for a given $\mathbf{p}(\mathbf{x},\mathbf{u})$

Given Eq. (49), we construct a composite variable,

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

where $K_{1}$ and $K_{2}$ are positive diagonal matrices to make the solution
to $K_{1}\mathbf{e}+K_{2}\dot{\mathbf{e}}=0$ converging.

We can solve for $K_{1}$ and $K_{2}$ that plugging Eq. (50) into Eq. (49)
yields a dynamic model in the latent space in the form of

$A\dot{\mathbf{y}}+B\mathbf{y}+\mathbf{u}=0$ (51)

where $A$ and $B$ are positive-definite diagonal matrices.

The resulting solutions are as follows:

  * •

When $K_{p}=0$, we can use $\mathbf{y}=\mathbf{\dot{e}}$, yielding

$\Lambda_{d}\dot{\mathbf{y}}+K_{d}\mathbf{y}+\mathbf{u}=0$ (52)

  * •

Otherwise, we can use $\mathbf{y}=K_{1}\mathbf{e}+K_{2}\dot{\mathbf{e}}$ where

$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}}}$
(53)

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$, which should be realized in the
task space controller design. The resulting dynamics becomes
$A\dot{\mathbf{y}}+B\mathbf{y}+\mathbf{u}=0$ with

$A_{ii}=\frac{\lambda_{d,i}}{2K_{p,i}}(K_{d,i}\pm\sqrt{K_{d,i}^{2}-4K_{p,i}\lambda_{d,i}}),B_{ii}=\lambda_{d,i}$
(54)

####  X-D3 Stability constraints

Provided $A$ and $B$ are positive-definite diagonal matrices, the derivatives
$\frac{\partial\mathbf{f}}{\partial\mathbf{y}}=-A^{-1}B$ and
$\frac{\partial\mathbf{f}}{\partial\mathbf{u}}=-A^{-1}$ are negative-definite
diagonal matrices. Hence in Theorem 2, we can use identity $T_{y}$ and
$T_{a}$, resulting in the inequality condition that if there exists $\alpha>0$
such that

$\displaystyle-\frac{\partial\pi^{i}}{\partial
s_{1i}}A^{-1}_{ii}-A^{-1}_{ii}B_{ii}<-\alpha,\ \ \
-\frac{\partial\pi^{i}}{\partial s_{2i}}A^{-1}_{ii}<-\alpha$ (55)

the system is contracting.

Given positive-definite $A$ and $B$, the inequalities Eq. (55) becomes
positive network Jacobians.

###  X-E Robustness

####  X-E1 Robustness to unknown deterministic perturbation

The robustness to unknown deterministic perturbation can be found in Theorem
2.4 in [30], which states the robustness of a contracting system as follows:
Let $\xi_{0}$ denote a trajectory of a contracting system
$\dot{\mathbf{y}}=\mathbf{g}(\mathbf{y},t)$ with the convergence rate $\beta$,
and $\xi_{1}$ a trajectory of its perturbed system
$\dot{\mathbf{y}}=\mathbf{g}(\mathbf{y},t)+\mathbf{d}(\mathbf{y},t)$. The
distance between $\xi_{0}$ and $\xi_{1}$ exponentially converges to a bounded
error ball as long as $\Theta\mathbf{d}$ is bounded, where $\Theta$ is the
smooth coordinate transformation that constructs the contraction metric
$M=\Theta^{T}\Theta$. Specifically, if
$\exists\underline{m},\bar{m}\in\mathbb{R}_{>0}$ and
$\exists\bar{d}\in\mathbb{R}_{\geq 0}$ s.t.
$\bar{d}=\sup_{\mathbf{y},t}\lVert\mathbf{d}(\mathbf{y},t)\rVert$ and
$\underline{m}I\preceq M(\mathbf{y},t)\preceq\bar{m}I$, $\forall x,t$, then we
have the following relation:

$\lVert\xi_{1}(t)-\xi_{0}(t)\rVert\leq\frac{V_{l}(0)}{\sqrt{\underline{m}}}e^{-\beta
t}+\frac{\bar{d}}{\beta}\sqrt{\frac{\bar{m}}{\underline{m}}}(1-e^{-\beta t})$
(56)

where $V_{l}(0)$ denotes the path integral at time zero between $\xi_{1}(0)$
and $\xi_{0}(0)$.

Given a contracting system
$\dot{\mathbf{y}}=\mathbf{f}(\mathbf{y},\mathbf{u},t)$ with the convergence
rate $\beta$, we consider the perturbed dynamic

$\displaystyle\dot{\mathbf{y}}$
$\displaystyle=\mathbf{f}(\mathbf{y},\mathbf{u},t)+\mathbf{d}(\mathbf{y},\mathbf{u},t)$
(57) $\displaystyle\mathbf{u}$
$\displaystyle=T_{a}\bm{\pi}(\mathbf{s}_{1},\mathbf{s}_{2})=T_{a}\bm{\pi}(\mathbf{z},\smallint\mathbf{z}dt)$

where $\mathbf{d}$ represents unknown, bounded, and deterministic
disturbances.

Defining $\mathbf{\bar{y}}=[\mathbf{y}^{T},\mathbf{u}^{T}]^{T}$, the perturbed
dynamic can be represented as

$\mathbf{\dot{\bar{y}}}=\mathbf{\bar{f}}(\mathbf{\bar{y}},t)+\mathbf{\bar{d}}$
(58)

where

$\displaystyle\mathbf{\bar{f}}(\mathbf{\bar{y}},t)$ (59) $\displaystyle=$
$\displaystyle\begin{bmatrix}\mathbf{f}(\mathbf{y},\mathbf{u},t)\\\
\dot{T}_{a}\bm{\pi}+(T_{a}\frac{\partial\bm{\pi}}{\partial\mathbf{s}_{1}}\dot{T}_{y}+T_{a}\frac{\partial\bm{\pi}}{\partial\mathbf{s}_{2}}T_{y})\mathbf{y}+(T_{a}\frac{\partial\bm{\pi}}{\partial\mathbf{s}_{1}}T_{y})f(\mathbf{y},\mathbf{u},t)\end{bmatrix}$

and

$\mathbf{\bar{d}}=\begin{bmatrix}\mathbf{d}(\mathbf{y},\mathbf{u},t)\\\
(T_{a}\frac{\partial\bm{\pi}}{\partial\mathbf{s}_{1}}T_{y})\mathbf{d}(\mathbf{y},\mathbf{u},t)\end{bmatrix}$
(60)

Therefore, applying Theorem 2.4 in [30], the distance between the trajectories
of the contracting system and its perturbed dynamics is converging to a
bounded error ball with the center at the equilibrium in the following way,

$\displaystyle\lVert\xi_{1}(t)-\xi_{2}(t)\rVert\leq C_{1}e^{-\beta
t}+C_{2}\sup_{\mathbf{y},\mathbf{u},t}\lVert\Theta_{M}\mathbf{\bar{d}}\rVert\frac{1-e^{-\beta
t}}{\beta}$ (61)

where $\Theta_{M}$ is in Eq. (1), $C_{1}$ and $C_{2}$ are constants determined
by the bounds of $M=\Theta_{M}^{T}\Theta_{M}$ and the initial conditions
$\xi_{0}(0)$ and $\xi_{1}(0)$, provided with bounded $T_{a}$, $T_{y}$, and
$\partial\bm{\pi}/\partial\mathbf{s}_{1}$.

####  X-E2 Robustness to model errors

Let $\hat{\mathbf{f}}$ denote the model used to compute the constraints and
$\mathbf{f}$ the real dynamics. Applying the stability theorem in Section.
X-B, we can obtain the following result in robustness.

Given $\hat{T}_{y}$, $\hat{T}_{a}$ that are computed with $\hat{\mathbf{f}}$,
the system remains stable at the presence of model errors, if there exists
$\beta>0$ such that $\forall\mathbf{x}$, $\forall t>0$,

$(F_{1,real}^{T}+F_{1,real})\preceq(\hat{F}_{1}^{T}+\hat{F}_{1})\prec-[\beta+max(\nu^{+},0)]I$
(62)

where

$\displaystyle F_{1,real}$
$\displaystyle=\begin{bmatrix}\hat{T}_{y}\frac{\partial\mathbf{f}}{\partial\mathbf{y}}\hat{T}_{y}^{-1}&\hat{T}_{y}\frac{\partial\mathbf{f}}{\partial\mathbf{u}}\hat{T}_{a}\\\\[8.0pt]
\frac{\partial\bm{\pi}}{\partial\mathbf{s}_{1}}\hat{T}_{y}\frac{\partial\mathbf{f}}{\partial\mathbf{y}}\hat{T}_{y}^{-1}+\frac{\partial\bm{\pi}}{{\partial\mathbf{s}_{2}}}&\frac{\partial\bm{\pi}}{\partial\mathbf{s}_{1}}\hat{T}_{y}\frac{\partial\mathbf{f}}{\partial\mathbf{u}}\hat{T}_{a}\end{bmatrix}$
(63) $\displaystyle\hat{F}_{1}$
$\displaystyle=\begin{bmatrix}\hat{T}_{y}\frac{\partial\hat{\mathbf{f}}}{\partial\mathbf{y}}\hat{T}_{y}^{-1}&\hat{T}_{y}\frac{\partial\hat{\mathbf{f}}}{\partial\mathbf{u}}\hat{T}_{a}\\\\[8.0pt]
\frac{\partial\bm{\pi}}{\partial\mathbf{s}_{1}}\hat{T}_{y}\frac{\partial\hat{\mathbf{f}}}{\partial\mathbf{y}}\hat{T}_{y}^{-1}+\frac{\partial\bm{\pi}}{{\partial\mathbf{s}_{2}}}&\frac{\partial\bm{\pi}}{\partial\mathbf{s}_{1}}\hat{T}_{y}\frac{\partial\hat{\mathbf{f}}}{\partial\mathbf{u}}\hat{T}_{a}\end{bmatrix}$

and $\nu^{+}$ is the upper bound of the eigenvalues of $F_{2}^{T}+F_{2}$ in
Eq. (21).

###  X-F Implementation in the PPO example

The dynamics of the peg-robot is

$\tau_{x}\dot{x}+x=u_{x},\ \ \tau_{z}\dot{z}+z=u_{z}$ (64)

where $\tau_{x}=0.0437$, $\tau_{z}=0.01$, $x\in[1,3]$, and $z\in[-3,1]$. Given
the surface profile $g(x)=K_{1}\sin(x)+K_{2}\cos(x)$ with randomly sampled
coefficient $K_{1}\in[-0.01,0.01]$ and $K_{2}\in[-0.01,0.01]$, the contact
force can be approximated by

$f=K_{sur}\min{(z-g(x),\ 0)}$ (65)

where $K_{sur}$ represents the stiffness of the surface material, randomly
sampled from $K_{sur}\in[1,31]$. The desired $x_{d}$ and $f_{d}$ are also
randomly sampled at the beginning of each trajectory.

The implementation includes (a) computation of $T_{y}$, $T_{a}$,
$\Lambda_{ii}$, and $R_{ii}$ and (b) simplification into constraints on MLP
weights.

##### Computation of $T_{y}$, $T_{a}$, $\Lambda_{ii}$, and $R_{ii}$

When the peg is in contact with the surface, the dynamic of $x$ and $f$ can be
represented as

$\displaystyle\dot{x}=$
$\displaystyle-\frac{1}{\tau_{x}}x+\frac{1}{\tau_{x}}u_{x}$ (66)
$\displaystyle\dot{f}=$ $\displaystyle K_{sur}(\dot{z}-\frac{dg}{dx}\dot{x})$
$\displaystyle=$
$\displaystyle-\frac{1}{\tau_{z}}f+[\frac{K_{sur}}{\tau_{x}}\frac{dg}{dx}x-\frac{K_{sur}}{\tau_{z}}g(x)]$
$\displaystyle-\frac{K_{sur}}{\tau_{x}}\frac{dg}{dx}u_{x}+\frac{K_{sur}}{\tau_{z}}u_{z}$

Defining $\mathbf{x}=[x,f]^{T}$, the nonlinear system model
$\dot{\mathbf{x}}=\mathbf{f}(\mathbf{x},\mathbf{u},t)$ becomes

$\displaystyle\dot{\mathbf{x}}=\frac{d}{dt}\begin{bmatrix}x\\\
f\end{bmatrix}=$ $\displaystyle\begin{bmatrix}-\frac{1}{\tau_{x}}x\\\\[7.0pt]
[\frac{K_{sur}}{\tau_{x}}\frac{dg}{dx}x-\frac{K_{sur}}{\tau_{z}}g(x)]-\frac{1}{\tau_{z}}f\end{bmatrix}$
(67) $\displaystyle+\begin{bmatrix}\frac{1}{\tau_{x}}&0\\\\[7.0pt]
-\frac{K_{sur}}{\tau_{x}}\frac{dg}{dx}&\frac{K_{sur}}{\tau_{z}}\end{bmatrix}\mathbf{u}$

where $g(x)$ is the surface profile.

Instead of accurately computing the partial derivatives
$\partial\mathbf{f}/\partial\mathbf{x}$ and
$\partial\mathbf{f}/\partial\mathbf{u}$, we approximate the derivative with

$\displaystyle\frac{\partial\mathbf{f}}{\partial\mathbf{x}}$
$\displaystyle\approx\begin{bmatrix}-\frac{1}{\tau_{x}}&0\\\\[7.0pt]
\frac{K_{sur}}{\tau_{x}}\frac{dg}{dx}-\frac{K_{sur}}{\tau_{z}}\frac{g(x)}{x}&-\frac{1}{\tau_{z}}\end{bmatrix}$
(68) $\displaystyle\frac{\partial\mathbf{f}}{\partial\mathbf{u}}$
$\displaystyle=\begin{bmatrix}\frac{1}{\tau_{x}}&0\\\\[7.0pt]
-\frac{K_{sur}}{\tau_{x}}\frac{dg}{dx}&\frac{K_{sur}}{\tau_{z}}\end{bmatrix}$

We then compute $T_{y}$, $T_{a}$, $\Lambda_{ii}$, and $R_{ii}$ time-stepwise
according to Theorem 2.

##### Simplification into constraints on MLP weights

Because $\Lambda_{11}=-1/\tau_{x}$ and $\Lambda_{22}=-1/\tau_{z}$ are
negative, the constraints become that the signs of network Jacobians are
opposite to $R_{ii}$ (Eq. (16)).

We use a parallel sequence of 3-layered MLPs with hyperbolic tangent function
as the activation function and with positive weights. To satisfy Eq. (16), a
conditional multiplication is added: if $R_{ii}>0$, the output of the control
policy is multiplied by -1.

TABLE I: Hyperparameters in Peg Maze Parameter | TD3/C-TD3 | HIRO
---|---|---
Gaussian noise std | HL:1.0, LL:0.2 | HL:1.0, LL:1.0
gradient clipping | 0.1 | 0.1
target update period | 100 | 100
train steps | HL:100, LL:200 | HL:100, LL:200
TABLE II: Peg Push Task Settings Geometry | Length (y-axis) | Width (x-axis) | Depth (z-axis)
---|---|---|---
Box | 36 mm | 36 mm | 20 mm
Sliding cover | 12 mm | 36 mm | 4 mm
Slot | - | 4 mm | 7 mm
Obstacle | 2 mm | 2 mm | 2 mm
Target | 1 mm | 1 mm | 12.5* mm
Peg | - | radius=0.8 mm | 20 mm
TABLE III: Parameters of TD3-based Algorithms critic learning rate | 0.001
---|---
target learning rate | 0.0001
adam beta1 | 0.9
adam beta2 | 0.999
high-level dqda clipping | None
low-level dqda clipping | 0.1
target q clipping | None
target tau | 0.005
target update period | 50
action regularizer | None
observation normalization | None
reward scaling | 1
high-level discount factor | 0.5
low-level discount factor | 1
high-level exploration noise | 0.1$\mathcal{N}(0,0.1)$ mm
low-level exploration noise | $\mathcal{N}(0,0.1)$
batch size | 128
buffer size | 65536
samples to buffer per update | 1024
train step per update | 100
Translation Kp | 20
Translation Kd | 30
Orientation Kp | 30
Orientation Kd | 40
high-level critic network. | (400, 300) with tanh
low-level critic network. | (400, 300) with tanh
high-level actor network. | (400, 300) with tanh
low-level actor network. | (200, 200) with $W>0$ and tanh

###  X-G Peg Maze Example

The goal is sampled between $x\in(0.41,0.415)$, $y\in(-0.025,-0.02)$, and
$z\in(0.005,0.015)$ with the success defined as the L2 norm is less than
0.005. The robot starts almost from with the same position with small random
changes in joint position.

As for the operational space controller, the force-based control is applied,
that is, Section 3.3 in [15]. The task space control uses $K_{p,trans}=40$,
$K_{d,trans}=20$, , $K_{p,rot}=20$, and $K_{p,rot}=10$. The joint position
controller uses $K_{p}=1$. Because of the high-level action in joint space is
very small, we used $K_{d}=0$ and the low-level action from HIRO is forced to
have the same sign with the error for stability guarantee. For HIRO, the high-
level action in joint space is within $[-0.1,0.1]$ that is small enough to
avoid the end effector to move unstably. In TD3/C-TD3, we used the task space
directly as the the high-level action space, that is, $x\in(0.40,0.43)$,
$y\in(-0.03,-0.015)$, and $z\in(-0.005,0.015)$.

Hyperparameters are in Table I and others are the defaults from HIRO [17].
Negative norm of error is used as rewards.

###  X-H Peg-Push Example

The geometry is listed in Table. II and the depth of target refers to the
distance from the target to the top surface of the box. The initial peg
position is randomly sampled in a $8\times 8\times 8mm^{3}$ volume on the
right of the sliding cover. The contact parameters in Mujoco are solref =
‘0.01 1’ and solimp=‘0.9 0.95 0.001’. Control frequency is 50 Hz and the high-
level action frequency is 2.5 Hz (20 environment steps). The trajectory length
is set to 20 s.

The high-level reward is
$-tanh(50||\mathbf{e}_{trans}||)-0.1||\mathbf{e}_{rot}||$ where the error is
with respect to the goal (the target orientation is perpendicular to the
table); the low-level reward is the negative sum of the translation error norm
and the rotation error norm with respect to the high-level action.

As for the operational space control, the acceleration-based control is used,
that is, Section 3.2 in [15]. The high-level action provides the 20-step
increment within the range $\pm 5$ mm for translation and $\pm 0.5$ radian for
orientation. The low-level action is limited to $\pm 1$.

Hyperparameters and controller gains are listed in Table. III. Otherwise,
default parameters are used.