# Final state sensitivity in noisy chaotic scattering

Alexandre R. Nieto<EMAIL_ADDRESS>Nonlinear Dynamics, Chaos and
Complex Systems Group, Departamento de Física, Universidad Rey Juan Carlos,
Tulipán s/n, 28933 Móstoles, Madrid, Spain  Jesús M. Seoane  Nonlinear
Dynamics, Chaos and Complex Systems Group, Departamento de Física, Universidad
Rey Juan Carlos, Tulipán s/n, 28933 Móstoles, Madrid, Spain  Miguel A.F.
Sanjuán  Nonlinear Dynamics, Chaos and Complex Systems Group, Departamento de
Física, Universidad Rey Juan Carlos, Tulipán s/n, 28933 Móstoles, Madrid,
Spain  Department of Applied Informatics, Kaunas University of Technology,
Studentu 50-415, Kaunas LT-51368, Lithuania

###### Abstract

The unpredictability in chaotic scattering problems is a fundamental topic in
physics that has been studied either in purely conservative systems or in the
presence of weak perturbations. In many systems noise plays an important role
in the dynamical behavior and it models their internal irregularities or their
coupling with the environment. In these situations the unpredictability is
affected by both the chaotic dynamics and the stochastic fluctuations. In the
presence of noise two trajectories with the same initial condition can evolve
in different ways and converge to a different asymptotic behavior. For this
reason, even the exact knowledge of the initial conditions does not
necessarily lead to the predictability of the final state of the system.
Hence, the noise can be considered as an important source of unpredictability
that cannot be fully understood using the conventional methods of nonlinear
dynamics, such as the exit basins and the uncertainty exponent. By adopting a
probabilistic point of view, we develop the concepts of probability basin,
uncertainty basin and noise-sensitivity exponent, that allow us to carry out
both a quantitative and qualitative analysis of the unpredictability on noisy
chaotic scattering problems.

###### pacs:

05.45.Ac,05.45.Df,05.45.Pq

##  I Introduction

Unpredictability is a fundamental topic in nonlinear science, and it is a
consequence of the sensitive dependence to initial conditions, inherent to
chaotic systems. However, there are several ways to understand
unpredictability, and for each of them many concepts, methods and tools have
been developed. The unpredictability can be defined as the difficulty to
predict the evolution of the trajectories. With this point of view, the
expansion entropy [1], the Kolmogorov-Sinai entropy [2, 3], the topological
entropy [4], and other measures have been developed. Nevertheless, in certain
physical situations such as chaotic scattering problems [5], we are interested
in the asymptotic behavior rather than the evolution of the trajectories. In
this case, the unpredictability is understood as the difficulty to predict the
final state of a trajectory that starts with a particular initial condition.
The asymptotic behavior can be a fixed point or a chaotic attractor in
dissipative systems, a leak in chaotic Hamiltonian maps or an opening in the
potential in open Hamiltonian systems. The most common tools to analyze the
unpredictability in this kind of systems are the basins of attraction [6] and
the exit basins [7]. A basin of attraction of a dissipative system is defined
as the set of initial conditions that are attracted to a certain attractor. On
the other hand, an exit basin of a Hamiltonian system is the set of initial
conditions that escape through a particular exit of the system. When two
different exits (or attractors) coexist, two basins appear, separated by a
smooth or a fractal boundary. The knowledge of the structure of the basin
boundary, together with other characteristics of the basins, allows to
understand and quantify the unpredictability of the system. For this purpose,
the uncertainty exponent [8] and the basin entropy [9] are some of the most
powerful tools when working on fundamental aspects of physical systems. On the
other hand, in the case of applied physics and some fields of engineering,
many methods to characterize the dynamical integrity [10] of the basins of
attraction have been developed. Some examples are the anisometric local
integrity measure (ALIM) [11], the integrity factor (IF) [12], and the local
integrity measure (LIM) [13].

Recently, much work has been done on the unpredictability of chaotic
scattering. In particular, current research has shifted the focus from purely
conservative systems to the effects of relativistic corrections [14], small
perturbations as dissipation [15], and periodic forcing [16, 17]. However, the
noise is ubiquitous in nature, so that important lines of research have arisen
to elucidate the noise effects in chaotic scattering [18, 19, 20, 21, 22].
Accordingly, in this manuscript we carry on a research on noisy chaotic
scattering, mainly focusing our attention on developing numerical and
theoretical techniques that allow to characterize the final state sensitivity.

Our research arises from the observation that some of the tools and methods
used to study unpredictability in chaotic scattering are useless in the
presence of stochastic fluctuations. As good examples, we can mention the exit
basins and the uncertainty exponent. These tools are based on the high
sensitivity of the asymptotic behavior to initial conditions or parameters.
However, if we consider a small amount of noise the system is not
deterministic anymore and the situation changes drastically. In this scenario,
two exactly identical initial conditions can evolve in a different way and
converge to a different asymptotic state. Hence, in presence of noise we are
not interested anymore in the neighborhood of the initial condition but in the
initial condition itself.

Next, we develop some tools and numerical methods that can help to understand
and quantify the effects of noise on the unpredictability of chaotic
scattering systems. With this goal in mind, and without loss of generality, we
have used two paradigmatic open Hamiltonian systems. Nevertheless, we expect
that the main results could be general in chaotic scattering problems in which
the noise models the effect of internal irregularities or the coupling of the
system with the environment. We presume potential applications to several
fields of physics such as celestial mechanics [14, 23] and plasma physics [24,
25], among others. In particular, both the exit basins and the noise play an
important role in the magnetic behavior in tokamaks [26, 27, 28]. Other
possible applications include chaotic scattering problems in different fields
of science, such as medicine [29], biology [30] or chemistry [31, 32].

The structure of this paper is as follows. In Sec. II, we describe the models
of this work, the Hénon-Heiles system and the four-hill potential, both in
presence of a source of additive uncorrelated Gaussian noise. The description
of the concepts and methods to compute the probability basin and the
uncertainty basin is carried out in Sec. III and Sec. IV, respectively. In the
latter, we also focus our attention on the relation between the structure of
the uncertainty basin and the stable manifold of the chaotic saddle. In Sec.
V, we develop the concept of noise-sensitivity exponent as a tool to quantify
the final state sensitivity to noise. Finally, in Sec. VI, we present the main
conclusions of this work.

##  II Description of the models

One of the models that we have used to illustrate the ideas and methods of
this manuscript is the Hénon-Heiles system [33]. This Hamiltonian system
appeared for the first time in $1964$ as a model of a galactic potential. The
Hamiltonian is characterized by a nonlinear axisymmetric potential and is
given by

${\cal{H}}=\frac{1}{2}(\dot{x}^{2}+\dot{y}^{2})+\frac{1}{2}(x^{2}+y^{2})+x^{2}y-\frac{1}{3}y^{3},$
(1)

where $x$ and $y$ denote the space coordinates, and $\dot{x}$ and $\dot{y}$
are the the momentum coordinates.

The system has been extensively studied due to the fractal exit basins and the
wide variety of dynamical behaviors existing for different values of the
energy [34, 35]. To visualize the system in the context of the
unpredictability analysis, in Fig. 1 we show two exit basins for different
values of the energy $E=0.25$ and $E=0.45$. When the energy of the system is
increased, the fractal basin boundaries become thinner and less fractal, with
a consequent reduction on the unpredictability. This decreasing can be
quantified by using both the uncertainty exponent and the basin entropy.

Figure 1: Exit basins of the Hénon-Heiles system in physical space with energy
(a) $E=0.25$ and (b) $E=0.45$. The colors red, green and blue refer to initial
conditions leading to the three exits of the potential: Exit 1 ($y\to\infty$),
Exit 2 ($x,y\to-\infty$) and Exit 3 ($x\to\infty,y\to-\infty$).

For the purposes of this research, we have included in the equations of motion
a source of additive uncorrelated Gaussian noise, so the equations of motion
become [22]

$\displaystyle\dot{p}$ $\displaystyle=-x-2xy+\sqrt{2\xi}\eta_{x}(t)$ (2)
$\displaystyle\dot{q}$ $\displaystyle=-y-x^{2}+y^{2}+\sqrt{2\xi}\eta_{y}(t),$

where $\xi$ is the intensity of the noise and $\eta_{x}(t)$, $\eta_{y}(t)$ are
white noise processes with variance $\sigma^{2}=2\xi$ and mean $\mu=0$ .

For comparison purposes, we have also used the four-hill system [36, 37],
whose potential consists of four hills located at $(x,y)=(\pm 1,\pm 1)$ and
its Hamiltonian is given by:

${\cal{H}}=\frac{1}{2}(\dot{x}^{2}+\dot{y}^{2})+x^{2}y^{2}e^{-(x^{2}+y^{2})}.$
(3)

Under the same considerations that we have mentioned in the case of the Hénon-
Heiles, the equations of motion in presence of noise read

$\displaystyle\dot{p}$
$\displaystyle=2xy^{2}(x^{2}-1)^{-(x^{2}+y^{2})}+\sqrt{2\xi}\eta_{x}(t)$ (4)
$\displaystyle\dot{q}$
$\displaystyle=2x^{2}y(y^{2}-1)^{-(x^{2}+y^{2})}+\sqrt{2\xi}\eta_{y}(t).$

To visualize the system, we represent the exit basins for different energies
$E=0.01$ and $E=0.1$ in Fig. 2. In the low energy case, the basin boundary
occupies an important part of the phase space. While in the high energy case
the basin boundary is thinner, letting space to extensive safe regions of high
predictability. If we increase even more the energy until $E=1/e^{2}$, the
trajectories evolve on the top of the hills and the scattering becomes non-
chaotic.

Figure 2: Exit basins of the four-hill system in the physical space with
energy (a) $E=0.01$ and (b) $E=0.10$. The colors green, yellow, blue and red
refer to initial conditions leading to the four exits of the potential: Exit 1
($y\to\infty$), Exit 2 ($x\to-\infty$), Exit 3 ($y\to-\infty$) and Exit 4
($x\to\infty$).

We have solved the stochastic differential equations for the Hénon-Heiles and
the four-hill system by using the stochastic Heun method [38], as was
previously used in Refs. [22, 39]. In both systems, the numerical methods gave
stable and convergent solutions. To ensure the effectiveness of the method, we
have also compared the results by using the numerical schemes of the Euler-
Maruyama method and the stochastic fourth-order Runge-Kutta method.

##  III Probability basins

Perhaps the most obvious and at the same time interesting effect of noise on
the escape dynamics is that identical initial conditions can escape through
different exits and spend a different time in the scattering region. To
illustrate this, we represent three trajectories with the same initial
condition escaping through different exits of the potential of the Hénon-
Heiles system in Fig. 3(a). In addition, the escape times obtained for
different launchings of the same initial condition are represented in Fig.
3(b). This figure reminds us the very usual scattering function. However, here
no coordinate is varied and the apparently disordered escape times are due to
the effects of the noise. Even with a perfectly known initial condition, we
cannot say by which exit the particle will escape and we can only establish
that the trajectory will remain in the scattering region during a certain
range of time $T\in[19,116]$.

Figure 3: Escape dynamics of the noisy Hénon-Heiles system with $E=0.25$ and
$\xi=10^{-5}$, showing (a) three trajectories in the physical space launched
from the same initial conditions escaping through different exits, and (b)
escape time distribution of $100$ launchings with the same initial condition.

In this situation, where the main characteristics of the scattering process
(asymptotic behavior and escape time) may change in different simulations, it
is necessary to adopt a probabilistic approach in order to understand the
escape dynamics. Hence, the escape of a trajectory should be understood as a
probability of escape. Additionally, the escape times of an initial condition
should be analyzed in terms of the average and the variance of the individual
escape times. Regarding unpredictability, one of the most useful portraits is
the structure of the exit basins. However, in presence of noise this
representation lacks meaning since for a given initial condition, the exit is
not uniquely defined. So, in noisy systems we do not have a unique exit basin
representation, we might have a different one for each different simulation we
carry out. As it has been reported in Refs. [22, 39], the exit basins appear
smeared or blurred, and the basin boundary becomes shaded off due to the
effects of noise. Hence, a deeper understanding of the system is not possible
by using the exit basins. Nevertheless, we can construct a similar
representation that we call probability basin, in which every initial
condition in the grid does not have an associated asymptotic behavior, but a
probability of escaping through a particular exit. Consequently, we will have
as many probability basins as there are exits in the system. For the
computation of the probability basins, we start fixing a value of the noise
intensity and establishing a grid of $P$ initial conditions in the phase space
region that we want to analyze. Then, we compute $N$ times the trajectory of
the very same initial condition and we detect the exit through which they
escape. Once we know the exits of all $N$ trajectories, we can establish the
probability of each exit. By repeating this procedure for all $P$ initial
conditions we will obtain the probability basins.

To illustrate the above, in Fig. 4 we use a color code map to represent the
probability basin of the Exit 1 of the Hénon-Heiles system with energy
$E=0.25$ and different values of the noise intensity. We do not show the
probability basins of the Exits $2$ and $3$ because, due to the symmetry of
the basin, they are a simply $2\pi/3$ rotation of the probability basin of the
Exit 1. In panel $(a)$ and $(c)$ the almost negligible noise intensity
$\xi=10^{-10}$ generates mainly extreme values in the probability, indicating
an almost deterministic behavior. However, even a very small amount of noise
of intensity $\xi=10^{-6}$ can generate an important source of uncertainty in
the vicinity of the deterministic basin boundary. The highly predictable
regions become uncertain and the probabilistic nature of the escapes manifests
along an important part of the physical space. Hence, we can conclude that the
effects of noise make even more difficult the prediction of the final state of
the system.

Figure 4: Color map showing the probability basin for the Exit 1 of the Hénon-
Heiles system with energy $E=0.25$, and different noise intensities (a)
$\xi=10^{-10}$ and (b) $\xi=10^{-6}$. Panels (c) and (d) are magnifications of
the probability basin of panels (a) and (b), respectively. The yellow color
indicates maximum probability of finding the Exit 1, while dark blue implies
zero probability. Intermediate colors refer to intermediate probabilities, as
shown in the color bar. To generate these figures, we have computed $100$
times the exit of each initial condition in a $400\times 400$ grid and we have
represented the probability of the Exit 1.

As we increase the noise intensity, the uncertainty and the blurring of the
probability basins also increase. To illustrate this, we represent in Fig. 5
the evolution of the segment $x=0$ of the probability basin of the Exit 1 when
the noise intensity varies from $10^{-10}$ to $10^{-3}$. We have not
considered higher noise intensities because they imply the dominance of the
stochastic behavior and the trajectories are simply random walks. In the
previous figure, we can observe how the regions where particles surely escape
through Exit 1 decrease in size as noise increases. For $\xi=10^{-3}$ the only
region where this happens is the one near to $y=0$, which is formed by initial
conditions above the Lyapunov orbits [40] that do not take place in the
chaotic scattering process.

Figure 5: Evolution of the segment $x=0$ of the probability basin of the Exit
1 of the Hénon-Heiles system with increasing noise. The energy of the system
is $E=0.25$, while the noise intensity is varied from $\xi=10^{-10}$ to
$\xi=10^{-3}$. The color code is as described in the caption of Fig. 4. The
base of the log is $10$.

##  IV Uncertainty basin

The probability basin representation is a useful tool to visualize the
probability distribution of the exits along the phase space. Nevertheless, in
many contexts we are not specifically interested in the probability, but on
the unpredictability of the escape. For this reason, in deterministic chaotic
scattering the boundary of the exit basins plays a central role when analyzing
and quantifying the unpredictability of a system. It was not in vain that a
large number of tools and concepts have been developed just to understand the
basin boundary and its sometimes rich fractal structure. However, in the case
of noisy chaotic scattering the basin boundary is not well defined since the
set of points that define it changes for different simulations. For this
reason it is necessary to change again our point of view and understand the
uncertain initial conditions as the set of initial conditions that can change
for different launchings. We call this set uncertainty basin, $U_{B}$. We
recall again that no uncertainty in the initial condition is considered here
and the source of uncertainty is all due to the influence of the noise on the
chaotic dynamics. In order to construct the uncertainty basin, we follow a
similar method as in the case of the probability basins. First, we define a
grid of $P$ initial conditions in all the phase space and we launch each of
them $N$ times, labeling an initial condition as uncertain if the exit through
which it escapes changes at least in one simulation. To illustrate the result,
we simply associate a different color to the certain and the uncertain initial
conditions. In Fig. 6 we show in yellow the uncertainty basin of the Hénon-
Heiles with energy $E=0.25$ and different noise intensities (a)
$\xi=10^{-10}$, (b) $\xi=10^{-8}$, (c) $\xi=10^{-6}$ and (d) $\xi=10^{-5}$. It
can be observed at naked eye that the uncertainty basin is located in the
vicinity of the stable manifold of the chaotic saddle ($W_{s}$) (i.e. the
boundary of the exit basins) and grows with increasing noise. This means that
some initial conditions that were predictable in the deterministic system due
to the large distance to the basin boundary become uncertain due to the
effects of noise. If we consider an initial condition sufficiently close to
the stable manifold of the chaotic saddle, a small perturbation can move it
from one basin to another, so it is intuitive that the uncertainty basin lies
close to $W_{s}$. Moreover, we can establish that

$\lim_{\xi\to 0}U_{B}=W_{s},$ (5)

which constitutes an accurate method to determine the structure of $W_{s}$.

Figure 6: Uncertainty basins of the Hénon-Heiles system with energy $E=0.25$,
and different noise intensities (a) $\xi=10^{-5}$, (b) $\xi=10^{-6}$, (c)
$\xi=10^{-8}$ and (d) $\xi=10^{-10}$. To generate these figures we have
computed the exit of each initial condition $100$ times and labeled the
initial condition as certain (black) or uncertain (yellow).

The uncertainty basin is a faithful portrait of the unpredictability of the
system in the presence of a particular noise intensity. By simply computing
and exploring the uncertainty basin, we can conclude which regions are
uncertain and which ones remain predictable under the effects of noise.

Even if the uncertainty basin is a set that emerges from the fine structure of
the stable manifold of the chaotic saddle, as we increase the noise its once
rich and complex structure becomes blurred, with the consequently lost of
fractality. All those “threads” that make up the boundary widen in the
uncertainty basin and end up coming together. As a consequence, they define
extensive regions of smooth geometry. This result is shown in Fig. 7, where
several magnifications on the exit basin and the uncertainty basin are
represented. Meanwhile, the geometry of the basin boundary of the exit basin
is not simplified when magnified. As a matter of fact, the smooth nature of
the uncertainty basin manifest on a finite scale. One of the most significant
implications of the absence of fractality in the uncertainty basin is that its
size is independent of its resolution, which allows us to establish an
absolute fraction of uncertain initial conditions.

Figure 7: (a,b,c) Magnifications in the exit basins of the Hénon-Heiles in the
physical space with energy $E=0.25$ (See Fig. 1) and (d,e,f) same regions but
in the uncertainty basin with noise intensity $\xi=10^{-6}$ (See Fig. 6(c)).
It is clear that the uncertainty basin grows from the boundary of the exit
basin, connecting different regions and generating a smooth geometry.

##  V Quantifying the final state sensitivity

In the previous sections we have introduced the concepts of probability basin
and uncertainty basin, in order to give a graphic and qualitative
understanding of the unpredictability on noisy chaotic scattering. However, in
many physical problems we need a quantitative measure of the unpredictability
in order to analyze the effect of a parameter or to compare different systems.
In particular, even if two systems have the same fraction of uncertain initial
conditions for a particular value of the noise, we do not have any guarantee
that both have the same sensitivity to the noise, i.e., the fraction of
uncertain initial conditions will increase in the same amount when the noise
intensity increases. From the perspective of the prediction of the system, it
is preferable that the uncertainty grows slowly with increasing noise. In this
section, we develop a method to quantify the final state sensitivity to noise.

In deterministic systems one of the most common and useful measures is the
uncertainty exponent, that can be obtained by using the uncertainty algorithm
[8]. The method consists of launching several initial conditions under some
uncertainty $\delta$. Then, we say that an initial condition is uncertain
under an uncertainty $\delta$ if the asymptotic behavior is affected by the
error. The fraction of uncertain initial conditions is then related to
$\delta$ through a power law

$f_{u}(\delta)\sim\delta^{\alpha},$ (6)

where $\alpha$ is the uncertainty exponent.

In noisy chaotic scattering system we have two sources of uncertainty: the
error in the initial condition and the noise. For this reason, the fraction of
uncertain initial conditions, $f(\delta,\xi)$, does not satisfy $f(0,\xi)=0$
since some initial conditions are uncertain even without considering an error.
These initial conditions are those that change the exit for different
launchings. These initial conditions will be uncertain $\forall\delta$, so
they do not increase or decrease when changing $\delta$. This makes Eq. (6)
invalid in noisy systems. In fact, since the uncertainty exponent is defined
in the limit $\delta\to 0$, the fraction of uncertain initial conditions for a
fixed and not negligible noise intensity will always be
$f(\delta,\xi)=f(\xi)=k$, where $k$ is a constant. This implies that
$\alpha=0$ in noisy chaotic scattering systems.

Here, we are interested in the effects of noise, so we consider $\delta=0$ and
we set that an initial condition is uncertain if it can escape through
different exits for different launchings. Under this assumption the fraction
of uncertain initial conditions obeys the law

$f_{u}(\xi)\sim\xi^{\beta},$ (7)

where $\beta>0$ is a magnitude that characterizes the intrinsic uncertainty of
the noisy system, and we call it noise-sensitivity exponent.

As $\beta$ approaches to zero a reduction in the noise intensity has less
effect in decreasing the fraction of uncertain initial conditions, so we say
that the system has high sensitivity to noise. In the limit $\beta=0$ a
reduction in the noise intensity does not decrease the fraction of uncertain
initial conditions at all. On the other hand, for high values of $\beta$ the
fraction of uncertain initial conditions tends to zero for low noise
intensities and very large noise intensities are needed in order to increase
the unpredictability. Hence, we say that the system has low sensitivity to the
noise. We have computed the parameter $\beta$ for several open Hamiltonian
systems and their result ranges from $0.05$ to $0.5$. Furthermore, we have
observed that extreme values appear only where the basin boundary occupies
almost all the phase space ($\beta<0.05$) or in basins that have a unique
smooth boundary ($\beta>0.5$).

In order to obtain $\beta$, we begin computing $f_{u}(\xi)$ for different
values of $\xi$. In particular, for a fixed value of the noise intensity we
launch $P$ initial conditions $N$ times, and we compute the fraction of
uncertain initial conditions. We repeat this procedure for several values of
$\xi$ and we proceed to a $\log-\log$ representation that will satisfy

$\log_{10}{f_{u}(\xi)}=\beta\log_{10}{\xi}+\mbox{C},$ (8)

where $C$ is a constant and the slope of the straight line is $\beta$.

The quality of the linear fit turns out very well for all the simulations that
we have carried out (linear correlation coefficient $r>0.999$). As an example,
we show the resulting linear fit for the case of the Hénon-Heiles system with
energy $E=0.25$ in Fig. 8. The noise-sensitivity exponent is estimated to be
$\beta=0.109$. If we increase the energy to $E=0.45$ the result is
$\beta=0.2309$, so that the final state sensitivity due to the noise is higher
for $E=0.25$. Certainly, this case can be easily understood due to the big
difference in the size of the basin boundaries (see Fig. 1).

Figure 8: Logarithm of the fraction of uncertain initial conditions in
function of the logarithm of the noise intensity of the Hénon-Heiles system
with energy $E=0.25$. The parameter $\beta$ is estimated to be $\beta=0.109$.
To generate this figure, we have used $40$ different values of the noise
intensity. For every noise intensity the fraction of uncertain initial
conditions has been computed after launching $100000$ initial conditions
($100$ times each one)

Another interesting application of the noise-sensitivity exponent is to
compare the sensitivity to the noise for different systems. To illustrate
this, we have chosen the four-hill system with energies $E=0.01$ and $E=0.1$.
The exit basins of the conservative system for these values of the energy have
been shown in Fig. 2. In these cases the noise-sensitivity exponent is
estimated to be $\beta=0.078$ ($E=0.01$) and $\beta=0.2902$ ($E=0.1$). Hence,
for the lower energy value the system is more sensitive to the noise than the
Hénon-Heiles for $E=0.25$, while for the higher energy value the four-hill
system is less sensitive than the Hénon-Heiles with $E=0.45$.

##  VI Conclusions

In summary, our research reveals that the noise, even if the intensity is very
weak, has important implications on the predictability of chaotic scattering
systems. We have shown that in order to have a deeper understanding of the
effects of noise, some concepts and tools of nonlinear dynamics and chaos
should be revisited. In particular, the usual ways to represent the exit
basins and to compute the uncertainty exponent are meaningless in the presence
of noise. Throughout this work we have adopted a probabilistic point of view
and developed the concepts and methods to compute the probability basin,
uncertainty basin and noise-sensitivity exponent. Under all these concepts
underlies the idea that an initial condition is uncertain if it can change its
exit for different launchings. The probability and uncertainty basins allow a
qualitative analysis of the probability distribution of the escapes and the
structure of the uncertain regions, respectively. On the other hand, the
noise-sensitivity exponent gives a measure to quantify the final state
sensitivity in presence of noise. This tool can be useful to analyze how the
sensitivity to the noise manifests itself for different dynamical systems, as
well as to see its behavior when a parameter of the system varies.

We expect that the results shown in this work could be helpful in providing
new concepts and tools to future research concerning chaotic scattering
problems.

## Acknowledgments

We dedicate this paper to the memory of Prof. Tito Arecchi with whom we had
the pleasure of collaborating in different nonlinear dynamics projects. This
work has been financially supported by the Spanish State Research Agency (AEI)
and the European Regional Development Fund (ERDF) under Project No.
PID2019-105554GB-I00.

## References

  * [1] Hunt BR, Ott E. Defining chaos. Chaos 2015;25:097618.
  * [2] Kolmogorov AN. New metric invariant of transitive dynamical systems and endomorphisms of Lebesgue spaces. Dokl Akad Nauk SSSR 1959;119:861-64.
  * [3] Sinai YG. On the notion of entropy of a dynamical system. Dokl Akad Nauk SSSR 1959;124:768-71.
  * [4] Adler RL, Konheim AG, McAndrew MH. Topological entropy. Trans Amer Math Soc 1965;114:309-19.
  * [5] Seoane JM, Sanjuán MAF. New developments in classical chaotic scattering. Rep Prog Phys 2013;76:016001.
  * [6] Nusse HE, Yorke JA. Basins of attraction. Science 1996;271:1376-80.
  * [7] Contopoulos G. Order and Chaos in Dynamical Astronomy. Berlin: Springer; 2002.
  * [8] Grebogi C, McDonald SW, Ott E, Yorke JA. Final state sensitivity: An obstruction to predictability. Phys Lett A 1983;99:415-8.
  * [9] Daza A, Wagemakers A, Georgeot B, Guéry-Odelin D, Sanjuán MAF. Basin entropy: a new tool to analyze uncertainty in dynamical systems. Sci Rep 2016;6:31416.
  * [10] Thompson JMT. Chaotic phenomena triggering the escape from a potential well. Proc R Soc Lond A 1989;421:195-225.
  * [11] Belardinelli P, Lenci S, Rega, G. Seamless variation of isometric and anisometric dynamical integrity measures in basin’s erosion. Commun Nonlinear Sci Numer Simul 2018;56:499-507.
  * [12] Rega G, Lenci S. Identifying, evaluating, and controlling dynamical integrity measures in non-linear mechanical oscillators. Nonlinear Anal Theory Methods Appl 2005;63;902-14.
  * [13] Soliman MS, Thompson JMT. Integrity measures quantifying the erosion of smooth and fractal basins of attraction. J Sound Vib 1989;135:453-75.
  * [14] Bernal JD, Seoane JM, Vallejo JC, Huang L, Sanjuán MAF. Influence of the gravitational radius on asymptotic behavior of the relativistic Sitnikov problem. Phys Rev E 2020;102:042204.
  * [15] Seoane JM, Sanjuán MAF, Lai YC. Fractal dimension in dissipative chaotic scattering. Phys Rev E 2007;76:016208.
  * [16] Blesa F, Seoane JM, Barrio R, Sanjuán MAF. Effects of periodic forcing in chaotic scattering. Phys Rev E 2014;89:042909.
  * [17] Nieto AR, Seoane JM, Alvarellos JE, Sanjuán MAF. Resonant behavior and unpredictability in forced chaotic scattering. Phys Rev E 2018;98:062206.
  * [18] Rodrigues CS, Moura APS, Grebogi C. Random fluctuation leads to forbidden escape of particles. Phys Rev E 2010;82:026211.
  * [19] Silva RS, Manchein C, Beims MW. Exploring conservative islands using correlated and uncorrelated noise. Phys Rev E 2018;97:022219.
  * [20] Mills P. The influence of noise on a classical chaotic scatterer. Commun Nonlinear Sci Numer Simul 2006;11:899-906.
  * [21] Altmann EG, Englert A. Noise-enhanced trapping in chaotic scattering. Phys Rev Lett 2010;105:244102.
  * [22] Seoane JM, Sanjuán MAF. Exponential decay and scaling laws in noisy chaotic scattering. Phys Lett A 2008;372:110-6.
  * [23] de Assis SC, Terra MO. Escape dynamics and fractal basin boundaries in the planar Earth–Moon system. Celest Mech Dyn Astr 2014;120:105-30.
  * [24] Mathias AC, Viana RL, Kroetz T, Caldas IL. Fractal structures in the chaotic motion of charged particles in a magnetized plasma under the influence of drift waves. Phyica A 2017;469:681-94.
  * [25] Mathias AC, Kroetz T, Caldas IL, Viana RL. Chaotic magnetic field lines and fractal structures in a tokamak with magnetic limiter. Chaos Solitons Fractals 2017;104:588-98.
  * [26] Viana RL, Da Silva EC, Kroetz T, Caldas IL, Roberto M, Sanjuán MAF. Fractal structures in nonlinear plasma physics. Phil Trans R Soc A 2011;369:371-95.
  * [27] Turri G, Coda S, Moret JM, Martin Y, Sauter O. The effect of MHD noise on the vertical observer in tokamaks. Plasma Phys Control Fusion 2008;50:035012.
  * [28] Punjabi A, Ali H. An accurate symplectic calculation of the inboard magnetic footprint from statistical topological noise and field errors in the DIII-D. Phys Plasmas 2011;18:022509.
  * [29] Schelin AB, Károlyi G, de Moura APS, Booth NA, Grebogi C. Chaotic advection in blood flow. Phys Rev E 2009;80:016213.
  * [30] Scheuring I, Czárán T, Szabó P, Károlyi G, Toroczkai Z. Spatial Models of Prebiotic Evolution: Soup Before Pizza? Orig Life Evol Biosph 2003;33:319-55.
  * [31] Ezra GS, Waalkens H, Wiggins S. Microcanonical rates, gap times, and phase space dividing surfaces. J Chem Phys 2009;130:164118.
  * [32] Kawai S, Bandrauk AD, Jaffé C, Bartsch T, Palacián J, Uzer T. Transition state theory for laser-driven reactions. J Chem Phys 2007;126:164306.
  * [33] Hénon M, Heiles C. The applicability of the third integral of motion: Some numerical experiments. Astron J 1964;69:73-9.
  * [34] Aguirre J, Vallejo JC, Sanjuán MAF. Wada basins and chaotic invariant sets in the Hénon-Heiles system. Phys Rev E 2001;64:066208.
  * [35] Barrio R, Blesa F, Serrano S. Fractal structures in the Hénon-Heiles Hamiltonian. Europhys Lett 2008;82:10003.
  * [36] Bleher S, Grebogi C, Ott E. Bifurcation to chaotic scattering. Physica D 1990;46:87-121.
  * [37] Zotos EE. Elucidating the escape dynamics of the four hill potential. Nonlinear Dyn 2017;89:135-51.
  * [38] Kloeden PE, Platen E, Schurz H. Numerical Solution of Stochastic Differential Equations. Berlin: Springer; 1994.
  * [39] Bernal JD, Seoane JM, Sanjuán MAF. Weakly noisy chaotic scattering. Phys Rev E 2013;88:032914.
  * [40] Contopoulos G. Asymptotic curves and escapes in Hamiltonian systems. Astron Astrophys 1990;231:41-55.


# Learning Agile Locomotion and Adaptive Behaviors via RL-augmented MPC

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

In the context of legged robots, adaptive behavior involves adaptive balancing
and adaptive swing foot reflection. While adaptive balancing counteracts
perturbations to the robot, adaptive swing foot reflection helps the robot to
navigate intricate terrains without foot entrapment. In this paper, we manage
to bring both aspects of adaptive behavior to quadruped locomotion by
combining RL and MPC while improving the robustness and agility of blind
legged locomotion. This integration leverages MPC’s strength in predictive
capabilities and RL’s adeptness in drawing from past experiences. Unlike
traditional locomotion controls that separate stance foot control and swing
foot trajectory, our innovative approach unifies them, addressing their lack
of synchronization. At the heart of our contribution is the synthesis of
stance foot control with swing foot reflection, improving agility and
robustness in locomotion with adaptive behavior. A hallmark of our approach is
robust blind stair climbing through swing foot reflection. Moreover, we
intentionally designed the learning module as a general plugin for different
robot platforms. We trained the policy and implemented our approach on the
Unitree A1 robot, achieving impressive results: a peak turn rate of 8.5 rad/s,
a peak running speed of 3 m/s, and steering at a speed of 2.5 m/s. Remarkably,
this framework also allows the robot to maintain stable locomotion while
bearing an unexpected load of 10 kg, or 83% of its body mass. We further
demonstrate the generalizability and robustness of the same policy where it
realizes zero-shot transfer to different robot platforms like Go1 and AlienGo
robots for load carrying. Code is made available for the use of the research
community at https://github.com/DRCL-USC/RL_augmented_MPC.git

##  I Introduction

In the quest for the practical deployment of quadruped robots in real-world
scenarios, the integration of adaptive behavior into their motion remains a
challenge. This adaptive behavior consists of two essential dimensions: 1)
real-time adaptation to external perturbations, and 2) self-adjustments such
as foot reflection when a robot’s foot gets stuck in an obstacle. Current
advancements in legged mobility predominantly lean on Model Predictive Control
(MPC) and Reinforcement Learning (RL). MPC, which employs real-time
optimization over a set horizon to compute the optimal control sequence, often
requires substantial computational resources and careful parameter tuning. RL,
while notable for exceptional adeptness at navigating uneven terrains and
unexpected disturbances, demands extensive offline computation, and careful
reward tuning, and often produces policies tailored to specific robots. To
combine the benefits of both MPC and RL, we present an innovative approach
synthesizing the strengths of model-based control and reinforcement learning.
Our central objective is to bolster agility, robustness, and adaptive behavior
in blind locomotion through the integration of stance foot control and swing
foot reflection using RL.

Figure 1: Experiment result highlights. a) High-speed steering in place; b)
High-speed running; c) High-speed running and steering; d) Generalization of
the same policy across different robot platforms; e) Blind stair climbing with
swing foot reflection.

MPC, as validated by multiple studies [1, 2, 3, 4, 5, 6] has gained traction
in the legged robot community for its capability to handle the hybrid nature
of quadrupedal locomotion under constraints. Central to the MPC paradigm is
its prediction based on simplified dynamics, offering a future state
estimation while preserving real-time computational feasibility. Yet, these
simplified models inherently come with model uncertainties. For instance, the
single rigid dynamics (SRB) model overlooks the resultant dynamics from leg
momentum and external disturbance. Further complexities arise when translating
these optimized trajectories into joint-level commands. Existing methodologies
adopt a hierarchical control structure for this conversion, using techniques
like Jacobian mapping [7], control barrier functions[8], and joint-level whole
body control[9, 2]. The problem of addressing uncertainties in legged motion
has been approached with adaptive control methods [10, 11, 12, 13, 14, 15,
16], adjusting the control parameters online. In addition, an implicit
limitation of the MPC framework is its inherent decoupling of stance foot
control from swing foot control due to the intricate modeling challenges of
their interplay.

Reinforcement learning offers a tantalizing alternative [17, 18, 19, 20, 21,
22, 23, 24], implicitly deciphering the dynamics interplay between the stance
and swing control for all kinds of locomotion. In this paradigm, agents
continually engage with environments, iteratively refining their action
strategies based on the reward, resulting in the mastery of complex terrains
and adaptive behavior attuned to environmental dynamics. Nevertheless,
deploying RL in real-world scenarios raises legitimate concerns about its
generalizability and safety. The aforementioned challenges underscore the
urgency for a control framework evolution, one that concurrently addresses
model uncertainties and optimizes swing foot reflection with regularized
motions.

Researchers integrate MPC and reinforcement learning to combine the benefits
of RL and model-based control. In [25], the RL framework utilizes model-based
optimal control to generate reference motion and then leverages motion
imitation technique[26] to learn versatile legged motion. RL [27, 28] is also
utilized to learn parameters or dynamics in the MPC problem. [29] proposes an
RL-based control of the accelerations of an SRB model which allows robust sim-
to-real transfer. [30] leverages RL to infer the set of unmodeled dynamics for
the RMPC framework for adaptive locomotion. Additionally, [31] proposed an
online supervised learning technique to derive a residue model to address the
model uncertainties for model-based controller. What sets our research apart
from these studies is our innovative synthesis of stance foot control and
swing foot reflection by leveraging RL, enabling adaptive balancing and
adaptive foot reflection within one unified framework.

In this paper, we present a novel RL-augmented MPC framework tailored to
enhance blind locomotion for quadruped robots. By leveraging RL, we synthesize
stance foot control and swing foot reflection from the linear MPC framework
[1], specifically addressing the inherent issues of model uncertainty and the
pre-defined swing foot trajectory. By tackling the dual challenges of adapting
to model uncertainty and optimizing foot reflection, we successfully
demonstrated improved agility, robustness, and adaptive behavior in blind
legged locomotion. Notably, our research introduces robot-agnostic action and
observation spaces to guarantee the policy’s generalizability across various
robot platforms. Our proposed framework has the following contributions:

  * •

We introduce a novel RL-augmented MPC framework designed for adaptive blind
quadruped locomotion, encompassing high-speed movement, uncertain dynamics
adaptation, and reactive obstacle traversal.

  * •

Our contribution uniquely combines stance foot force control with swing foot
reflection, addressing model uncertainties and bridging foot swing and force
control, overcoming inherent challenges in the nominal MPC framework.

  * •

Our framework provides a universal RL module for MPC, realizing zero-shot
transfer across various robot platforms, and showcasing state-of-the-art
performance on the Unitree A1, Go1, and AlienGo robots.

The paper is organized as follows: Sec. II presents our novel RL augmented MPC
framework. Then, experimental validation is presented in Sec. III. Sec. IV
provides concluding remarks.

Figure 2: System architecture of the proposed framework. The high-level
module, framed in blue, includes the adaptive behavior policy and locomotion
control module, operating at 33Hz. The low-level module, running at 1kHz,
includes leg control (using Jacobian and IK), state estimation, and the
robot’s hardware. The $F/m$ block normalizes the MPC force command into
accelerations as a robot-agnostic input to the adaptive behavior policy.

##  II Proposed Approach

###  II-A System Overview

Illustrated in Fig. 2 is the overall system architecture, which is built upon
[7]. The user provides velocity commands to the robot, while an event-driven
finite state machine determines the gait schedule. Within the locomotion
control module, the MPC is in charge of stance foot control. In contrast, the
swing foot control determines the desired foot positions $p_{f}$. Force
commands $F$ are converted into joint torques using the Jacobian, and
concurrently, the desired foot positions are mapped to corresponding joint
angles through inverse kinematics. A Kalman filter facilitates state
estimation, delivering proprioception data to both the locomotion control and
the adaptive behavior policy.

Central to this system is our innovative adaptive behavior policy. Its primary
aim is to impose supplementary actions onto the baseline MPC framework to
bring adaptive behavior to the robot while ensuring performance across
multiple robot platforms. This policy processes past commands, proprioception,
acceleration from MPC force commands and heuristic foot placement. The result
is dynamic compensation (explained in Sec. II-B) for the MPC and an offset
joint angle $\Delta q$ for swing foot reflection (explained in Sec. II-C). The
adaptive behavior policy (explained in Sec. II-D), which learns to synthesize
both the dynamics compensation essential for force control and the reaction
for swing foot trajectory, given gait schedule and velocity commands. More
than mere compensation, our policy amplifies the agility and robustness of the
locomotion. Importantly, it achieves broad generalizability across different
robot platforms without resorting to domain randomization.

###  II-B MPC with Dynamics Compensation

To address the challenges of model uncertainties while retaining the
generalizability across different robot platforms, we build upon the linear
MPC setup in [1]. Model uncertainties inherently sprout from the
simplifications when abstracting from full-order dynamics, for instance, leg
dynamics, and joint-level torque mapping. Furthermore, unpredictable external
perturbations such as varying loads and terrain conditions introduce unknown
forces and moments acting on the robot’s body, contributing to the model’s
uncertainty. To capture these unknown dynamics, we introduce time-varying,
locally-linear acceleration terms to incorporate into the linearized
continuous-time state space equation. The dynamics compensation terms are
encapsulated $\Delta\bm{\alpha}$ and $\Delta\bm{a}$, which represent angular
and linear accelerations respectively in the continuous-time state space
equation:

$\frac{d}{dt}\left[\begin{array}[]{c}\bm{\Theta}\\\ \bm{p}\\\ \bm{\omega}\\\
\bm{\dot{p}}\\\
\end{array}\right]=\left[\begin{array}[]{cccc}\bm{0}_{3}&\bm{0}_{3}&\bm{R}_{z}(\psi)&\bm{0}_{3}\\\
\bm{0}_{3}&\bm{0}_{3}&\bm{0}_{3}&\bm{I}_{3}\\\
\bm{0}_{3}&\bm{0}_{3}&\bm{0}_{3}&\bm{0}_{3}\\\
\bm{0}_{3}&\bm{0}_{3}&\bm{0}_{3}&\bm{0}_{3}\\\
\end{array}\right]\left[\begin{array}[]{c}\bm{\Theta}\\\ \bm{p}\\\
\bm{\omega}\\\ \bm{\dot{p}}\\\ \end{array}\right]+\\\
\left[\begin{array}[]{ccc}\bm{0}_{3}&...&\bm{0}_{3}\\\
\bm{0}_{3}&...&\bm{0}_{3}\\\
\bm{\hat{I}^{-1}}[\bm{r}_{1}]_{\times}&...&\bm{\hat{I}^{-1}}[\bm{r}_{n}]_{\times}\\\
\bm{1}_{3}/m&...&\bm{1}_{3}/m\\\
\end{array}\right]\left[\begin{array}[]{c}\bm{F}_{0}\\\ .\\\ .\\\
\bm{F}_{n}\end{array}\right]+\left[\begin{array}[]{c}0_{3\times 1}\\\
0_{3\times 1}\\\ \Delta\bm{\alpha}\\\ \Delta\bm{a}+\bm{g}\end{array}\right]$
(1)

where $\bm{\Theta}$ represents the robot’s orientation as a vector of Euler
angels $[\phi,\theta,\psi]^{T}$, $R(\psi)$ is the rotation matrix
corresponding to the yaw angle $\psi$, $\bm{p}$ and $\bm{\dot{p}}$ is the COM
position and velocity of the robot, $\omega$ is the angular velocity of the
robot, $r_{i}$ is the vector from the robot’s COM to foot $i$, $F_{i}$ is the
ground reaction force for leg $i$, $I$ is the inertia, $m$ is the mass of the
robot, and $\bm{g}$ is the gravity term.

Equation (1) can be rewritten with an auxiliary state to represent the
dynamics into a convenient state-space form:

$\displaystyle\dot{\bm{x}}_{c}(t)=\bm{A}_{c}(\psi,\Delta\bm{\alpha},\Delta\bm{a})\bm{x}_{c}(t)+\bm{B}_{c}(\bm{r}_{1...n},\psi)\bm{u}(t)$
(2)

where

$\bm{x}_{c}(t)=\left[\begin{array}[]{cccccc}\bm{\Theta}&\bm{p}&\bm{\omega}&\bm{\dot{p}}&1\end{array}\right]^{T}\in\mathbb{R}^{13}\\\
\bm{A}_{c}(t)=\left[\begin{array}[]{ccccc}\bm{0}_{3}&\bm{0}_{3}&\bm{R}_{z}(\psi)&\bm{0}_{3}&\bm{0}_{3\times
1}\\\ \bm{0}_{3}&\bm{0}_{3}&\bm{0}_{3}&\bm{I}_{3}&\bm{0}_{3\times 1}\\\
\bm{0}_{3}&\bm{0}_{3}&\bm{0}_{3}&\bm{0}_{3}&\Delta\bm{\alpha}\\\
\bm{0}_{3}&\bm{0}_{3}&\bm{0}_{3}&\bm{0}_{3}&\Delta\bm{a}+\bm{g}\\\
\bm{0}_{1\times 3}&\bm{0}_{1\times 3}&\bm{0}_{1\times 3}&\bm{0}_{1\times
3}&0\\\ \end{array}\right]\\\ \in\mathbb{R}^{13\times 13}$ (3)

Then, this formulation is discretized and formulated as a QP problem as in
[1].

Opting for accelerations over forces and moments offers a broader view of
disturbances. This is particularly pertinent considering robots vary in their
ability to withstand external forces and moments. Consequently, we choose a
metric that stands independent of the unique mass characteristics specific to
each robot, recognizing that these attributes play a pivotal role in
adaptability. It’s noteworthy that even if the inertia of robots is usually
minimal, any oversight in compensating moments can critically impair
controller efficacy. Angular acceleration $\bm{\Delta\alpha}$, in this
respect, offers a more intuitive and efficient mechanism to modulate the
robot’s orientation. To delve deeper into this nuance: the process of
translating force/moment to acceleration inherently mandates knowledge of the
robot’s mass and inertia. This perspective inherently considers the robot’s
mass and inertia, allowing our approach to seamlessly apply across various
robotic platforms regardless of mass and inertia differences. Moreover, this
design choice also allows a compact formulation of the optimization problem as
the size of all the matrices remains the same as in [1], which also
facilitates onboard computation.

###  II-C Adaptive Foot Swing Reflection

The swing foot in legged robots involves two main components: 1) foot
placement and 2) swing trajectory. The foot placement determines the contact
location which is crucial for stance control, while the swing trajectory would
help the foot reflection to overcome obstacles. Our goal transcends the
conventional scope of adaptive control that solely addresses model
uncertainties. Instead, we seek to attain adaptive behavior in both the foot
placement and swing trajectory, responding reactively to external
disturbances, such as external loads or varying terrain.

In the baseline MPC framework, foot placement follows a predetermined
heuristic[2] based on velocity commands, feedback, and stance time:

$\bm{p}_{heuristic,i}=\bm{p}_{hip,i}+\frac{T_{stance}}{2}\bm{v}+k(\bm{v}-\bm{v}_{cmd})\\\
+\frac{1}{2}\sqrt{\frac{z_{0}}{||g||}}\bm{v}\times\bm{\omega}_{cmd}$ (4)

where $\bm{p}_{hip,i}$ is the hip location in the world frame for leg $i$,
$T_{stance}$ is the scheduled stance phase time, $\bm{v}$ is the velocity of
the robot’s COM, $\bm{v_{cmd}}$ is the velocity command, $z_{0}$ is the
nominal height of locomotion, $\bm{\omega}_{cmd}$ is the yaw rate command and
in this setup, we used a $k$ of $0.03$. We then employ a pre-defined Bezier
curve to interpolate the foot swing trajectory, outputting the $p_{f}$ as the
desired foot location to the swing foot.

Our methodology places a heightened emphasis on adaptive swing reflection.
Unlike traditional approaches that manually adjust trajectories, our system
utilizes an offset joint angle, $\Delta q$, to modulate the nominal swing
trajectory. This doesn’t just modify foot placement; it also dynamically
adjusts its swing trajectory over discrete obstacles. By prioritizing adaptive
swing control, our system offers a more holistic and responsive solution,
synchronously adjusting both the trajectory shape and final foot placement in
real time, all under the guidance of one unified adaptive behavior policy.
When combined with dynamic compensation, our action space enhances the
robustness and agility of legged locomotion.

###  II-D Learning to Synthesize Stance Control and Swing Control for Adaptive
Behavior

In this section, we present our novel learning framework that integrates the
strengths of the traditional MPC control approach with the dynamic
adaptability offered by Reinforcement Learning (RL). While the baseline MPC
framework excels in forward-looking predictions, RL is capable of reasoning
over past experiences. Our primary goal transcends merely addressing model
uncertainties; we strive to synthesize these decoupled control realms (stance
foot control and swing foot control) using RL. This method unravels the
intricate connections between stance foot and swing foot controls. This
synthesis means that when force optimization is constantly evolving, enriched
by insights from the swing foot’s heuristic and past proprioception data. In
parallel, the swing foot’s trajectory and placement are fine-tuned based on
cues from the force optimizations and past proprioception data. This seamless
integration and reciprocal adaptation ensure that the robot exhibits adaptive
behavior under different conditions, underscoring the power and efficacy of
our proposed approach.

####  II-D1 Action Space

Considering that the MPC problem intrinsically incorporates the robot’s mass
properties, we can design robot-agnostic action space to account for model
uncertainty while ensuring generalizability. Our learning module computes both
dynamics compensation components - namely $[\Delta\bm{\alpha},\Delta\bm{a}]$
\- and swing foot reflection in the form of joint angle offset $\Delta q$ from
the nominal swing trajectory. In other words, our adaptive behavior policy
seeks to derive supplementary actions that can be seamlessly layered onto the
locomotion controls. This design ensures improved robustness and agility in
locomotion performance and substantially facilitates the generalizability of
the framework.

####  II-D2 Observation Space

To ensure generalizability, the observation space must remain independent of
the robot’s mass properties, because the MPC problem inherently considers
them. At every time $t$. the policy obtains an observation and performs
supplementary actions to the nominal MPC control framework. As presented in
Fig. 2, the observation for the policy takes a history window of 5 MPC
horizons, capturing parameters including the joint angle and velocities
$\bm{q}$ and $\bm{\dot{q}}$, linear and angular velocities of the robot
$\bm{v}_{com}$ and $\bm{\omega}_{com}$, planned contact boolean of every foot
from the given gait schedule $\bm{s}_{\phi}$, the actual contact state of each
foot from the contact sensor data $\bm{s}_{actual}$, desired COM state from
user input $\bm{v}_{des}$ and $\bm{\omega}_{des}$, the heuristic foot
placement of every foot $\bm{p}_{heuristic}$ and the acceleration from force
commands of MPC optimization $\bm{F}/m$. Similar to dynamics compensation, the
ground reaction force command is expressed in terms of acceleration. This
consideration is pivotal, given that robots come with diverse mass properties,
prompting the MPC to produce force commands on varying scales. By ensuring the
observation space remains agnostic to a robot’s unique internal attributes, we
substantially facilitate the generalizability of our framework.

####  II-D3 Training

In this paper, we employ PPO[32] for training and use the Unitree A1 robot in
the simulation. The reward function at time $t$ is designed to ensure velocity
tracking of the robot while minimizing the energy cost.

$\displaystyle
R(t)=w_{1}r_{survival}+w_{2}r_{velocity}+w_{3}r_{energy}+w_{4}r_{height}$ (5)

where

$r_{survival}=1\\\
r_{velocity}=||\bm{v}_{des}-\bm{v}_{COM}||+||\bm{\omega}_{des}-\bm{\omega}||\\\
r_{energy}=\sum_{i=1}^{12}||\tau_{i}\dot{q}_{i}||t\\\
r_{height}=0.02-||z_{COM}-z_{des}||$ (6)

and $w_{i}$ are the corresponding weight factors for each reward term. We use
the same reward function by varying weights for the 3 applications we validate
our approach. In this work, we approximate the policy using MLPs with hidden
layers of [256, 32, 256] neurons with tanh as activation function. The
training of the MLP is performed offline with numerical simulation by
Pybullet[33].

###  II-E Learning with MPC in Simulation

In the architectural design of our training methodology, the MPC setup holds a
central position. The integration of MPC becomes a computation bottleneck,
prominently manifested in the form of optimization-induced latency. With the
QPOASES[34] solver, the MPC problem’s computation time averages 1$ms$. In
comparison, a simulation step in Pybullet is computed in less than 0.1$ms$.
The implication here is clear: updating the ground reaction force at each
simulation step would drastically decelerate simulation and consequently slow
done training. Therefore, our strategy is to update the MPC problem at a less
frequent rate, while the lower-level joint commands receive updates much more
frequently. As illustrated in Fig. 2, we update the MPC problem every 30 $ms$
(MPC horizon time) to generate the desired ground reaction force command while
the Jacobian maps this force command to joint torques every 1 $ms$. We also
have our policy to be set up with an update rate of 33 $Hz$. This
configuration aligns with our hardware experiments, negating the need to
constantly run the MPC solver for updating the ground reaction force at every
control step. Thanks to this efficient setup, our approach not only expedites
the training but also ensures the sim2real of the framework as it mirrors the
exact setup used in our hardware experiments.

##  III Experimental Validation

In this section, we detail the experimental validation of our framework,
highlighting the enhanced robustness, agility, and adaptive behavior in blind
locomotion. Our approach is rigorously tested across various computation
platforms and robotic systems, with all policy learning conducted offline
using the Unitree A1 robot model. For the comparative study, we maintained
uniformity by adhering to the baseline MPC specifications, which include gait
patterns, foot swing height, MPC weights, and joint PD gains. Throughout the
experiments, the robot was maneuvered using a joystick by the author. A
seamless transfer of all policies to the actual hardware was made possible
thanks to the robustness of the MPC framework and the generalizability of the
learning framework.

###  III-A High Speed Locomotion

The primary aim here is to validate our RL-augmented MPC methodology during
challenging high-speed maneuvers—both running and turning. Due to their stark
dynamical distinctions, these activities provide a rigorous testbed. We use a
flying trot gait for these high-speed maneuvers. The learning domain comprises
dynamics compensation factors and foot placement offsets, delineated as The
angular acceleration within $\pm[2.0,10.0,2.0]rad/s^{2}$, the linear
acceleration bounded by $\pm[4.0,2.0,3.0]m/s^{2}$ and the joint angle set to
$\pm 0.3rad$.

####  III-A1 High Speed Turning

We validated our approach with high-speed turning in place and we are able to
achieve state-of-the-art performance in terms of the turning rate. As
presented in Fig. 3, we ramp up the yaw rate command to $\pm 7rad/s$.
Impressively, the A1 robot hit a peak turn rate of 8.5 $rad/s$ in the counter-
clockwise direction, maintaining an average yaw rate close to $7rad/s$, while
the baseline MPC couldn’t survive this command and failed immediately. As
demonstrated in the support video, the policy doesn’t merely rely on adjusting
foot placement to enhance motion. A significant contribution comes from the
angular acceleration compensation in the roll direction as shown in Fig. 1.
This intelligent adjustment by the policy facilitates smoother turns, enabling
us to register the swiftest turn ever recorded on the Unitree A1 robot.

####  III-A2 High-Speed Running and Steering

Furthermore, we conducted tests involving high-speed running and steering. In
experiments depicted in Fig. 4, when velocity commands were increased to
$3.5m/s$, the baseline MPC failed within a few running steps due to the model
uncertainty at high speed. Contrarily, our learned policy enabled the robot to
withstand disturbances at high speeds, achieving a peak velocity of around
$3m/s$ based on state estimation data. Further tests involved steering at high
speed. We first accelerated the robot to $2.5m/s$ and then applied a yaw rate
command of $0.5rad/s$. Fig. 5 suggests that our learned policy surpassed the
baseline MPC in velocity tracking. Prior to turning, our RL-augmented MPC
tracked the intended velocity of $2.5m/s$, while the baseline lagged, peaking
at roughly $2m/s$. Integrating translation and rotation dynamics, our RL-
augmented MPC deftly steered at $2.5m/s$ as highlighted in Fig. 1. Meanwhile,
the baseline MPC decelerates to $1m/s$ to prevent failure.

Figure 3: Yaw rate plot of high speed turning policy from IMU data Figure 4:
Linear velocity in the body frame of high-speed running policy Figure 5:
Linear velocity in body frame and yaw rate of high-speed running and steering
policy

###  III-B Walking with Significant Model Uncertainty

In this section, we highlight our framework’s effectiveness in managing model
uncertainty and enhancing the robustness of locomotion. In the training setup,
we randomize the velocity commands for the robot in body frame of
$\bm{v}_{x}\in{[-1,1]m/s}$, $\bm{v}_{y}\in{[-0.5,0.5]m/s}$ and
$\bm{\omega}_{z}\in{[-2.0,2.0]rad/s}$ given a trotting gait of 0.3$s$ gait
cycle. We choose trotting as the test gait as it is more difficult than
standing or quasi-static walking. Adding to this complexity, random external
forces and moments are applied, challenging the robot further. Within this
setting, agents learn dynamics compensation and foot placement offsets, with
angular acceleration constrained to $\pm[4.0,10.0,2.0]rad/s^{2}$, linear
acceleration at $\pm[4.0,2.0,8.0]m/s^{2}$ and the joint angle set to $\pm
0.3rad$.

Notably, even though the policy was primarily trained on flat terrain, it
showcased remarkable adaptability on soft, uneven soil terrain, when the A1
robot carries an additional 5$kg$ load, as depicted in Fig. 1. Moreover, the
generalizability of our approach is evident as the policy learned on the A1
robot realizes zero-shot transfer to different robotic platforms like Go1 and
AlienGo while adapting effortlessly to different gait cycle timings (see
support video). This generalizability is credited to our selected action space
and observation space, which is independent of the robot’s internal
properties.

Our framework also excels in compensating for external moments. The exemplary
adaptive behavior of our framework is evidenced in Fig. 6, where we introduce
loads that impose additional moments on the robot’s body. Despite these
challenges, our adaptive behavior policy effectively mitigates the
uncertainties. This performance can be attributed to the incorporation of
angular acceleration as the key dynamics compensation term. Furthermore, our
framework’s computational efficiency stands out, comfortably running on a
Raspberry Pi 4 board on the Go1 robot. The policy inference and MPC
optimization respectively clock in at approximately 1$ms$ and 3$ms$
respectively on the Raspberry Pi 4 board.

TABLE I: Load carrying capacity of the same policy over different robots with trotting gait Robot | Maximum Load Capacity
---|---
AlienGo | 10$kg$
Go1 | 7$kg$
A1 | 10$kg$

(a) Baseline MPC

(b) RL-Augmented MPC

Figure 6: Comparison between baseline and proposed approach in terms of pitch
compensation. The Go1 robot is carrying a load of 5$kg$ and the A1 robot is
carrying a load of 10$kg$. They all run the same adaptive behavior policy
trained on the A1 robot.

(a) Baseline MPC

(b) RL-Augmented MPC

Figure 7: Comparison between baseline and proposed approach for blind stair
climbing. The stair height is approximately 13cm and the foot swing height is
all set to 8cm.

###  III-C Blind Walking over Discrete Terrain

We also validated our framework in the realm of adaptive foot swing
reflection, specifically on discrete terrain. We trained the policy to
traverse over randomly generated discrete terrains ranging from 8 to 12 $cm$
in height, with the robot’s foot swing height target fixed at $8cm$. The
action space for this scenario are angular acceleration within
$\pm[4.0,10.0,2.0]rad/s^{2}$, the linear acceleration at
$\pm[2.0,2.0,2.0]m/s^{2}$ and the joint angle set to $\pm 0.3rad$.

In this setup, dynamic compensation handles unexpected ground contact, and
adaptive foot behavior prevents foot entrapment. Impressively, this policy
seamlessly transitioned to blind stair climbing of a stair height of $13cm$.
As presented in Fig. 7, while the standard MPC often led to the robot’s foot
being trapped, our refined framework adopted adaptive swing trajectories,
facilitating immediate foot reaction upon contact with the discrete
obstacle(see support video). This emergent behavior is learned by the policy,
and the MPC ensures the robot’s stability and movement. The seamless
integration is a testament to the effectiveness of our RL-augmented MPC
approach for adaptive behavior.

##  IV Conclusion

In this research, we have presented an integration of Reinforcement Learning
(RL) and Model Predictive Control (MPC), improving the agility and robustness
of legged locomotion. Through the integration of MPC’s forward-looking
predictions and RL’s ability to reason on past experiences, our framework has
unveiled marked improvements in a robot’s dexterity to traverse intricate
landscapes and handle unexpected perturbations. Notably, the adoption of
adaptive swing foot reflection showcases how the blend of these two
methodologies can lead to real-world locomotion improvements. The extensive
experimental validations presented further underscore the robustness and
generalizability of our proposed RL-augmented MPC framework, setting a new
benchmark for state-of-the-art robotic control systems. We envision that this
synthesized approach will pave the way for future research such as perceptive
locomotion, fostering even more resilient and adaptive behavior for legged
robots.

## References

  * [1] J. Di Carlo, P. M. Wensing, B. Katz, G. Bledt, and S. Kim, “Dynamic locomotion in the mit cheetah 3 through convex model-predictive control,” in _2018 IEEE/RSJ international conference on intelligent robots and systems (IROS)_. IEEE, 2018, pp. 1–9.
  * [2] D. Kim, J. Di Carlo, B. Katz, G. Bledt, and S. Kim, “Highly dynamic quadruped locomotion via whole-body impulse control and model predictive control,” _arXiv preprint arXiv:1909.06586_ , 2019.
  * [3] Y. Ding, A. Pandala, C. Li, Y.-H. Shin, and H.-W. Park, “Representation-free model predictive control for dynamic motions in quadrupeds,” _IEEE Transactions on Robotics_ , vol. 37, no. 4, pp. 1154–1171, 2021.
  * [4] R. Grandia, F. Jenelten, S. Yang, F. Farshidian, and M. Hutter, “Perceptive locomotion through nonlinear model-predictive control,” _IEEE Transactions on Robotics_ , 2023.
  * [5] W. Chi, X. Jiang, and Y. Zheng, “A linearization of centroidal dynamics for the model-predictive control of quadruped robots,” in _2022 International Conference on Robotics and Automation (ICRA)_. IEEE, 2022, pp. 4656–4663.
  * [6] N. Rathod, A. Bratta, M. Focchi, M. Zanon, O. Villarreal, C. Semini, and A. Bemporad, “Model predictive control with environment adaptation for legged locomotion,” _IEEE Access_ , vol. 9, pp. 145 710–145 727, 2021.
  * [7] G. Bledt, M. J. Powell, B. Katz, J. Di Carlo, P. M. Wensing, and S. Kim, “Mit cheetah 3: Design and control of a robust, dynamic quadruped robot,” in _2018 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_. IEEE, 2018, pp. 2245–2252.
  * [8] R. Grandia, A. J. Taylor, A. D. Ames, and M. Hutter, “Multi-layered safety for legged robots via control barrier functions and model predictive control,” in _2021 IEEE International Conference on Robotics and Automation (ICRA)_. IEEE, 2021, pp. 8352–8358.
  * [9] F. Farshidian, M. Neunert, A. W. Winkler, G. Rey, and J. Buchli, “An efficient optimal planning and control framework for quadrupedal locomotion,” in _2017 IEEE International Conference on Robotics and Automation (ICRA)_. IEEE, 2017, pp. 93–100.
  * [10] C. Cao and N. Hovakimyan, “L 1 adaptive controller for a class of systems with unknown nonlinearities: Part i,” in _2008 American Control Conference_. IEEE, 2008, pp. 4093–4098.
  * [11] Q. Nguyen and K. Sreenath, “L 1 adaptive control for bipedal robots with control lyapunov function based quadratic programs,” in _2015 American Control Conference (ACC)_. IEEE, 2015, pp. 862–867.
  * [12] Q. Nguyen, A. Agrawal, W. Martin, H. Geyer, and K. Sreenath, “Dynamic bipedal locomotion over stochastic discrete terrain,” _The International Journal of Robotics Research_ , vol. 37, no. 13-14, pp. 1537–1553, 2018.
  * [13] K. Sreenath, H.-W. Park, I. Poulakakis, and J. W. Grizzle, “Embedding active force control within the compliant hybrid zero dynamics to achieve stable, fast running on mabel,” _The International Journal of Robotics Research_ , vol. 32, no. 3, pp. 324–345, 2013.
  * [14] M. V. Minniti, R. Grandia, F. Farshidian, and M. Hutter, “Adaptive clf-mpc with application to quadrupedal robots,” _IEEE Robotics and Automation Letters_ , vol. 7, no. 1, pp. 565–572, 2021.
  * [15] M. Sombolestan, Y. Chen, and Q. Nguyen, “Adaptive force-based control for legged robots,” in _2021 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_. IEEE, 2021, pp. 7440–7447.
  * [16] M. Sombolestan and Q. Nguyen, “Adaptive force-based control of dynamic legged locomotion over uneven terrain,” _arXiv preprint arXiv:2307.04030_ , 2023.
  * [17] J. Lee, J. Hwangbo, L. Wellhausen, V. Koltun, and M. Hutter, “Learning quadrupedal locomotion over challenging terrain,” _Science robotics_ , vol. 5, no. 47, p. eabc5986, 2020.
  * [18] G. B. Margolis, G. Yang, K. Paigwar, T. Chen, and P. Agrawal, “Rapid locomotion via reinforcement learning,” _arXiv preprint arXiv:2205.02824_ , 2022.
  * [19] L. Krishna and Q. Nguyen, “Learning multimodal bipedal locomotion and implicit transitions: A versatile policy approach,” _arXiv preprint arXiv:2303.05711_ , 2023.
  * [20] V. Makoviychuk, L. Wawrzyniak, Y. Guo, M. Lu, K. Storey, M. Macklin, D. Hoeller, N. Rudin, A. Allshire, A. Handa, _et al._ , “Isaac gym: High performance gpu-based physics simulation for robot learning,” _arXiv preprint arXiv:2108.10470_ , 2021.
  * [21] Y. Ji, Z. Li, Y. Sun, X. B. Peng, S. Levine, G. Berseth, and K. Sreenath, “Hierarchical reinforcement learning for precise soccer shooting skills using a quadrupedal robot,” in _2022 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_. IEEE, 2022, pp. 1479–1486.
  * [22] A. Kumar, Z. Fu, D. Pathak, and J. Malik, “Rma: Rapid motor adaptation for legged robots,” _arXiv preprint arXiv:2107.04034_ , 2021.
  * [23] X. Cheng, A. Kumar, and D. Pathak, “Legs as manipulator: Pushing quadrupedal agility beyond locomotion,” _arXiv preprint arXiv:2303.11330_ , 2023.
  * [24] G. Feng, H. Zhang, Z. Li, X. B. Peng, B. Basireddy, L. Yue, Z. Song, L. Yang, Y. Liu, K. Sreenath, _et al._ , “Genloco: Generalized locomotion controllers for quadrupedal robots,” in _Conference on Robot Learning_. PMLR, 2023, pp. 1893–1903.
  * [25] D. Kang, J. Cheng, M. Zamora, F. Zargarbashi, and S. Coros, “Rl+ model-based control: Using on-demand optimal control to learn versatile legged locomotion,” _arXiv preprint arXiv:2305.17842_ , 2023.
  * [26] X. B. Peng, E. Coumans, T. Zhang, T.-W. Lee, J. Tan, and S. Levine, “Learning agile robotic locomotion skills by imitating animals,” _arXiv preprint arXiv:2004.00784_ , 2020.
  * [27] J. Sacks and B. Boots, “Learning to optimize in model predictive control,” in _2022 International Conference on Robotics and Automation (ICRA)_. IEEE, 2022, pp. 10 549–10 556.
  * [28] Y. Yang, T. Zhang, E. Coumans, J. Tan, and B. Boots, “Fast and efficient locomotion via learned gait transitions,” in _Conference on Robot Learning_. PMLR, 2022, pp. 773–783.
  * [29] Z. Xie, X. Da, B. Babich, A. Garg, and M. v. de Panne, “Glide: Generalizable quadrupedal locomotion in diverse environments with a centroidal model,” in _International Workshop on the Algorithmic Foundations of Robotics_. Springer, 2022, pp. 523–539.
  * [30] A. Pandala, R. T. Fawcett, U. Rosolia, A. D. Ames, and K. A. Hamed, “Robust predictive control for quadrupedal locomotion: Learning to close the gap between reduced-and full-order models,” _IEEE Robotics and Automation Letters_ , vol. 7, no. 3, pp. 6622–6629, 2022.
  * [31] Y. Sun, W. L. Ubellacker, W.-L. Ma, X. Zhang, C. Wang, N. V. Csomay-Shanklin, M. Tomizuka, K. Sreenath, and A. D. Ames, “Online learning of unknown dynamics for model-based controllers in legged locomotion,” _IEEE Robotics and Automation Letters_ , vol. 6, no. 4, pp. 8442–8449, 2021.
  * [32] J. Schulman, F. Wolski, P. Dhariwal, A. Radford, and O. Klimov, “Proximal policy optimization algorithms,” _arXiv preprint arXiv:1707.06347_ , 2017.
  * [33] E. Coumans and Y. Bai, “Pybullet, a python module for physics simulation for games, robotics and machine learning,” http://pybullet.org, 2016–2021.
  * [34] H. Ferreau, C. Kirches, A. Potschka, H. Bock, and M. Diehl, “qpOASES: A parametric active-set algorithm for quadratic programming,” _Mathematical Programming Computation_ , vol. 6, no. 4, pp. 327–363, 2014.