# The inverse problem for fractal curves solved with the dynamical approach
method

Luiz Bevilacqua111Corresponding author. E-mail<EMAIL_ADDRESS>Coimbra
Institute. Graduate School of Engineering, COPPE - Universidade Federal do Rio
de Janeiro, Brazil Laboratory for Interdisciplinary Numerical Modeling, LIMON
- Universidade Federal de Juiz de Fora, Brazil  Marcelo M. Barros  Laboratory
for Interdisciplinary Numerical Modeling, LIMON - Universidade Federal de Juiz
de Fora, Brazil  Graduate School of Civil Engineering, PEC - Universidade
Federal de Juiz de Fora, Brazil  Department of Structural Engineering, ETU -
Universidade Federal de Juiz de Fora, Brazil  Felipe C. V. Venturelli
Laboratory for Interdisciplinary Numerical Modeling, LIMON - Universidade
Federal de Juiz de Fora, Brazil  Graduate School of Civil Engineering, PEC -
Universidade Federal de Juiz de Fora, Brazil

###### Abstract

The purpose of the present paper is to present the main applications of a new
method for the determination of the fractal structure of plane curves. It is
focused on the inverse problem, that is, given a curve in the plane, find its
fractal dimension. It is shown that the dynamical approach extends the
characterization of a curve as a fractal object introducing the effects of
mass density, elastic properties, and transverse geometry. The dynamical
dimension characterizes material objects and suggests that biological
characterization can be much more complete with the methodology presented
here.

Keywords: fractal curves, inverse problem, dynamical approach.

##  1 Introduction

The geometry of curves in a plane is usually analyzed with the help of
mathematical measures theories. Particularly important in recent years is the
theory introduced by Hasdorff aimed at exploring the fractal dimension of a
certain class of curves, namely fractal curves [1]. Let us call fractal a
sequence of curves, called fractal sequence, composed by curves $F_{i}$ with
length $L_{i}$ such that $\lim_{i\rightarrow\infty}L_{i}\rightarrow\infty$,
that can be packed into a two-dimensional box $B$ and it will always be
possible to find a finite neighborhood $\Omega\in B$ such that
$\lim_{i\rightarrow\infty}F_{i}\cap\Omega=\emptyset$. That is, a fractal curve
no matter its length will never fill the box. The success of fractal curves
can partly be explained by their practical applications in representing the
geometry of objects found in nature [2, 3, 4]. The association of biological
functions with fractal dimensions has stimulated research efforts to improve
knowledge in this area [5]. Despite considerable advances, a specific issue
needs to be improved namely the consideration of complementary variables in
real objects such as mass, thickness, and physical properties. This paper
introduces a method capable of taking into account material and geometric
properties that cannot be achieved with the classical approach to fractal
objects.

The method explored in this paper is derived from the dynamical behavior of
simple mass-spring oscillators whose spring is an elastic wire bent with the
same shape as the fractal curve. A typical harmonic oscillator is shown in
Fig.1. The wire is fixed at the left end and carries a mass $m$ at the right
end. Three different displacements can be independently applied to the mass,
as shown in the figure. Therefore, each spring-mass oscillator exhibits three
independent natural frequencies. Each frequency can be excited by an initial
displacement induced by a moment $M$, a vertical force $V$ or a horizontal
force $H$. The frequency response of simple oscillators assembled as proposed
here leads to a deeper source of information as compared with purely geometric
approaches.

Figure 1: Simple oscillator corresponding to the second term of the Koch
triadic. The elastic spring is bent due to: (a) the action of a moment, (b)
the action of a horiontal force or (c) the action of a vertical force, as
initial conditions.

Consider a curve obtained by assembling a collection of elements all with the
same size, to obtain a single chain with no bifurcations. The general
expression for the period corresponding to the $k^{\text{th}}$ term in the
sequence for any of the three inicial conditions in Fig.1 can be shown to be
[6, 7]:

$\frac{T_{k}^{2}}{T_{0}^{2}}=N_{k}\frac{\lambda_{k}}{L_{0}}\Omega_{k},$ (1)

where $T_{0}^{2}=m_{0}h_{0}^{2}L_{0}/E_{0}I_{0}$ is a normalization constant
obtained with material and geometric characteristics of a reference
oscillator. For the initial excitation induced by a moment $M$, the factor
$\Omega_{k}=1$ and, the equation (1) can be rewritten:

$\log\left(\frac{T_{k}}{T_{0}}\right)=\frac{1}{2}\left(\log
N_{k}+\log\frac{\lambda_{k}}{L_{0}}\right),$

or:

$\log\left(\frac{T_{k}}{T_{0}}\right)=\frac{1}{2}(-F_{M}+1)\log\frac{\lambda_{k}}{L_{0}}\qquad\text{where}\qquad\log
N_{k}=-F_{M}\log(\lambda_{k}/L_{0}).$

Now from the relation above we have
$N_{k}(\lambda_{k})^{F_{M}}=(L_{0})^{F_{M}}$ and clearly $F_{M}$ is the
fractal dimension of the curve. Therefore, the equation (1) can be reduced to
the following expression:

$\log\left(\frac{T_{k}}{T_{0}}\right)=\frac{1-F_{M}}{2}\log\frac{\lambda_{k}}{L_{0}}.$
(2a)

The equation (2a) is a fundamental equation to determine the fractal dimension
of plane curves with the dynamical response of simple oscillators.

The other two independent initial displacements may be used to determine the
complementary periods of the harmonic oscillators, namely the displacement
induced by a horizontal force and the displacement induced by a vertical
force, Fig.1(b-c). The expressions are similar, except that the factor
$\Omega_{k}\neq 1$, then the general expression for the relative period is:

$\log\left(\frac{T_{k}}{T_{0}}\right)=\frac{1-F_{H,V}}{2}\log\frac{\lambda_{k}}{L_{0}}+\log\Omega_{k}.$
(2b)

The extra term $\Omega_{k}$ in (2b) comes from the position and orientation of
the elements $\lambda_{k}$ in the geometry. Indeed, the bending moment
distribution induced by the moment $M$, Fig. 1(a), does not depend on the
position and orientation of the elements $\lambda_{k}$. Therefore, the bending
moment is uniformly distributed along the spring length leading to
$\Omega_{k}=1$. For the other cases, corresponding to the initial displacement
induced by a horizontal or vertical force, the bending moment distribution and
consequently the displacement of the mass $m$ depends on the position and
orientation of the elements $\lambda_{k}$ in the spring geometry. This non-
homogeneous distribution of the bending moment leads to the inclusion of a
disturbance term $\Omega_{k}$ in the determination of the period of
oscillation. It can be shown that for regular fractal curves the term
$\Omega_{k}\rightarrow 1$ for large values of $k$ [6, 7, 8, 9]. Therefore the
angular coefficients of the straight lines corresponding to the functions
$\log(T_{k}/T_{0})\times\log(\lambda_{k}/L_{0})$ for the three initial
conditions $M,H$ and $V$ are related with the Hausdorff dimension as
$D_{M}=D_{H}=D_{V}=(1-F)/2$ since all three angular coefficients coincide for
the cases under consideration, $F_{M}=F_{H}=F_{V}=F$.

Therefore, the fractal dynamical dimension depends on the elastic energy
distribution in the oscillator spring. The period of oscillation excited by a
moment $M$ is independent of the orientation of the elements $\lambda_{k}$.
Now, the elastic energy stored in the system generated by the bending moment
is a function of the position and orientation of each element in the curve.
This property is particulary important for detecting random assemblages of
fractal curves. The random assemblages as shown in Fig.2 can be characterized
with the dynamical fractal method.

Figure 2: Different generation laws, in order to build a: (a) regular
sequence, self-similar; (b) non-regular sequence self-similar and (c) random
sequence, non self-similar.

The dynamical responses to initial excitations induced by a moment and by a
horizontal force for a set of curves generated with the triadic initiator
using the random not self-similar pattern, Fig.2(c), are shown in Fig.3. All
terms in the sequence are composed by the assemblage of a well-defined
geometry in decreasing scales but with different orientations in space.
Fig.2(a) represents the classical sequence of the Koch triadic. The Fig.2(b)
corresponds to a self-similar sequence of terms with a peculiar law of
formation. It is similar to the Koch triadic without the symmetry displayed by
the classical Koch triadic. The Fig.2(c) is a non self-similar sequence, built
with the triadic initiator properly reduced to fit the order of the term in
the sequence but with a random orientation.

The dynamical approach proposed here leads to the following conclusions:

  1. 1.

All curves belong to the Koch triadic family since the response to an initial
displacement generated by the moment leads to a straight line as given by
equation (2a) with angular coefficient $-0.1309$ corresponding to
$F_{M}=D=1.2618$ which is the dimension of the Koch triadic Fig.3.

  2. 2.

The response corresponding to an initial displacement induced by a horizontal
force is a randomly distributed set of points for the case 2(c). A straight
line obtained with the least square method applied to this set of points is a
straight line with the same slope as in the previous case Fig.3. Therefore,
this sample consists of a random assemblage of a sequence of Koch triadic
curves properly downscaled.

Figure 3: Log of the normalized period versus the log of the relative length
of the reference element for five random triadic sequences. (a) Initial
condition imposed by a moment at the free end. (b) Initial condition imposed
by a horizontal force at the free end.

Classical methods to determine the fractal dimension are not able to find the
fine structure of fractal curves. Additionally, the dynamical characterization
of fractal curves provides a class of information not accessible to other
methods, namely the material properties of the objects represented by the
curves. More precisely, the fractal dynamic dimension can characterize a
material object by combining its length with cross section geometry, specific
mass and elastic properties. Particularly important is the consideration of
the wire cross section geometry which plays an important role in the geometry
of natural objects as trees and blood vessels.

To analyze the contribution of mass, material properties and cross section
variations to the fractal dimension consider the equation (1) and let the
mass, the elastic modulus, and the moment of inertia of the cross section vary
according to the following power laws:

$m_{k}/m_{0}=(\lambda_{k}/L_{0})^{\nu},\qquad
E_{k}/E_{0}=(\lambda_{k}/L_{0})^{\gamma},\qquad
I_{k}/I_{0}=(\lambda_{k}/L_{0})^{\mu}.$

Introducing the above relations in equation (1) and after similar operations
as before the following equation is obtained:

$\log\left(\frac{T_{k}}{T_{0}}\right)=\frac{1-\overline{F}_{M,H,V}}{2}\log\frac{\lambda_{k}}{L_{0}}+\log\Omega_{k}\qquad\text{where}\quad\overline{F}_{M,H,V}=F_{M,H,V}+\mu+\gamma-\nu$
(2c)

Clearly the fractal dynamical dimension depends on a larger group of
parameters as compared with the classical geometric approach. This is a
particularly important result. Indeed, with an appropriate variation of the
cross section the dynamical fractal dimension of a certain sequence may vanish
despite of the fact that the classical geometric dimension is different from
zero (Fig.4).

Figure 4: The dynamical fractal dimension may consider additional geometric
properties of the elements in a sequence as for the case shown above
$\overline{F}_{M,H,V}=(1+\mu-D)/2$.

##  2 The inverse problem

The basic ideas about the new approach to analyzing the complex fractal
dimension of material objects were introduced in the previous section.
However, a critical question for practical applications concerns the inverse
problem. That is, given an element find its fractal dimension, if any. Since
the dynamical dimension $\overline{F}_{M,H,V}$ depends on four parameters
there is no unique solution. The classical fractal dimension in the sense of
the Hausdorff metric is obtained if $\mu=0,\gamma=0$ and $\nu=0$. It is not
the purpose of this paper to discuss the uniqueness problem posed by the
addition of three new parameters to the classical fractal dimension. The main
purpose here is to use the dynamic dimension to obtain the fractal dimension
of a given sequence.

With the dynamical approach it is possible to determine the fractal dimension
of a given sample by exploring variables that are not accessible to classical
methods. In [6, 7, 8, 9] it was shown that the dynamical dimension of a
sequence of subsamples obtained by successively cutting the previous one from
the original sample led to its fractal dimension.

In this section we present other ways to determine the fractal dimension of a
given curve.

###  2.1 The hidden sequence

Consider a plane curve $C$ as shown in Fig.5(a). There are two fundamental
hidden sequences. The first consisting of a collection of curves
$C_{k}(k=1,2,\dots,m)$ obtained by taking points $P_{k,1},P_{k,2},\dots
P_{k,N_{k}}$ on $C$ such that the distance
$\overline{P_{k,i}P_{k,i+1}}=R_{k}$. For each value $R_{k}$ coresponds a term
of the hidden sequence as shown in Fig.5(b,c,d). The second consists of a
collection of curves $C_{k}$ defined by points on the curve $C$ such that the
length $S_{k}$ measured on the curve $C$ is the same for any two consecutive
points $P_{k,i},P_{k,i+1}$, as shown in the Fig.5.

The dynamical dimension is obtained from the hidden sequences. The results are
shown in Fig.5(e) for the terms in a hidden sequence with constant segment
length $R_{k}$ and in the Fig.5(f) for the terms in the hidden sequence with
constant length along the curve $S_{k}$. The results are compared. The first
option $R_{k}=constant$ for the given element of the hidden sequence displays
a larger dispersion compared to the second criterion taking constant $S_{k}$.
However, the solution obtained with $R_{k}=constant$ leads to a better
approximation with $F_{R}\approx 1.26$ representing a relative error equal to
0.1%.

Figure 5: The sample (a) is decomposed in a sequence of curves each with a
characteristic element length $R_{k}$ in (b,c,d). Dynamical behaviors
associated with the constant radius approach (e), and with the constant length
along the curve approach (f). Illustration of the constant along length
approach to defining the element length $l_{k,i}$ (g)

Now the expressions (2a) and (2b) do not apply for the case $S_{k}=constant$
because for these groupings the elements $l_{k}$ are not constant. Let us call
the elements of this type of grouping $l_{k,i}$. The expression (1) must be
rewritten:

$\frac{T_{k}^{2}}{T_{0}^{2}}=\sum_{i=1}^{N_{k}}\frac{l_{k,i}}{L_{0}}.$

The $l_{k,i}$ elements are assumed to have all the same material properties
and the same cross section dimensions. Now let
$l_{k,i}=\overline{\alpha}_{i}l_{k}$:

$\frac{T_{k}^{2}}{T_{0}^{2}}=\sum_{i=1}^{N_{k}}\overline{\alpha}_{i}\frac{l_{k}}{L_{0}}=\left(\sum_{i=1}^{N_{k}}\frac{\overline{\alpha}_{i}}{N_{k}}\right)N_{k}\frac{l_{k}}{L_{0}}.$

With $\sum_{i=1}^{N_{k}}\overline{\alpha}_{i}/N_{k}=\alpha_{k}$ the expression
above reads:

$\frac{T_{k}^{2}}{T_{0}^{2}}=\alpha_{i}N_{k}\frac{\lambda_{k}}{L_{0}},$

from which follows

$\log\left(\frac{T_{k}}{T_{0}}\right)=\frac{1}{2}(1-F^{*})\log\frac{l_{k}}{L_{0}}\qquad\text{where}\quad
F^{*}=\frac{\log(\alpha_{k}N_{k})}{l_{k}/L_{0}}.$

It is also possible to write:

$\log\left(\frac{T_{k}}{T_{0}}\right)=\frac{1}{2}[1-(F-\delta)]\log\frac{l_{k}}{L_{0}},$
(3)

where

$\delta=\frac{\log(\alpha_{k})}{\log(l_{k}/L_{0})}\qquad\text{and}\quad
F=\frac{\log(N_{k})}{\log(l_{k}/L_{0})}.$

Since $\alpha_{k}\leq 1$ and $\delta\geq 0$ the angular coefficient of the
line given in (3), which is the dynamical dimension, is smaller than the
corresponding classical fractal dimension $D$. The approach taking the
approximation sequence constructed with points $P_{i}$, $P_{i+1}$ on the
sample leads to a shorter total length since $\overline{P_{i}P_{i+1}}\leq
S_{i}$ for all $i$, where $S_{i}$ is the length measured on the curve.
Therefore, $\alpha_{k}$ is less than one for all cases and $\delta$ is always
positive. This means that $F^{*}$ is always less than the fractal dimension
$F$. For the present case we find $F^{*}=1.20$ leading to the value of
$\delta=0.03$.

The function $\log(T_{k}/T_{0})=f\left(\log(\lambda_{k}/L_{0})\right)$
obtained from the virtual sequence with the first option $R_{k}=constant$,
$k=1,2,\dots,n$, presents a greater dispersion as compared with the second
criterion taking $S_{k}=constant$. However, the solution obtained taking
$R_{k}=constant$ leads to a good approximation to the Hausdorff dimension,
namely $D=F\approx 1.26$ with a relative error equal to 0.1%. The
approximation with $S_{k}=constant$ leads to $F^{*}\approx 1.20$. But $F^{*}$
does not represent the Hausdorff fractal dimension, it is a peculiar fractal
dimension associated with the dynamical behavior of such particular
oscillators sequence.

###  2.2 Variation of the physical properties of the hidden sequence

Instead of determining the fractal dimension of a given sample it is often
convenient to find the upper and lower bounds of the hidden dimension.

Now the periods of the oscillators corresponding to the hidden sequence depend
on the moment of inertia of the wire cross-section. If the inertia of the
cross-section varies according to the law
$I_{k}=I_{0}(\lambda_{k}/L_{0})^{\mu}$ where $\lambda_{k}$ is the length of
the characteristic element, the period can be determined by the equation:

$\log\frac{T_{k}}{T_{0}}=\frac{(1-\mu-D)}{2}\log\left(\frac{\lambda_{k}}{L_{0}}\right)+\log\Omega_{k}$
(4)

This corresponds to an approximated straight line with angular coefficient
$(1-\mu-D)/2$. The fractal dimension $D$ is a fixed value depending on the
geometry of the plane curve. Let us determine the evolution of
$\log(T_{k}/T_{0})\times\log(R_{k}/R_{0})$ for different values of $\mu$ as
shown in Fig.6. For $\mu=-0.4$ the function (4) has a positive slope and for
$\mu=0.1$ it has a negative slope.

Figure 6: Period as function of the characteristic element $\lambda_{k}$ for
two sequences consisting of oscillators with different natural frequencies.

Now if for a given value $\mu=\mu^{*}$ the slope is equal to zero,
$(1-\mu^{*}-D)=0$, the fractal dimension is determined, $D=1-\mu^{*}$. From
the experiments above it is possible to say that $1.4>D>1.1$. Therefore, it is
possible to determine the upper and lower bounds of the fractal dimension. The
above results provide a first approximation of the range where the fractal
dimension should be found. More accurate results can be obtained by adjusting
the values of $\mu$ and observing the change in sign of the angular
coefficient. The specific value of the exponent $\mu^{*}$ defining the
characteristic of the wire cross-section for each sequence defines the
characteristic transition between dynamical responses of the sequence of
oscillators. That is, if the exponent that governs the geometry of the cross-
sections keeps the periods constant for the sequence of oscillators, the
corresponding fractal dimension is obtained from $\mu^{*}$, provided that the
other properties are kept constant.

###  2.3 Variation of the mass of the oscillators

Instead of defining a hidden sequence it is possible to find the fractal
dimension of a plane curve through the dynamical response of the given sample
by modifying the characteristics of the oscillator. The fractal dynamical
dimension depends on the mass attached at the free end of the harmonic
oscillator as shown in equation (2c). Therefore, the inverse problem can be
easily solved by analyzing the responses of a sequence where all terms
preserve the same geometry, but the mass varies according to an appropriate
law.

Consider the curve in Fig.7. It is a curve assembled with the Koch initiator
on a given scale with the characteristic segment given by $\lambda_{0}$. Now
for each reference mass $m_{k}$ ($k=1,2,\dots,n$) let us build a sequence of
masses $m_{k,\nu}$ for each $\nu=\nu_{1},\nu_{2},\dots,\nu_{m}$ according to
the power law:

$\frac{m_{k,\nu}}{m_{0}}=\left(\frac{\lambda_{k}}{L_{0}}\right)^{\nu},\quad\nu\geq
0.$ (5)

For the initial displacement imposed by a moment the relative period is given
by:

$\log\frac{T_{k}}{T_{0}}=\frac{(1+\nu-D)}{2}\log\left(\frac{\lambda_{0}}{L_{0}}\right)$
(6)

The responses of each series of oscillators for five characteristic values of
$m_{k}$ are given as a function of $\nu$ in Fig.7. They all cross the line
$\log(T_{k}/T_{0})=0$ at $\nu=0.1918$. Since for this point $(1+\nu-D)=0$
according to (6) the fractal dimension of the curve is determined, namely
$D=1.1918$.

Figure 7: (a) Determination of the geometric fractal dimension of a curve
assembled with the initiator of the Koch triadic at different scales. (b)
Variation of the relative period with the mass coefficient $\nu$. Each line
corresponds to a given order of the geometry, namely $k=1,2,3,4,5$.

This method is very accurate and uses the sample itself to determine its
fractal dimension. Besides, the solution exposed here suggests that the
fractal dimension can be determined experimentally. It is not difficult to
build a harmonic oscillator with a spring bent according to the geometry
determined by the fractal curve. The dynamical responses of the oscillator for
a set of masses $m_{\nu},\nu=\nu_{1},\nu_{2},\dots,\nu_{m}$, varying according
to a power law as given in (5) can be used to determine the fractal dimension
with equation (6).

##  3 Conclusion

In the previous sections, the application of a new method to determine the
fractal characteristics of fractal curves was shown. The new method uses the
dynamical behavior of simple oscillators to determine the particular
characteristics of fractal curves. The elastic element of the oscillator is a
linear spring simulating the geometry of the fractal curve. As the linear
oscillator on a plane has three degrees of freedom the information obtained
from the dynamical behavior is equally more complete than the information
obtained with the geometric approach. This fact is clear from the results
obtained for the case of random fractals as explained above.

The inverse problem, that is, the identification of possible fractal
characteristics in a given curve can be explored with the dynamical method
with different approaches since three material characteristics are involved in
the process. The geometry of the wire cross-section together with its
respective specific mass and elastic properties play a decisive role in
determining the fundamental frequencies of the oscillators. Therefore, the
fractal dynamic dimension is more suitable for characterizing material objects
than the classical approach.

To solve the problem posed by the influence of the cross-section on the
dimension of real objects it is suggested to analyze the cross-section
variation as an independent variable. This, however, is not satisfactory,
since the coupling of the two independent dimensions is incomplete. The
dynamical approach automatically maintains coupling providing consistent
characterization of material objects. Networks of blood vessels and branch
bifurcations in plants are characteristic examples of these material objects
[10].

A considerable advantage of the dynamical approach consists of determining the
fractal dimension from a single sample in a fractal sequence, provided that
the given sample contains sufficient information. With the sample it is
possible to build a sequence of curves, that we call the hidden sequence, with
an appropriate criterion as shown in the section 2.1.

The solution of the identification problem can also be obtained by modifying
other variables that interfere with the dynamical response, such as mass. This
approach can lead to the solution with the help of the given sample coupled to
an appropriated sequence of masses composing a corresponding sequence of
oscillators. With this method the original sample remains undisturbed, as
shown in section 2.3.

Therefore, it is also possible to obtain the fractal dimension of a given
sample through experimental procedures with the help of a linear oscillator.
The oscillator is assembled with an elastic spring in the shape of the sample.
The sequence of the fundamental frequencies, the first harmonic, corresponding
to an appropriated sequence of masses discloses the fractal dimension.

The method presented here can be easily extended to the case of curves in
three-dimensional space. The dynamical approach to analyzing fractal curves is
just beginning. It is rich enough to open up a consistent theoretical and
practical investigation of the particular geometric formation of curves and
surfaces. The dynamical dimension is richer than the geometric fractal
dimension and should be seen as an independent theory. Real objects can be
characterized more consistently through the dynamical dimension as it
incorporates in the fractal dimension not only the curve geometry but also
complementary material characteristics such as density, elastic properties and
the geometry of the cross section.

##  4 Acknowledgements

We acknowledge the CNPq (National Research Council) and the ANP (Brazilian
National Agency for Petroleum, Natural Gas and Biofuels) for the support given
to the authors through the Senior Research Grant and the contract ANP
5850.0104954.17.9: ”Theoretical and experimental research on the application
of physical methods for the mitigation of mineral incrustations”,
respectively.

## References

  * [1] Falconer, K.J., Fractal Geometry: Mathematical Foundations and Applications, Wiley, 1990.
  * [2] Mandelbrot, B.B., The fractal geometry of nature, W.H. Freeman and Company, 1982.
  * [3] Feder, J., Fractals, Plenum Press, 1988.
  * [4] Turcotte, D.L., Fractals in Geology and Geophysics, Pure and Applied Geophysics, 131, p.171-196, 1989.
  * [5] Bassingthwaighte, J.B., Liebovitch, L.S., West, B.J., Fractal physiology, Oxford University Press, 1994.
  * [6] Bevilacqua, L., Barros, M.M., Dynamical fractal dimension: direct and inverse problems, IUTAM Symp. Dyn. Cont. Nonlin. Syst. Uncert., p.127-136, 2007
  * [7] Bevilacqua, L., Barros, M.M., Galeão, A.C.N.R, Geometry, Dynamics and Fractals, J. Braz. Soc. Mech. Sci. Engin., 30(1), p.11-21, 2008.
  * [8] Bevilacqua, L., Barros, M.M., Dynamical characterization of mixed fractal structures, J. Mech. Mat. Struct., 6(1-4), p.51-69, 2011.
  * [9] Bevilacqua, L. Barros, M.M,, Silva, G.M.A, Dynamical characterization of fractal objects: determination of the fine fractal topology using the energy cover, https://arxiv.org/abs/1408.2448, 2014.
  * [10] Barros, M.M., Bevilacqua, L., Elastic fractal trees: a correspondence among geometry, stress, resilience and material quantity, J. Braz. Soc. Mech. Sci. Engin., 37, p.1479–1483, 2015.


# Bayesian Optimisation for Robust Model Predictive Control
under Model Parameter Uncertainty

Rel Guzman, Rafael Oliveira, and Fabio Ramos The authors are with the School
of Computer Science, the University of Sydney, Australia. Rafael Oliveira is
also with the Australian Research Council’s centre in Data Analytics for
Resources and Environments (DARE) and the University of Sydney’s Brain and
Mind Centre, and Fabio Ramos is also with NVIDIA, USA. {rel.guzmanapaza,
rafael.oliveira<EMAIL_ADDRESS>
###### Abstract

We propose an adaptive optimisation approach for tuning stochastic model
predictive control (MPC) hyper-parameters while jointly estimating probability
distributions of the transition model parameters based on performance rewards.
In particular, we develop a Bayesian optimisation (BO) algorithm with a
heteroscedastic noise model to deal with varying noise across the MPC hyper-
parameter and dynamics model parameter spaces. Typical homoscedastic noise
models are unrealistic for tuning MPC since stochastic controllers are
inherently noisy, and the level of noise is affected by their hyper-parameter
settings. We evaluate the proposed optimisation algorithm in simulated control
and robotics tasks where we jointly infer control and dynamics parameters.
Experimental results demonstrate that our approach leads to higher cumulative
rewards and more stable controllers.

##  I Introduction

Stochastic model predictive control (MPC) methods have been successfully used
in many applications, from steady-state control to path planning in robotics
[1]. These methods rely on a model of the system dynamics to obtain optimal
control strategies, or actions, by weight-averaging sampled trajectories at
each time step. There have been several variations of these methods, including
path integral and cross-entropy approaches [2]. Unfortunately, because
stochastic MPC relies on a dynamics model, its robustness degrades when such a
model is different from the true dynamics. The controller could misrepresent
the environment, wrongly optimising a performance function on an alternative
dynamical system, which might not translate to optimal actions on the real
system. The broad term for this issue is reality gap, which appears in sim-to-
real transfer in the context of robotic simulators [3, 4] and system
identification as typical cases. Despite the reliance on an accurate model,
model-based control and model-based reinforcement learning typically require
fewer interactions with the environment, which makes them more suitable for
real-world robotics applications than model-free approaches [5].

In parallel developments, Bayesian optimisation (BO) as a black-box
optimisation method [6] has been applied to control and robotics to optimise
hyper-parameters in environment control [7, 8] and real robotic manipulators
[9]. The conventional homoscedastic BO approach, however, is based on the
strict assumption of constant observation noise. Heteroscedastic noise models
can better quantify uncertainty in practical control problems by considering
input-dependent noise. In [10], varying noise levels were discovered for a
stochastic MPC hyper-parameter $\sigma_{\epsilon}$ that determines control
variance in a model predictive path integral (MPPI) [2] controller for
benchmark control problems.

As a way to deal with the reality gap, domain randomisation has been applied
to randomise simulators to expose the robot to different scenarios instead of
learning controllers on a single simulated scenario [11]. Domain randomisation
could be applied at any component of the reality gap, and it has been used to
simulate environments [12] and to optimise the robot model parameters [4, 13].

In this paper, we address the reality gap problem in stochastic MPC by
deriving a BO framework to learn controllers which are robust under model
uncertainty. Our framework deals with robustness since it does not have to be
precise regarding the dynamics model parameters. We perform adaptive domain
randomisation to automatically estimate probability distributions over model
parameters solely based on task performance data while jointly optimising the
controller’s hyper-parameters. To account for non-homogeneous response noise,
we extend the framework in [10] which employs heteroscedastic BO to tune MPPI
controllers. Our contributions can be summarised as follows:

  * •

a framework for tuning stochastic MPC while jointly estimating probability
distributions of the transition model parameters, which is based only on
observing rewards;

  * •

an analysis on whether capturing the uncertainty of dynamic parameters leads
to better performance;

  * •

experimental results on benchmark simulated classic control and robotic
problems.

##  II Background

###  II-A Stochastic Model Predictive Control

Model predictive control (MPC) consists of optimising robot actions over a
horizon $T$ based on the idea of solving inner cost optimisation problems over
predicted trajectories. MPC returns a next optimal action $a^{*}$ that is sent
to the system actuators. A stochastic MPC method models disturbances as random
variables. At each time step $t$, stochastic MPC generates sequences of
perturbed actions $V_{t}=\\{a_{i}^{*}+\epsilon_{i}\\}_{i=t}^{t+T}$ where
$\epsilon_{i}\sim\mathcal{N}(0,\sigma_{\epsilon}^{2})$, based on a roll-over
sequence of optimal actions $\\{a_{i}^{*}\\}_{i=t}^{t+T}$ that start as 0.
Each action results in a state produced by a transition or dynamics model
$s_{t+1}=f\left(s_{t},a_{t}\right)$. Action sequences result in a state
trajectory $\mathrm{{S}}_{t}=\\{s_{t+i}\\}_{i=1}^{T}$. Each one has a cost
determined by a cumulative function $C$ with instant cost $c$ and terminal
cost $q$:

$\displaystyle
C(\mathrm{{S}}_{t})=q\left(s_{t+T}\right)+\sum_{i=1}^{T-1}c(s_{t+i})~{}.$ (1)

A stochastic method known as model predictive path integral (MPPI) and its
variations [2] provide optimal actions for the entire horizon working under an
information-theoretic approach. After $M$ rollouts, MPPI updates a sequence of
optimal actions and weights:

$a_{i}^{*}\leftarrow
a_{i}^{*}+\sum_{j=1}^{M}w(\mathrm{{V}}_{t}^{j})\epsilon_{i}^{j}\,,$ (2)
$w(\mathrm{{V}}_{t})=\frac{1}{\eta}\exp\left(-\frac{1}{\lambda}\left(C(\mathrm{{S}}_{t})+\frac{\lambda}{\sigma_{\epsilon}^{2}}\sum_{i=t}^{t+T}a_{i}^{*}\cdot
v_{i}\right)\right)\,,$ (3)

where $j\in\\{1,\ldots,M\\}$ and $\eta$ is a normalisation constant. For the
hyper-parameter temperature $\lambda\in\mathbb{{R}}^{+}$, $\lambda\to 0$ leads
to a single trajectory having higher probability of occurrence, and there is
also the control variance $\sigma_{\epsilon}^{2}$ that results in more varying
and forceful actions when it increases [2]. In this way, both $\lambda$ and
$\sigma_{\epsilon}$ control exploration and exploitation of trajectories. The
other hyper-parameters are horizon $T$ and number of roll-outs $M$.

###  II-B Gaussian Processes

Gaussian processes are non-parametric probabilistic models to approximate
functions [14]. Considering a function $g:{\mathcal{{X}}}\to\mathbb{{R}}$, GPs
are specified by a mean function
${\boldsymbol{\mathbf{m}}}=m({\boldsymbol{\mathbf{{{\boldsymbol{\mathbf{X}}}}}}})$
and a covariance function
$k:{\mathcal{{X}}}\times{\mathcal{{X}}}\to\mathbb{{R}}$. Given a design matrix
of $n$ points ${\boldsymbol{\mathbf{X}}}$, function evaluations
$g({\boldsymbol{\mathbf{X}}})=[g({\boldsymbol{\mathbf{x}}}_{1}),\dots,g({\boldsymbol{\mathbf{x}}}_{n})]^{\mathsf{T}}$
where $\\{{\boldsymbol{\mathbf{x}}}_{i}\\}_{i=1}^{n}\subset{\mathcal{{X}}}$
are Gaussian:

$g({\boldsymbol{\mathbf{{{\boldsymbol{\mathbf{X}}}}}}})\sim\mathcal{N}({\boldsymbol{\mathbf{m}}},{\boldsymbol{\mathbf{{K}}}})~{},$
(4)

where ${\boldsymbol{\mathbf{{K}}}}$ is a $n$-by-$n$ covariance matrix
$[{\boldsymbol{\mathbf{{K}}}}]_{i,j}=k({\boldsymbol{\mathbf{x}}}_{i},{\boldsymbol{\mathbf{x}}}_{j})$.
By denoting a test point as ${\boldsymbol{\mathbf{x}}}_{*}$, we condition
$g_{*}=g({\boldsymbol{\mathbf{x}}}_{*})$ on the observations
${\boldsymbol{\mathbf{y}}}=g({\boldsymbol{\mathbf{X}}})+{\boldsymbol{\mathbf{\nu}}}$,
where
${\boldsymbol{\mathbf{\nu}}}\sim\mathcal{N}({\boldsymbol{\mathbf{0}}},{\boldsymbol{\mathbf{{\Sigma}}}}_{\nu})$
represents observation noise, producing a Gaussian posterior distribution
$g_{*}|{\boldsymbol{\mathbf{X}}},{\boldsymbol{\mathbf{y}}},{\boldsymbol{\mathbf{x}}}_{*}\sim\mathcal{N}(\bar{g}_{*},\operatorname{var}(g_{*}))$.
The posterior infers $g$ at unobserved locations and has the closed form:

$\displaystyle\bar{g}_{*}$
$\displaystyle=m({\boldsymbol{\mathbf{x}}}_{*})+k({\boldsymbol{\mathbf{x}}}_{*},{\boldsymbol{\mathbf{{{\boldsymbol{\mathbf{X}}}}}}})({\boldsymbol{\mathbf{{K}}}}+{\boldsymbol{\mathbf{{\Sigma}}}}_{\nu})^{-1}({\boldsymbol{\mathbf{y}}}-m({\boldsymbol{\mathbf{{{\boldsymbol{\mathbf{X}}}}}}}))$
(5) $\displaystyle\operatorname{var}(g_{*})$
$\displaystyle=k({\boldsymbol{\mathbf{x}}}_{*},{\boldsymbol{\mathbf{x}}}_{*})-k({\boldsymbol{\mathbf{x}}}_{*},{\boldsymbol{\mathbf{{{\boldsymbol{\mathbf{X}}}}}}})({\boldsymbol{\mathbf{{K}}}}+{\boldsymbol{\mathbf{{\Sigma}}}}_{\nu})^{-1}k({\boldsymbol{\mathbf{{{\boldsymbol{\mathbf{X}}}}}}},{\boldsymbol{\mathbf{x}}}_{*})~{},$
(6)

Under a homoscedastic noise assumption,
${\boldsymbol{\mathbf{{\Sigma}}}}_{\nu}=\sigma^{2}_{\nu}{\boldsymbol{\mathbf{{I}}}}$.
In this paper, however, we assume a heteroscedastic, i.e. input-dependent,
noise formulation where
$[{\boldsymbol{\mathbf{{\Sigma}}}}_{\nu}]_{ij}=k_{\nu}({\boldsymbol{\mathbf{x}}}_{i},{\boldsymbol{\mathbf{x}}}_{j})$
and $k_{\nu}:{\mathcal{{X}}}\times{\mathcal{{X}}}\to\mathbb{{R}}$ is a
positive-definite covariance function.

###  II-C Bayesian Optimisation

We use BO [6] as it has been applied in robotics and control to optimise
expensive black-box functions. Given a search space ${\mathcal{{X}}}$, BO
commonly uses a Gaussian process (GP) (Sec. II-B) as a surrogate model
$\mathcal{M}$ to internally approximate an objective function
$g:{\mathcal{{X}}}\to\mathbb{R}$, which is commonly a reward maximisation
problem
${\boldsymbol{\mathbf{x}}}^{*}\in\operatornamewithlimits{argmax}_{{\boldsymbol{\mathbf{x}}}\in{\mathcal{{X}}}}g({\boldsymbol{\mathbf{x}}})$.
Given a set of collected observations $\mathcal{D}_{1:t}$, BO constructs the
surrogate model that provides a posterior distribution over the objective
function $g$. This posterior is used to construct the acquisition function
$h$, which measures both performance and uncertainty of unexplored points. The
next step is optimising the acquisition function, obtaining a sample
$({\boldsymbol{\mathbf{x}}}_{t},y_{t})$. We use the acquisition function known
as upper confidence bound (UCB) [15] due to its simplicity in balancing
exploration and exploitation with the balance factor $\delta>0$. A high value
of $\delta$ leads to more exploration:

$h({\boldsymbol{\mathbf{x}}}):=\mu({\boldsymbol{\mathbf{x}}})+\delta\sigma({\boldsymbol{\mathbf{x}}})~{},\quad{\boldsymbol{\mathbf{x}}}_{t}\in\operatornamewithlimits{argmax}_{{\boldsymbol{\mathbf{x}}}\in{\mathcal{{S}}}}h({\boldsymbol{\mathbf{x}}},\mathcal{M},\mathcal{D}_{1:t})~{},$
(7)

where $\mu({\boldsymbol{\mathbf{x}}})$ is the posterior mean and
$\sigma({\boldsymbol{\mathbf{x}}})$ is the posterior variance of the objective
function $g({\boldsymbol{\mathbf{x}}})$. Both are obtained with a surrogate
model $\mathcal{M}$ trained with $\mathcal{D}_{1:t}$. The next sampling point
${\boldsymbol{\mathbf{x}}}_{t}$ is obtained by solving Eq. 7.

##  III Methodology

This section describes our heteroscedastic Bayesian optimisation method for
adaptive MPC. We start with a description of the transition model and the
adaptation problem we consider and then proceed with adaptive MPC and our
approach to model heteroscedastic noise.

###  III-A Transition Model and Adaptation

We use a dynamics model with model parameters $\theta$ and a deterministic
state transition model $s_{t+1}=f\left(s_{t},a_{t},\theta\right)$ where a
given state $s_{t+1}$ depends only on the previous state $s_{t}$ and the
applied action $a_{t}$. The system state $s\in{\mathcal{{S}}}$ lies in a
continuous state space ${\mathcal{{S}}}$, and $a\in\mathcal{A}$ is an action
specified by a discrete action space. Upon the execution of an action, the
system transitions to the next state $s_{t+1}$, producing an instantaneous
reward $r_{t}:\mathcal{S}\times\mathcal{A}\to\mathbb{{R}}$ that measures the
system performance at a given state and action.

While stochastic MPC optimises actions, we add model parameter randomisation
at each trajectory rollout in order to adapt the controller to the transition
model in real-time. For instance, in the case of the Pendulum problem, the
parameter mass $\theta=m$ could be beta-distributed,
$m\sim\operatorname{Beta}(\alpha,\beta)$, and we would obtain mass
distributions as in Fig. 1b. In such a case, we would have an adaptive
transition model $s_{t+1}=f\left(s_{t},a_{t},m\right)$ and the objective would
be to adapt a prior mass distribution $p(m)$ to a robust mass probability
distribution with expected value around the true mass.

(a) Gamma distributions

(b) Beta distributions

Figure 1: Examples of distributions for different parameters.

Figure 2: Rewards for model parameter distributions.

To handle model uncertainty, we consider model parameters as random variables
with values within a given range. The closer the range is from the true
parameter, the higher the reward as it appears in Fig. 2, which shows
cumulative rewards for possible distributions of the pendulum mass and length
parameterised by $\mu$ and $\sigma$. The red regions present the highest
cumulative rewards. We estimate these optimal parameter distributions via
Bayesian posterior inference.

We define a random vector of model parameters $\boldsymbol{\theta}$ and adapt
them to the heteroscedastic MPC controller. Each model parameter follows a
probability distribution parameterised by $\boldsymbol{\psi}$:

$\boldsymbol{\theta}\sim
p_{\boldsymbol{\theta}}(\boldsymbol{\theta};\boldsymbol{\psi}),\qquad\mathbf{s}_{t+1}=f\left(\mathbf{s}_{t},\mathbf{a}_{t}+\boldsymbol{\epsilon}_{t},\boldsymbol{\theta}\right)~{},$
(8)

where $\mathbf{s}_{t}$ is a vector of states obtained at time $t$, and
$\mathbf{a}_{t}+\boldsymbol{\epsilon}_{t}$ is a vector of perturbed actions.
$\boldsymbol{\theta}$ become transition model inputs. Finally, optimal actions
found by MPC are sent to the system using the adaptive transition model $f$.

###  III-B Adaptive Model Predictive Control

We do not directly aim at finding parameters that match the observed dynamics
as we would do in a system identification approach [16]. Instead, we adapt
model parameter distributions to the controller, which means those resulting
distributions may not be close to their true values. We make sure that we are
using the right MPC hyper-parameters by adapting them to the transition model
parameters. In the case of MPPI [2], the hyper-parameters are the temperature
$\lambda$, noise $\sigma_{\epsilon}$, horizon $T$, and the number of rollouts
$M$ as described in Sec. II-A. They can be collectively described as
$\boldsymbol{\phi}$. In this work, we only work with $\lambda$ and
$\sigma_{\epsilon}$. Fig. 3 shows cumulative rewards for grid search hyper-
parameter combinations of Pendulum and Reacher simulators where optimal
regions are in red.

Figure 3: Rewards for MPPI hyper-parameter combinations.

We optimise the episodic cumulative reward $R$ by jointly optimising
$\boldsymbol{\psi}$, which parameterises the model parameters $\theta$, and
the controller hyper-parameters $\boldsymbol{\phi}$. The optimisation problem
is described as follows:

$\boldsymbol{\psi}^{\star},\boldsymbol{\phi}^{\star}=\operatornamewithlimits{argmax}_{\\{\boldsymbol{\psi},\boldsymbol{\phi}\\}}\,R(\boldsymbol{\psi},\boldsymbol{\phi})~{}.$
(9)

We define $R=\sum_{i=1}^{n_{s}}r_{i}$, where $n_{s}$ is the number of time
steps of an episode. We aim to maximise the expected cumulative reward
$g:=\mathbb{E}[R]$ of the controller under the real transition dynamics, which
is unknown. We approximate the expected cumulative reward $g$ empirically by
averaging over $n_{e}$ episodes.

In order to handle noisy heteroscedastic observations, we define
${\boldsymbol{\mathbf{x}}}=\\{\boldsymbol{\psi},\boldsymbol{\phi}\\}$ as a
variable to optimise the cumulative reward as performance measure $g$ with BO
posterior inference
$g_{*}|{\boldsymbol{\mathbf{X}}},{\boldsymbol{\mathbf{y}}},{\boldsymbol{\mathbf{x}}}_{*}\sim\mathcal{N}(\bar{g}_{*},\operatorname{cov}(g_{*}))$
where ${\boldsymbol{\mathbf{X}}}$ and ${\boldsymbol{\mathbf{y}}}$ is a set of
hyper-parameter observations and $({\boldsymbol{\mathbf{x}}}_{*},g_{*})$ are
the hyper-parameter setting and reward sampled by BO. This way, we are able to
find optimal MPC hyper-parameters adapted to the transition model.

###  III-C Parametric Heteroscedastic Noise Model

Cumulative rewards as a function of hyper-parameters for stochastic MPC
present varying noise [10]. For that reason, a heteroscedastic noise model for
BO optimisation should be considered. We use a heteroscedastic noise model
$\mathbf{\Sigma}_{\nu}$ for the GP defined in Eq. 4. To define
$\mathbf{\Sigma}_{\nu}$, we assume that episodes are executed independently,
which means that
$k_{\nu}({\boldsymbol{\mathbf{x}}},{\boldsymbol{\mathbf{x}}}^{\prime})=0$ for
${\boldsymbol{\mathbf{x}}}\neq{\boldsymbol{\mathbf{x}}}^{\prime}$. Then, we
model only
$k_{\nu}({\boldsymbol{\mathbf{x}}},{\boldsymbol{\mathbf{x}}})=\sigma_{\nu}^{2}({\boldsymbol{\mathbf{x}}})$.
Previously, [10] defined the parametric noise model:

$\displaystyle\sigma_{\nu}({\boldsymbol{\mathbf{x}}})$
$\displaystyle=z\cdot\exp\left({\boldsymbol{\mathbf{\beta}}}^{\mathsf{T}}{\boldsymbol{\mathbf{\rho}}}({\boldsymbol{\mathbf{x}}})\right)+\zeta~{},$
(10)

where ${\boldsymbol{\mathbf{\beta}}}\in\mathbb{{R}}^{m}$, $\zeta\geq 0$ and
${\boldsymbol{\mathbf{\rho}}}:{\mathcal{{X}}}\to\mathbb{{R}}^{m}$ is a feature
map. Small values of $z$ would produce functions that are closer to $\hat{g}$,
and when $z$ is large, the model will account for more outliers. The offset
term $\zeta$ represents minimum homoscedastic noise.

The exponential term ensures that $k_{\nu}$ is always positive definite and
includes a generalised linear model
${\boldsymbol{\mathbf{\beta}}}^{\mathsf{T}}{\boldsymbol{\mathbf{\rho}}}({\boldsymbol{\mathbf{x}}})$
that fits the expected cumulative reward. The feature map
${\boldsymbol{\mathbf{\rho}}}$ could be polynomial or kernel-based. The degree
of the polynomial depends on the cumulative reward function to fit. In our
experiments, we used a polynomial of degree 10, which was sufficient to model
the noise we observed in practice.

Figure 4: Expected reward function for the Half-Cheetah task in (a) using a
10-degree polynomial regression estimate $\hat{g}$, and the fitted noise
$\sigma_{\nu}$ in (b) seen as the green curve.

Having a parametric form for $k_{\nu}$, we learn noise parameters separately
in a two-stage regression problem. With a set of sample inputs
${\boldsymbol{\mathbf{x}}}_{i}\in{\mathcal{{S}}}$, we can approximate the
expected reward $\hat{g}$ as in Fig. 4 using a flexible generalised linear
regression model
${\boldsymbol{\mathbf{\rho}}}:{\mathcal{{X}}}\to\mathbb{{R}}^{m}$ by:

$\begin{split}\hat{g}({\boldsymbol{\mathbf{x}}})\approx{\boldsymbol{\mathbf{\alpha}}}^{\mathsf{T}}{\boldsymbol{\mathbf{\rho}}}({\boldsymbol{\mathbf{x}}})~{}.\end{split}$
(11)

With the estimate $\hat{g}$, we then fit the residuals
$|g({\boldsymbol{\mathbf{x}}})-\hat{g}({\boldsymbol{\mathbf{x}}})|$ with (10)
as a regression problem that leads to finding $\sigma_{\nu}$.

##  IV Related Work

The importance of physics and robotic simulation tools is undeniable because
they make it possible to generate large amounts of data for learning models
[17, 18]. However, there are two main problems with this approach. Simulators
do not fully replicate real hardware behaviour and complex environments, and
there is a limit on the amount of data to collect due to resource limitations.
These problems are related to the reality gap, or sim-to-real problem [3, 12],
and they become noticeable when transferring policies learned in simulators to
real systems.

###  IV-A Model-based Methods and Domain Randomisation

In model-based methods, a dynamics or transition model is fit by maximising
the likelihood of a set of collected trajectories. This optimisation does not
always correlate with optimising the system reward, as demonstrated in [19].

As a way to deal with the reality gap, rather than learning a policy on a
single simulated environment, domain randomisation randomises the simulator to
expose the policy to a wide range of environments at training time [11]. In
the literature, dynamics randomisation has been applied to find approximations
for robust dynamics in common control and robotics problems. For example,
there are generalised transition models that address environment contexts
where a robot’s dynamics could change due to some robotic part malfunctioning
[13]. Other approaches propose inferring simulation parameters based on data
instead of uniform parameter randomisation [3, 4].

Bayesian optimisation has also been applied to adapt domain parameter
distributions during learning, improving sim-to-real transfer for classic
control problems. [12, 20]. In particular, the approach in [12] also jointly
optimises controller parameters. However, they rely on external algorithms
[21] to optimise control policies and do not account for heteroscedasticity in
the reward signal.

###  IV-B Stochastic MPC

Stochastic model predictive control (MPC) has several variants applied from
steady-state control to path planning [2] to robotics [22]. This is due to
their capacity to deal with highly non-linear dynamics and optimise non-convex
reward or cost functions. Moreover, because of the similarities between model
predictive control and model-based reinforcement learning [23], approaches
exploiting both types of methods have been proposed. In [22], for instance,
stochastic MPC together with Gaussian processes allows the system to adapt to
disturbances using the GP to model uncertainty where little data is available
for a robotic arm.

###  IV-C Heteroscedastic Bayesian optimisation

Bayesian optimisation has been applied to model selection and high-dimensional
hyper-parameter tuning, where BO methods form some of the main sample-
efficient, derivative-free optimisation methods used in robotics [24]. There
have been hyper-parameter tuning approaches for MPC using BO, which optimise
back-off terms for thermal energy storage [25], showing considerably optimal
terms in a number of simulations. Other approaches propose MPC as a
performance-driven controller by optimising design parameters with BO [26].
Another method uses heteroscedastic BO for MPC [10], which optimises MPC
hyper-parameters with BO. The difference with this last one is that we add
transition model parameter distributions to the search space.

BO infers samples modelled by a Gaussian process, usually defined with
homoscedastic i.i.d. observation noise [14]. The problem resides when
modelling heteroscedastic noise. There are two main ways to do it:
interpreting the heteroscedasticity using a separate model or integrating the
heteroscedasticity within the GP [27].

One approach for handling heteroscedasticity is to use a separate model for
the noise, which could be less expensive for big data because the system would
need to fit a noise function beforehand [28]. Other approaches are useful for
real-time systems by optimising the noise in real-time. The approach in [27]
adds a diagonal term to the GP covariance matrix to model-independent, non-
identically distributed noise, which makes inference analytically intractable,
requiring approximate solutions. Another approach is to model input-dependent
noise variance with a separate GP and then perform variational inference to
estimate an approximate posterior GP [29]. That approach has been applied to
control learning [30]. However, in addition to the extra computational cost of
variational inference, the method in [30] does not evaluate uncertainty in the
dynamics.

##  V Experiments

We evaluate the performance of the proposed adaptive model predictive control
by solving simulated control and robotic tasks. In the first subsection, we
specify the tasks and then the dynamics model and controller settings used in
the experiments. Next, we evaluate the approach using optimisation methods to
support the use of adaptive MPC. Finally, we present results with robotic
tasks. More details can be found in the video submitted with the paper.

###  V-A Optimisation Variables

We evaluated adaptive MPC in OpenAI111OpenAI Gym: https://gym.openai.com
benchmark tasks: Cartpole, Pendulum, Half-Cheetah, Reacher, and Fetchreach
with dense reward functions. The same functions taken from [10] are applied to
compute instant costs for stochastic MPC.

First of all, different tasks have different constraints for model parameters.
Some parameters have to be positive, like the mass and length, and others can
be defined only within a specific range. We considered using a probability
distribution with positive support for $\boldsymbol{\theta}$. A beta
distribution has support $[0,1]$, which is inconvenient since one has to scale
sampled model parameters if values outside of this range are needed. We also
need to specify ranges for distribution parameters $\boldsymbol{\psi}$, e.g.
if we define model parameters as gamma distributed, we need to specify ranges
for $\alpha$ and $\beta$. We optimise the distribution mean $\mu$ and standard
deviation $\sigma$ and then derive $\alpha$ and $\beta$. We use gamma
distributions for the simplicity of such conversion. We define each model
parameter as gamma-distributed:
$\theta\sim\operatorname{Gamma}(\alpha,\beta)$, where
$\alpha=\frac{\mu^{2}}{\sigma^{2}}$ and $\beta=\frac{\mu}{\sigma^{2}}$.

Next, we define ranges for the parameterisation $\boldsymbol{\psi}$ and the
controller hyper-parameters
$\boldsymbol{\phi}=\\{\lambda,\sigma_{\epsilon}\\}$. These ranges were found
by narrowing down large enough intervals until noticeable optimal regions are
found. Then, we use heteroscedastic BO, to perform episodic cumulative reward
optimisation from Eq. 9 where we jointly optimise $\boldsymbol{\psi}$ and
$\boldsymbol{\phi}$.

###  V-B Experiments on Classic Control Tasks

We aim to find optimal hyper-parameter settings ${\boldsymbol{\mathbf{\phi}}}$
together with model distribution parameters ${\boldsymbol{\mathbf{\psi}}}$.
For Cartpole and Pendulum, $\mu_{l}$ and $\mu_{m}$ are the mean pole length
and mass. For Half-Cheetah and Reacher, we consider the scaling factor
$\kappa$ as a random variable with mean and standard deviation denoted by
subscripts. $\kappa$ multiplies the model parameters, e.g. $\kappa$ with mean
$\mu=1.0$ and $\sigma\to 0$ means that the model parameter corresponds exactly
to its true value. $\kappa_{m}$ and $\kappa_{d}$ are the scaling factors of
the masses and damping ratios for Half-Cheetah and Fetchreach, and
$\kappa_{d1}$ and $\kappa_{d2}$ are scaling factors for two Reacher damping
ratios.

Problem | $n_{e}$ | $T$ | $M$ | Distribution parameter range
---|---|---|---|---
Cartpole | 40 | 10 | 250 |  $\mu_{l}\in{[}0.2,1.5{]}$ $\mu_{m}\in{[}0.1,1.5{]}$ |  $\sigma_{l}\in{[}1e-5,0.1{]}$ $\sigma_{m}\in{[}1e-5,0.1{]}$
Pendulum | 15 | 20 | 400 |  $\mu_{l}\in{[}0.2,2.0{]}$ $\mu_{m}\in{[}0.2,2.0{]}$ |  $\sigma_{l}\in{[}1e-5,0.1{]}$ $\sigma_{m}\in{[}1e-5,0.1{]}$
Half-Cheetah | 25 | 14 | 10 |  $\kappa_{m,\mu}\in{[}0.6,1.2{]}$ $\kappa_{d,\mu}\in{[}0.6,1.4{]}$ |  $\kappa_{m,\sigma}\in{[}1e-5,0.1{]}$ $\kappa_{m,\sigma}\in{[}1e-5,0.1{]}$
Reacher | 20 | 12 | 18 |  $\kappa_{d1,\mu}\in{[}0.1,8.0{]}$ $\kappa_{d2,\mu}\in{[}0.1,8.0{]}$ |  $\kappa_{d1,\sigma}\in{[}0.001,0.6{]}$ $\kappa_{d2,\sigma}\in{[}0.001,0.6{]}$
Fetchreach | 120 | 3 | 12 | $\kappa_{d,\mu}\in{[}1.0,50.0{]}$ | $\kappa_{d,\sigma}\in{[}0.001,0.6{]}$
TABLE I: Ranges of parameters and scaling factors.

We use true parameter values of 1.0 and search spaces within
$\lambda\in{[}1e-5,2.5{]}$, $\sigma_{\epsilon}\in{[}1e-5,4.0{]}$ for the
controller hyper-parameters and distribution parameters with ranges shown in
Table I with $n_{s}=200$ for every task. Results in Fig. 5 show averaged
cumulative rewards for 50 iterations by running both BO versions, and
covariance matrix adaptation evolution strategy (CMA-ES) [31], which was
configured according to [10]. Each BO method had GP kernel hyper-parameters
and noise model globally optimised via maximum likelihood estimation from a
batch of 150 observations. We used an exponential kernel for both. The maximum
likelihood optimisation for finding optimal GP kernel settings is done using
L-BFGS-B with random initialisations. For both BO versions, we used the
exponential kernel
$k(x,x^{\prime})=\sigma_{n}^{2}\exp\left(-\frac{(x-x^{\prime})^{2}}{2\ell^{2}}\right)$.
We set $\Omega:=\\{\sigma_{\nu},\sigma_{n},\ell\\}$ for homoscedastic BO and
$\Omega:=\\{z,\sigma_{n},\ell\\}$ for heteroscedastic BO. For both BO
versions, we used a UCB acquisition function (7) with $\delta=2.0$ optimised
with the same global L-BFGS-B procedure.

We can see that heteroscedastic BO can outperform the others due to the noise
nature. Each optimisation starts at points that give a minimum reward. On
average, adaptive MPC with heteroscedastic BO can find optimal
${\boldsymbol{\mathbf{\psi}}}$ and ${\boldsymbol{\mathbf{\phi}}}$ in fewer
iterations and, in most cases, with less variance.

Figure 5: Averaged cumulative rewards per iteration. The shaded areas denote
two standard deviations. Each method started at a point with minimum expected
cumulative reward.

###  V-C Adapted Parameter Distributions

For Half-Cheetah, the optimised distribution of the scaling factor of the mass
$\kappa_{m}$ yields a reward of $701.29$, which is higher than the reward
obtained with true value. This can be explained because we optimise the reward
for the controller and not for system identification. In Fig. 6, we show the
cumulative reward of using the true parameters and the optimised distribution.
The optimal parameter distribution estimated shows the approximation to the
true value.

Figure 6: Mass parameter scaling factor $\kappa_{m}$ for Half-Cheetah. Figure
7: Optimised scaling factor $\kappa_{d}$ for Fetchreach.

In the case of Fetchreach, the damping ratio plays a major role. We found out
that even wrong dynamics model parameters can lead to better performance. In
Fig. 7, we start at an initial parameter distribution denoted in blue for the
damping ratio $\kappa_{d}$ that gives a maximum reward found via random
search. Then, we use adaptive MPC with heteroscedastic BO to optimise the
parameter distribution. We found that a high damping ratio allows faster
movements, improving reward performance with the true model. That is because,
for an overdamped model, the controller would be prone to apply higher torques
by choosing a larger control variance $\sigma_{\epsilon}$. However, that leads
to uncertain arm movements when the gripper is close to the target.

###  V-D Experiments on a Robotic Simulator

We used a simulated PANDA robotic arm222IssacGym:
https://developer.nvidia.com/isaac-gym with the task of reaching a yellow
target with a gripper in a single-obstacle environment using the MPPI-based
motion planning framework from [32]. We define $n_{e}=10$, $n_{s}=480$, and
the fixed hyper-parameters $T=20$, $M=150$, and $\sigma_{\epsilon}=0.5$.

Trajectory evaluations are done on the GPU, which helped overcome efficiency
issues of MPC. Then, for simplicity, the task only has one obstacle as in Fig.
8a. In this case, we adapt the controller hyper-parameter
$\phi=\lambda\in{[}0.01,2.0{]}$ to the environment by defining distribution-
based length $x$, width $y$ and height $z$ defined by
$\psi=\\{x_{\mu},x_{\sigma},y_{\mu},y_{\sigma},z_{\mu},z_{\sigma}\\}$, and
with true values (0.3, 0.1, 0.6) respectively. Their search space is limited
to 0.02 higher than each true size.

In Fig. 8b, we compare the averaged cumulative reward against the number of
iterations with each method after running 80 iterations. The heteroscedastic
behaviour of $\lambda$ is shown in Fig. 8c, and that is exploited by BOhetero.

(a) Robot task

(b) Reward optimisation

(c) Heteroscedastic behavior

Figure 8: Experiments with the PANDA simulator.

###  V-E Experiments on a Real Robot

We configured a reaching task using the Jaco 2 robot shown in Fig. 10 (left).
We define $n_{e}=8$, $n_{s}=1$, and the fixed hyper-parameters $T=1$,
$M=15000$, and $\sigma_{\epsilon}=1$. The action space is composed of the 6
joint angles. The objective is to reach a fixed red target from Jaco’s
predefined base position. The arm length $d_{1}$ and the front arm length
$d_{2}$ with true values (0.4100, 0.2073) are defined as randomised variables
with $\psi=\\{d_{1,\mu},d_{1,\sigma},d_{2,\mu},d_{2,\sigma}\\}$, and we adapt
them to $\phi=\\{\lambda\\}$ as with the PANDA robot. We define intervals of
$\pm 0.10$ the true lengths, and standard deviations
$d_{1,\sigma},d_{2,\sigma}\in(0.001,0.03)$. The reward is computed as the
negative distance from the gripper to the target.

Figure 9: Jaco robot (left) and reward optimisation (right).

Figure 10: Initial distribution (left) and best inference after 50 iterations
(right).

Fig. 10 shows the best approximation by BOhetero of the length variables to
their true values after 50 iterations as shown in Fig. 10 (right). Optimal
values of those variables should reach point estimates, but optimising point
estimates is out of the scope of this work. Similar optimal distributions are
inferred by both BO versions in this case.

##  VI Conclusions

This work addresses the problem of optimising the performance of stochastic
MPC in uncertain environments. We presented a BO framework that simultaneously
optimises the internal dynamics model and the controller itself. In this
context, we verified the heteroscedasticity in controller hyper-parameters and
dynamics model parameters for control tasks and robotic tasks since
heteroscedastic BO was able to do as well or better than the traditional BO.
Experiments were run with two versions of BO and a non-BO method for reward
optimisation, and results showed that having parameter and object dimension
distributions can lead to improved performance in a few iterations. For future
work, an approach could be to use a non-GP surrogate model for BO, such as a
Bayesian deep neural network. Another point would be to obtain theoretical
results on the effects of heteroscedasticity and on the performance of the
framework.

## References

  * [1] J. Pravitra, K. A. Ackerman, C. Cao, N. Hovakimyan, and E. A. Theodorou, “L1-Adaptive MPPI Architecture for Robust and Agile Control of Multirotors,” 2020. [Online]. Available: http://arxiv.org/abs/2004.00152
  * [2] G. Williams, P. Drews, B. Goldfain, J. M. Rehg, and I. A. Theodorou, “Information-Theoretic Model Predictive Control: Theory and Applications to Autonomous Driving,” _IEEE Transactions on Robotics_ , vol. 34, no. 6, pp. 1603–1622, 2018.
  * [3] X. B. Peng, M. Andrychowicz, W. Zaremba, and P. Abbeel, “Sim-to-Real Transfer of Robotic Control with Dynamics Randomization,” in _2018 IEEE International Conference on Robotics and Automation (ICRA)_ , Brisbane, Australia, 2018, pp. 3803–3810.
  * [4] F. Ramos, R. Carvalhaes Possas, and D. Fox, “BayesSim : adaptive domain randomization via probabilistic inference for robotics simulators,” in _Robotics: Science and Systems (RSS)_ , Freiburg im Breisgau, Germany, 2019\.
  * [5] T. Wang, X. Bao, I. Clavera, J. Hoang, Y. Wen, E. Langlois, S. Zhang, G. Zhang, P. Abbeel, and J. Ba, “Benchmarking Model-Based Reinforcement Learning,” pp. 1–25, 2019. [Online]. Available: http://arxiv.org/abs/1907.02057
  * [6] B. Shahriari, K. Swersky, Z. Wang, R. P. Adams, and N. De Freitas, “Taking the human out of the loop: A review of Bayesian optimization,” _Proceedings of the IEEE_ , vol. 104, no. 1, pp. 148–175, 2016.
  * [7] R. Antonova, A. Rai, T. Li, and D. Kragic, “Bayesian Optimization in Variational Latent Spaces with Dynamic Compression,” 2019. [Online]. Available: http://arxiv.org/abs/1907.04796
  * [8] M. Fiducioso, S. Curi, B. Schumacher, M. Gwerder, and A. Krause, “Safe contextual Bayesian optimization for sustainable room temperature PID control tuning,” _IJCAI International Joint Conference on Artificial Intelligence_ , vol. 2019-August, pp. 5850–5856, 2019.
  * [9] D. Dries, P. Englert, and M. Toussaint, “Constrained Bayesian optimization of combined interaction force/task space controllers for manipulations,” _Proceedings - IEEE International Conference on Robotics and Automation_ , no. May, pp. 902–907, 2017.
  * [10] R. Guzman, R. Oliveira, and F. Ramos, “Heteroscedastic Bayesian Optimisation for Stochastic Model Predictive Control,” _IEEE Robotics and Automation Letters_ , vol. 6, no. 1, pp. 1–1, 2020.
  * [11] J. Tobin, R. Fong, A. Ray, J. Schneider, W. Zaremba, and P. Abbeel, “Domain randomization for transferring deep neural networks from simulation to the real world,” _IEEE International Conference on Intelligent Robots and Systems_ , vol. 2017-September, pp. 23–30, 2017.
  * [12] F. Muratore, C. Eilers, M. Gienger, and J. Peters, “Data-efficient domain randomization with Bayesian optimization,” _IEEE Robotics and Automation Letters_ , vol. 6, no. 2, pp. 911–918, 2021.
  * [13] K. Lee, Y. Seo, S. Lee, H. Lee, and J. Shin, “Context-aware dynamics model for generalization in model-based reinforcement learning,” in _Proceedings of the 37th International Conference on Machine Learning, ICML 2020, 13-18 July 2020, Virtual Event_ , ser. Proceedings of Machine Learning Research, vol. 119. PMLR, 2020, pp. 5757–5766. [Online]. Available: http://proceedings.mlr.press/v119/lee20g.html
  * [14] C. E. Rasmussen and C. K. I. Williams., _Gaussian Processes for Machine Learning_ , 2006.
  * [15] D. D. Cox and S. John, “A statistical method for global optimization,” _Conference Proceedings - IEEE International Conference on Systems, Man and Cybernetics_ , vol. 1992-Janua, pp. 1241–1246, 1992.
  * [16] D. Romeres, G. Prando, G. Pillonetto, and A. Chiuso, “On-line bayesian system identification,” in _2016 European Control Conference (ECC)_ , 2016, pp. 1359–1364.
  * [17] E. Todorov, T. Erez, and Y. Tassa, “MuJoCo: A physics engine for model-based control,” in _IEEE International Conference on Intelligent Robots and Systems_ , 2012, pp. 5026–5033.
  * [18] K. Cranmer, J. Brehmer, and G. Louppe, “The frontier of simulation-based inference,” _Proceedings of the National Academy of Sciences_ , vol. 117, no. 48, 2020.
  * [19] N. Lambert, B. Amos, O. Yadan, and R. Calandra, “Objective mismatch in model-based reinforcement learning,” ser. Proceedings of Machine Learning Research, A. M. Bayen, A. Jadbabaie, G. Pappas, P. A. Parrilo, B. Recht, C. Tomlin, and M. Zeilinger, Eds., vol. 120. The Cloud: PMLR, 10–11 Jun 2020, pp. 761–770. [Online]. Available: http://proceedings.mlr.press/v120/lambert20a.html
  * [20] R. Oliveira, L. Ott, and F. Ramos, “No-regret approximate inference via Bayesian optimisation,” in _37th Conference on Uncertainty in Artificial Intelligence (UAI 2021)_ , 2021.
  * [21] J. Schulman, F. Wolski, P. Dhariwal, A. Radford, and O. Klimov, “Proximal Policy Optimization Algorithms,” _arXiv_ , 2017.
  * [22] A. Carron, E. Arcari, M. Wermelinger, L. Hewing, M. Hutter, and M. N. Zeilinger, “Data-Driven Model Predictive Control for Trajectory Tracking With a Robotic Arm,” _IEEE Robotics and Automation Letters_ , vol. 4, no. 4, pp. 3758–3765, 2019.
  * [23] D. Görges, “Relations between Model Predictive Control and Reinforcement Learning,” _IFAC-PapersOnLine_ , vol. 50, no. 1, pp. 4920–4928, 2017. [Online]. Available: https://doi.org/10.1016/j.ifacol.2017.08.747
  * [24] R. Oliveira, F. H. Rocha, L. Ott, V. Guizilini, F. Ramos, and V. Grassi, “Learning to Race Through Coordinate Descent Bayesian Optimisation,” in _Proceedings - IEEE International Conference on Robotics and Automation_ , 2018, pp. 6431–6438.
  * [25] J. Tobin, L. Biewald, R. Duan, M. Andrychowicz, A. Handa, V. Kumar, B. McGrew, A. Ray, J. Schneider, P. Welinder, W. Zaremba, and P. Abbeel, “Domain Randomization and Generative Models for Robotic Grasping,” in _IEEE International Conference on Intelligent Robots and Systems_ , 2018, pp. 3482–3489.
  * [26] D. Piga, M. Forgione, S. Formentin, and A. Bemporad, “Performance-oriented model learning for data-driven mpc design,” _IEEE Control Systems Letters_ , vol. 3, no. 3, pp. 577–582, 2019.
  * [27] H. Liu, Y.-S. Ong, and J. Cai, “Large-Scale Heteroscedastic Regression via Gaussian Process,” _IEEE Transactions on Neural Networks and Learning Systems_ , pp. 1–14, 2020.
  * [28] F. C. Pereira, C. Antoniou, J. A. Fargas, and M. Ben-Akiva, “A metamodel for estimating error bounds in real-time traffic prediction systems,” _IEEE Transactions on Intelligent Transportation Systems_ , vol. 15, no. 3, pp. 1310–1322, 2014.
  * [29] M. Lázaro-Gredilla and M. K. Titsias, “Variational heteroscedastic Gaussian process regression,” _Proceedings of the 28th International Conference on Machine Learning, ICML 2011_ , pp. 841–848, 2011.
  * [30] S. Kuindersma, R. Grupen, and A. Barto, “Variational Bayesian Optimization for Runtime Risk-Sensitive Control,” in _Robotics: Science and Systems (RSS)_ , N. Roy, P. Newman, and S. Srinivasa, Eds. Sydney, Australia: MIT Press, 2012.
  * [31] D. V. Arnold and N. Hansen, “Active covariance matrix adaptation for the (1+1)-CMA-ES,” in _Proceedings of the 12th annual conference on Genetic and evolutionary computation - GECCO ’10_. Portland, OR: ACM, 2010, p. 385.
  * [32] Bhardwaj, Mohak and Sundaralingam, Balakumar and Mousavian, Arsalan and Ratliff, Nathan and Fox, Dieter and Ramos, Fabio and Boots, Byron, “STORM: An Integrated Framework for Fast Joint-Space Model-Predictive Control for Reactive Manipulation,” in _Conference on Robot Learning (CoRL)_ , 2021\.