Accurate predictions of chaotic motion of a free fall disk

Tianzhuang Xu (徐天壮)1, Jing Li (李靖)1 ***Email address for correspondence:
<EMAIL_ADDRESS>Zhihui Li (李志辉)3,4 and Shijun Liao (廖世俊)1,2 †††Email
address for correspondence<EMAIL_ADDRESS>
1 Center of Advanced Computing, School of Naval Architecture, Ocean and Civil
Engineering, Shanghai Jiaotong University, 200240 Shanghai, China

2 School of Physics and Astronomy

Shanghai Jiaotong University, 200240 Shanghai, China

3 National Laboratory for Computational Fluid Dynamics, 100191 Beijing, China

4 China Aerodynamics Research and Development Center, 621000 Mianyang, China

Abstract It is important to know the accurate trajectory of a free fall object
in fluid (such as a spacecraft), whose motion might be chaotic in many cases.
However, it is impossible to accurately predict its chaotic trajectory in a
long enough duration by traditional numerical algorithms in double precision.
In this paper, we give the accurate predictions of the same problem by a new
strategy, namely the Clean Numerical Simulation (CNS). Without loss of
generality, a free fall disk in water is considered, whose motion is governed
by the Andersen-Pesavento-Wang model. We illustrate that convergent and
reliable trajectories of a chaotic free fall disk in a long enough interval of
time can be obtained by means of the CNS, but different traditional algorithms
in double precision give disparate trajectories. Besides, unlike the
traditional algorithms in double precision, the CNS can predict the accurate
posture of the free fall disk near the vicinity of the bifurcation point of
some physical parameters in a long duration. Therefore, the CNS can provide
reliable prediction of chaotic systems in a long enough interval of time.

Key words free fall; chaos; numerical noises; re-entry

AMS Subject Classifications 37D45, 34C28, 65P20, 65M12

##  1 Introduction

Free fall motion widely exists in nature and industry. For example, falling
leaves and feathers are the common phenomena for people to see in daily life.
A very important application is the spacecraft re-entry, since some objects
are massive in size and high-temperature resistant, and the re-entry may cause
structural, environmental and safety issues on the Earth’s surface. Recently,
the re-entry of China’s Tiangong-1 spacecraft draws much attention to
accurately predicting the uncontrolled trajectories[1]. Six degrees of freedom
(DoF) orbit model is the standard model describing this problem[2, 3].
Nevertheless, accurately modeling a free-falling process is quite challenging
by the 6 DoF model. The chaotic features inherent of this system strongly
depends on initial conditions (SDIC)[4] and, hence, is a great threat to
numerical calculation.

Unfortunately, traditional numerical methods are naturally not “clean”. More
or less noises (i.e., truncation errors and round-off errors) are yielded
during calculation. These noises will increase exponentially in chaotic cases,
due to the SDIC, which was first discovered by Poincaré[5] in 1890 and
developed by Lorenz[6] in 1963, i.e., the so-called “butterfly effect”. Thus,
for a chaotic dynamic system, a tiny variation of the initial condition can
result in significant differences between numerical trajectories after a long
time simulation[7, 8]. Furthermore, Lorenz [9] reported that it is also
sensitive to numerical algorithm. He found that the (maximum) Lyapunov
exponent alters between negative and positive values even when the time step
is very small. Teixeira et al. [10] further investigated the time step
sensitivity of three non-linear atmospheric models utilizing traditional
algorithms in double precision. They made a somewhat pessimistic conclusion
that “for chaotic systems, numerical convergence cannot be guaranteed
forever.”

In fact, the free fall problem itself has been experiencing a long and
complicated history until a widely accepted theory emerges. The qualitative
discussion of it can date back to Maxwell[11]. At that time, little was known
about the nature of the transitions between different modes. Many simple
shapes, such as disks, cylinders, polygons ,cones and even particles have been
studied via experiments or numerical simulations of Navier-Stokes equations,
and complex falling modes were discovered, including tumbling, fluttering,
steady, and chaotic postures[12, 13, 14, 15, 16, 17, 18]. Among these shapes,
the disk is the most well-studied subtopic in this field. A few analytical
models developed by researchers greatly enhanced the understanding of free
fall motion in fluids and promised the possibility to predict the falling
without solving the Navier-Stokes equations which are notorious for the high
computational cost. Kuznetsov[19] organized a comprehensive summary on that
topic.

A direct way to reveal the rules is the experiment. As for the free fall
motion of disks, it mainly focuses on the relationships between physical
parameters and falling modes. Willmarth et al. designed a series of
experiments and measured a phase diagram of fluttering, tumbling and steady
descent according to six related physical parameters[20]. By dimensional
analysis, three similarities were obtained, and, with small ratio between the
thickness and the diameter of the disk, falling modes only depend on the
dimensionless moment of inertia and the Reynolds number. Field et al. further
found a chaotic transition region between fluttering and tumbling[21]. Zhong
et al. recently conducted the most comprehensive experimental research on
free-fall thin disks. They closely studied the relationships between
fluttering (zigzag in their words) free-fall motion and Reynolds number. They
found a critical Reynolds number $Re_{cr}\approx 2000$. The oscillatory
amplitude is proportional to $Re$ below $Re_{cr}$ while invariant beyond
$Re_{cr}$ [22]. Zhong et al also researched the mechanism how two-dimensional
fluttering modes transform to three-dimensional spiral modes[23, 24]. Their
experimental results have been numerically confirmed with immersed boundary-
lattice Boltzmann flux solver[25, 26]. Based on the three-dimensional spiral
modes, Kim et al. further studied the free-fall motion of a pair of rigidly
linked disks. They discovered a mutative falling mode with two disks falling
in helical and conical motions[27]. Recently, Lee et al also studied the
bristled disk where they found that the bristled structure of disks could
strengthen the stability of free-fall motion[28].

Naturally, analytical models were built to explain the various modes and the
bifurcation. Kirchhoff made a pioneering distribution by deriving finite-
dimensional governing equations[29], called Kirchhoff equations, based on an
important fact that velocities of the solid body moving in an ideal
incompressible fluid can be decoupled from the field equations of the fluid
itself. Nevertheless, Kirchhoff equations are only related to ideal fluid
corresponding to conservative systems without considering dissipation. It can
describe the steady fall regimes and to sustain regular or chaotic
oscillations and rotations[30, 31], but still far away from the real motion.
Inherited from Kirchhoff equations, many researchers tried to introduce
appropriate amendments to build a better model, especially on replicating the
fluttering, tumbling, stable and chaotic falling modes and also manifesting
the bifurcation structure between fluttering and tumbling. There are two
models with significant importance, and we introduce them briefly. The first
is the Tanabe–Kaneko model[32]. This model firstly introduces the Joukouvsky
theorem to introduce the effect of circulation and implements a sign function
to relate the lift term with the kinematic information of a disk. Tanabe and
Kaneko explained that the introduction of Joukouvsky theorem and expressed the
circulation with a sign function may give rise to complex dynamics and chaos
during falling in a fluid due to gravity. Though there was some incorrectness
of the Tanabe–Kaneko formulation, including that they omitted the effect of
added mass and Archimedean buoyancy and there was some contradiction between
the coefficients, which was criticized after publication[33, 34], the Tanabe-
Kaneko model qualitatively gives a reasonable picture of possible regimes of
complex dynamics for a disk falling in a fluid. Andersen, Pesavento, and Wang
proposed a more elaborate model to describe the fall of a flat disk or a body
with an elliptic profile in a fluid through a finite-dimensional model[35,
36]. The Andersen-Pesavento-Wang model considers the problems in Tanabe–Kaneko
model and is more coherent with the experimental results and the numerical
results from the direct numerical simulation of the Navier-Stokes equations.

However, Andersen et al. paid more attention on the phenomenology of their
model but lacked close research in the simulation part. Without reliable
numerical simulation near the heteroclinic bifurcation region, they only
mentioned the possibility that the chaotic transition region found in
experiments could be the heteroclinic bifurcation in their model. Moreover, it
turned out in our simulation that the model cannot provide meaningful
prediction in chaos by traditional numerical methods. Thus, we aim to conduct
reliable numerical simulations to closely study the chaotic cases and
heteroclinic bifurcation regions. The motivation of this paper is to implement
a radical numerical strategy to empower the Andersen-Pesavento-Wang model with
the ability to accurately predict trajectories in extremely sensitive cases.

In the present paper, the impact of numerical noises are eliminated by a novel
approach. Liao [37] suggested a numerical strategy in 2009, namely the “Clean
Numerical Simulation” (CNS) [38, 39], to overcome the limitations mentioned
above of traditional algorithms in double precision. Employing the CNS,
reliable/convergent numerical simulations of chaotic dynamical systems can be
obtained in a controllable interval of time $0\leq t\leq T_{c}$, where $T_{c}$
is called the “critical predictable time”. Compared with the traditional
validated numerics methods like interval arithmetics[40], the CNS is a
practical numerical strategy. The implementation of MP makes it easier to use
and computationally cheaper than interval arithmetic while the convergence
checks to determine $T_{c}$ still practically ensure the reliability of
computational results. This method has been proved effective to calculate
reliable trajectories in many chaotic systems, including Lorenz equation[37],
three-body system[42, 43, 44], and also spatio-temporal chaos[45, 46]. For
example, by implementing the CNS, Li et al. successfully found more than 2000
new periodic orbits of the three-body problem which was pointed out by
Poincaré[5] as a classic chaotic system. Most of these periodic solutions are
inaccessible by traditional means[42, 43, 44], which illustrates the
usefulness of the CNS as a powerful tool for reliable investigation of chaotic
systems in physics with high fidelity. As for the spatio-temporal chaos, Lin
et al. [45] used the CNS to control numerical noises smaller than the micro-
level thermal fluctuations, by which it was proved that the inherent micro-
level thermal fluctuations are the root source of macroscopic randomness of
Rayleigh-Bérnard turbulent convection flows. Hu et al. [46] developed a more
efficient algorithm of the CNS to simulate the one-dimensional complex
Ginzburg-Landau equation(CGLE). It further exhibits that CNS method can
accurately maintain both the statistical features of spatio-temporal systems
and the symmetric characteristics of the solutions in which traditional
numerical treatment has failed.

We organize this manuscript as follows. The Andersen-Pesavento-Wang model and
the CNS strategy are briefly introduced in section 2. In the third section, we
demonstrate the sensitivity of the free fall problem to numerical noises and
the advantage of the CNS method by comprehensive comparisons from both chaotic
and periodic simulations. Finally, we close with discussions and concluding
remarks in the last section.

##  2 Mathematical model and numerical algorithm

###  2.1 Andersen-Pesavento-Wang model

Figure 1: The local reference frame $(x^{\prime},y^{\prime})$ fixed with the
disk and the global reference frame $(x,y)$ of a freely falling two-
dimensional desk, where $\theta$ denotes the rotation angle of the disk, $g$
is the acceleration due to gravity, $\Gamma$ is the circulation, respectively.

As shown in Fig. 1, the Andersen-Pesavento-Wang model is a comprehensive
analytical model to predict the trajectories of a freely falling two-
dimensional disk driven by gravity:

$\left\\{\begin{aligned}
I^{*}\dot{V}_{x^{\prime}}&=\left(I^{*}+1\right)\dot{\theta}V_{y^{\prime}}-\Gamma
V_{y^{\prime}}-\sin\theta-F_{x^{\prime}},\\\
\left(I^{*}+1\right)\dot{V}_{y^{\prime}}&=-I^{*}\dot{\theta}V_{x^{\prime}}+\Gamma
V_{x^{\prime}}-\cos\theta-F_{y^{\prime}},\\\
\frac{1}{4}\left(I^{*}+\frac{1}{2}\right)\ddot{\theta}&=-V_{x^{\prime}}V_{y^{\prime}}-\tau,\end{aligned}\right.$
(1)

with the coordinate transformation

$\left\\{\begin{aligned} \dot{x}&=V_{x^{\prime}}\cos\theta-
V_{y^{\prime}}\sin\theta,\\\
\dot{y}&=V_{x^{\prime}}\sin\theta+V_{y^{\prime}}\cos\theta,\\\
\end{aligned}\right.$ (2)

where the dot denotes the derivative with respect to the time $t$,
($x^{\prime},y^{\prime}$) is the local coordinate fixed with the disk, ($x,y$)
is the global (inertia) coordinate,
$\left(V_{x^{\prime}},V_{y^{\prime}}\right)$ is the velocity of disk in the
local coordinate, $\theta$ is the rotation angle of the disk, the circulation
$\Gamma$ is given by

$\Gamma=\frac{2}{\pi}\left(-\frac{C_{T}V_{x^{\prime}}V_{y^{\prime}}}{\sqrt{V_{x^{\prime}}^{2}+V_{y^{\prime}}^{2}}}+C_{R}\dot{\theta}\right),$
(3)

the viscous forces $\left(F_{x^{\prime}},F_{y^{\prime}}\right)$ and torque
$\tau$ are given by

$\left(F_{x^{\prime}},F_{y^{\prime}}\right)=\frac{1}{\pi}\left(A-B\frac{V_{x^{\prime}}^{2}-V_{y^{\prime}}^{2}}{V_{x^{\prime}}^{2}+V_{y^{\prime}}^{2}}\right)\sqrt{V_{x^{\prime}}^{2}+V_{y^{\prime}}^{2}}\left(V_{x^{\prime}},V_{y^{\prime}}\right),\hskip
14.22636pt\tau=\left(\mu_{1}+\mu_{2}|\dot{\theta}|\right)\dot{\theta},$ (4)

respectively. All variables are dimensionless. This model has seven
dimensionless parameters $I^{*},C_{T},C_{R},A,B,\mu_{1},\mu_{2}$, where
$I^{*}=({\rho_{s}b})/({\rho_{f}a})$ with $\rho_{s},\rho_{f}$ being the
densities of disk and fluid and $a,b$ the length of the semi-major and semi-
minor axis of the elliptical disk, other parameters are related with the
geometry of the falling disk. For an elliptical disk, [35, 36] gave

$C_{T}=1.2,C_{R}=\pi,A=1.4,B=1.0,\mu_{1}=0.2,\mu_{2}=0.2$ (5)

by means of fitting the viscous forces and torques of experimental results and
the direct numerical simulation of Navier-Stokes equation. For details, please
refer to [35, 36]. In this paper, we use the Andersen-Pesavento-Wang model
with the above-mentioned values of the parameters
$C_{T},C_{R},A,B,\mu_{1},\mu_{2}$. In this paper we only adjust the parameter
$I^{*}$ to produce different falling modes of the disk.

###  2.2 Clean Numerical Simulation (CNS)

Generally, the CNS is able to obtain long-term reliable results thanks to its
strategy to control the “numerical noises”, say, decrease the truncation
errors to a required level by implementing numerical schemes with extremely
high precision and control the round-off errors within a required range with
all physical/numerical variables and parameters in multiple-precision.

Truncation errors come from the discretization of continuous systems. The
numerical methods have the following general form:

$f(t+h)=f(t)+h\times RHS(t)$ (6)

where $RHS(t)$ denotes the right hand side. It varies according to the
numerical methods. For N-th order Runge-Kutta family method, N-step multi-step
method and N-th order Taylor series method, the right hand side have the
general forms:

$\left\\{\begin{aligned} RHS(t)&=\sum_{i=1}^{N}k_{i}f(t_{k_{i}})+O(h^{N})\\\
RHS(t)&=\sum_{i=1}^{N}k_{i}f(t_{i-N-1})+O(h^{N})\\\
RHS(t)&=\sum_{i=1}^{N}\frac{f^{(i)}(t)}{i!}h^{i-1}+O(h^{N})\end{aligned}\right.$
(7)

$O(h^{N})$ is the order of global truncation errors. The CNS is aimed to
reduce the truncation errors so small that it would not damage the long term
prediction, either by reducing time steps $h$ or increasing $N$ with high
order methods. The round-off errors invariably arise with data are stored in
computers in finite digits. We implemented the multiple-precision libraries
(MP, called the MPFR library in the C language)[41] to also reduce the round-
off errors to a small enough level.

By that, the numerical noises of the simulation are controlled arbitrarily
small. To determine the critical predictable time $T_{c}$, one would conduct
an additional simulation with even smaller numerical noises. In a temporal
dynamic system, we assume that numerical noises grow exponentially within an
interval of time $t\in[0,T_{c}]$:

$\mathcal{E}(t)=\mathcal{E}_{0}\exp(\kappa t),\quad t\in\left[0,T_{c}\right]$
(8)

where the constant $\kappa>0$ is called noise-growing exponent, which is
coherent with the largest Lyapunov exponent (LLE), $\mathcal{E}_{0}$. It
denotes the level of initial noises (i.e., truncation and round-off errors),
and $\mathcal{E}(t)$ is the level of evolving noises of numerical simulation.
Theoretically, a critical level of noise $\mathcal{E}_{c}$ determines the
critical predictable time $T_{c}$ by the equation:

$\mathcal{E}_{c}=\mathcal{E}_{0}\exp\left(\kappa T_{c}\right)$ (9)

It is obvious to tell from the above equation that the smaller initial noise
$\mathcal{E}_{0}$ promises a longer $T_{c}$. Since the true orbits are
impossible to get, the CNS implements a practical method to determine $T_{c}$.
Let $\Phi(t)$ be a numerical simulation reliable in $t\in[0,T_{c}]$ with the
initial noise $\mathcal{E}_{0}$, and $\Phi^{\prime}(t)$ be another simulation
(with the same physical parameters and the same initial conditions) in
$t\in[0,T_{c}]$ with the initial noise $\mathcal{E}^{\prime}_{0}$ that is
several orders of magnitude smaller than $\mathcal{E}_{0}$. According to the
hypothesis that numerical noises grow exponentially, there sure is
$T^{\prime}_{c}>T_{c}$ and $\Phi^{\prime}(t)$ in $t\in[0,T_{c}]$ should be
much closer to the true orbit than $\Phi(t)$. Therefore, we use the
$\Phi^{\prime}(t)$, a better simulation with less numerical noises, to decide
the $T_{c}$ of $\Phi(t)$. After obtaining the $T_{c}$ of $\Phi(t)$, we can
name safely that $\Phi(t)$ is “clean” numerical simulation (CNS) in
$t\in[0,T_{c}]$. Those, as mentioned earlier, provide us a heuristic
explanation of the strategy of the CNS. The CNS also applies to non-hyperbolic
chaotic systems.

We implement the strategy above to obtain the CNS results of the Andersen-
Pesavento-Wang model. Given the discontinuous term in equation (4), we use a
fixed step fourth order Runge-Kutta method with a strict time-step in multiple
precision. To determine $T_{c}$, an additional simulation with even smaller
time steps and more digits to store data in computer to guarantee that the
extra simulation contains even less numerical noises. For the formal
definition of $T_{c}$, we follow the form in the paper [37]. With formula that
$\left|1-\frac{u_{1}}{u_{2}}\right|>\delta,\quad$ at $t=T_{c}$ where
$\delta=1\%$ in this paper, we determine the exact $T_{c}$ of the Andersen-
Pesavento-Wang model. In the following manuscript, we would like to
demonstrate where and how numerical noises can harm the fidelity of trajectory
prediction in the Andersen-Pesavento-Wang model.

##  3 Accurate prediction given by the CNS

###  3.1 Chaotic trajectories

It was regraded that chaotic free-fall motion of disks in the water is long-
term unpredictable by numerical simulation. In this section, we address that
the CNS is able to provide long-term reliable prediction. We study the chaotic
case with $I^{*}=2.2$ reported by Andersen et al. [35, 36]. Four different
types of initial conditions are considered for generality. These four initial
conditions are listed in the form
$(x,y,\theta,V_{x^{\prime}},V_{y^{\prime}},\dot{\theta})$ in Table. 1 with
their corresponding schematic diagram shown in Fig. 2.

No. | angle of attack | rotation | $(x,y,\theta,V_{x^{\prime}},V_{y^{\prime}},\dot{\theta})$
---|---|---|---
1 |  |  | $(0,0,0,0,0.01,0)$
2 | $\surd$ |  | $(0,0,1,0,0.01,0)$
3 |  | $\surd$ | $(0,0,0,0,0.01,1)$
4 | $\surd$ | $\surd$ | $(0,0,1,0,0.01,1)$
Table 1: The initial conditions of the four cases of chaotic falling Figure 2:
The schematic diagram of the four types of initial conditions considered. The
four cases are correspondent to disks that fall from a static state with no
angle of attack, that fall from a static state with an angle of attack, that
fall from a rotational state with no angle of attack and that fall from a
rotational state with an angle of attack Figure 3: (a) The trajectories of the
same initial condition computed by the fourth order Runge-Kutta method in
double precision with different time steps. All trajectories are divergent
from each other, making these trajectories of little value as prediction. (b)
The trajectories of the same initial condition but computed by the CNS. The
black dotted trajectory is computed with smaller time steps and data stored in
more digits, which proves the red trajectory is “clean”.

Traditional methods are powerless as a prediction tool with chaos. In Case 1,
we computed trajectories by the fourth order Runge-Kutta method in double
precision with different time steps from $1\mathrm{e}{-2}$ to
$1\mathrm{e}{-7}$. The trajectories are plotted in Fig. 3, It is found that
all trajectories are divergent from each other after about $250UT$. Though
traditional numerical methods can obtain qualitatively correct chaotic
trajectories, these trajectories are of little use from the aspect of
prediction: it’s impossible to tell which one should be used as the predicted
trajectory. The CNS, on the other hand, is able to give true trajectories: the
trajectory with $h=1\mathrm{e}{-6}$ in 60-digits and that with
$h=1\mathrm{e}{-7}$ in 120-digits. Therefore, for a chaotic system, the
reliable solution is too sensitive to be found by traditional approach. The
good thing is that, CNS can perfectly circumvent this issue.

Figure 4:  (a) The linear relationships between $T_{c}$ and decimal digits
$K$; (b) The linear relationships between $T_{c}$ and the logarithm of steps.

Without loss of generality, we consider the relationships between $T_{c}$ and
numerical noises in case1. We have demonstrated that with $h=1\mathrm{e}{-6}$
in 60-digits the trajectory is “clean”. According to that, two sets of
numerical experiments are conducted to consider the effect of truncation and
round-off errors respectively. The first set is designed to contain only
truncation errors: thus all data are stored in 60-digits while the time steps
are larger than $h=1\mathrm{e}{-6}$. The second set is designed to contain
only round-off errors: thus the same time step $h=1\mathrm{e}{-6}$ is
implemented but data are stored in less digits, from 16-digits(double
precision) to around 25-digits. With the definition of $T_{c}$ we could
compute the $T_{c}$ of these cases by comparing the trajectories with the
“clean” trajectory computed with $h=1\mathrm{e}{-6}$ in 60-digits. With these
two sets of numerical experiments, we confirm that $T_{c}$ pf the Andersen-
Pesavento-Wang model also applies to the exponential growth law. There exist
linear relationships between $T_{c}$ and the logarithm of time steps $h$,
$T_{c}$ and the digits of data $K$, respectively. The quantitative results are
given by the following equations:

$T_{c}\approx\min\left\\{10.95K+96.11,-42.81\log_{10}{h}+112.56\right\\}$ (10)

which provide a rough estimation of the $T_{c}$ of case1. The exact $T_{c}$ is
decided by the minimum of the $T_{c}$ with truncation and round-off errors.
Generally, with smaller time steps $h$ and larger digits $K$, longer reliable
prediction could be obtained. Also it is able to estimate the time steps and
digits needed for a certain $T_{c}$ for the purpose of prediction.

Figure 5: Comparison between trajectories computed by traditional methods
versus that computed by the CNS. (a) (b) (c) (d) correspond to case 1 to 4
respectively.

In all four cases, the CNS is able to provide convergent prediction of the
exact falling trajectories of disks. Let us compare the differences between
the falling trajectories computed by the traditional methods and the CNS. The
details of how trajectories computed by traditional methods diverges from the
true trajectories of all four cases are demonstrated in Fig. 5. The
trajectories in red are computed by traditional methods with
$h=1\mathrm{e}{-6}$ while the trajectories in black are computed by the CNS
with the same time steps. In all our cases, the numerical noises first result
in an observable phase difference of posture angle $\theta$. Then the phase
difference continuously increases and finally develops into a different
falling path, marking the failure of prediction. Among the four cases, the
divergence happens at different times. The first case that disks fall from a
static state with no angle of attack has the latest divergence, because it
takes the longest time for the initial condition to develop into chaos. Except
for that, these four cases have the same qualitative phenomenon.

###  3.2 Heteroclinic bifurcation

In this section, we study the heteroclinic bifurcation region reported by
Andersen et al. in their paper on the analysis of transition between
fluttering and tumbling state[36]. They described it as a sharp transition
while the bifurcation region are sensitive to noises. In their work, Andersen
et al. succeeded in studying the heteroclinic bifurcation to the precision of
$1\mathrm{e}{-5}$. With the help of CNS, we approach the heteroclinic
bifurcation region in even higher precision of $1\mathrm{e}{-11}$ and provide
reliable prediction which was never reported before. It is worth noting that
only the CNS can predict the falling modes in that high precision as shown in
the following example.

Given the dynamical characteristics of the Andersen-Pesavento-Wang model, it
has four steady solutions, corresponding to two fixed points where disk falls
vertically, and gravity is balanced by drag:

$\left[\begin{array}[]{c}v_{x^{\prime}}\\\ v_{y^{\prime}}\\\ \theta\\\
\dot{\theta}\end{array}\right]=\left[\begin{array}[]{c}\mp\sqrt{\frac{\pi}{A-B}}\\\
0\\\ \frac{\pi}{2},\frac{3\pi}{2}\\\
0\end{array}\right]=\left[\begin{array}[]{c}\mp V\\\ 0\\\
\frac{\pi}{2},\frac{3\pi}{2}\\\ 0\end{array}\right]$ (11)

and two other fixed points where the face of disk is normal to the direction
of motion:

$\left[\begin{array}[]{c}v_{x^{\prime}}\\\ v_{y^{\prime}}\\\ \theta\\\
\dot{\theta}\end{array}\right]=\left[\begin{array}[]{c}0\\\
\mp\sqrt{\frac{\pi}{A+B}}\\\ 0,\pi\\\
0\end{array}\right]=\left[\begin{array}[]{c}0\\\ \mp W\\\ 0,\pi\\\
0\end{array}\right]$ (12)

Figure 6: Two characteristic trajectories near the $I^{*}_{c}$. The trajectory
with $I^{*}=1.2190$ fall in fluttering modes while the trajectory with
$I^{*}=1.2191$ fall in tumbling modes.

Andersen et al. discovered that the fluttering modes transform to tumbling
modes via a heteroclinic bifurcation at a specific $I^{*}_{c}$. That is with
$I$ close to $I^{*}_{c}$, all disks with $I^{*}<I^{*}_{c}$ fall in fluttering
modes while disks with $I^{*}>I^{*}_{c}$ fall in tumbling modes, as shown in
Fig. 6. So mathematically, $I^{*}_{c}$ is the critical parameter that decides
disks’ falling modes. Noticing that the periodicity of disks grows
exponentially longer as $I^{*}$ approaches $I^{*}_{c}$. Andersen et al.
explained that the phenomenon results from heteroclinic bifurcation like the
famous Silnikov’s phenomenon.

With the help of the CNS, it is able to obtain an $I^{*}$ which is extremely
close to $I^{*}_{c}$ by a dichotomization process based on the fact that
$I^{*}_{c}$ must lie between fluttering and tumbling modes. For example,
$I^{*}=1.2191466312021015$ is very close to $I^{*}_{c}$ with precision
$1\mathrm{e}{-11}$. Its true trajectory is in fluttering modes computed by the
CNS, plotted in Fig. 7. The CNS results can predict the falling modes even in
extremely high-precision near the heteroclinic bifurcation point.

However, the trajectories computed by the traditional methods with different
time steps ($1\mathrm{e}{-2}$, $1\mathrm{e}{-3}$ and $1\mathrm{e}{-4}$) behave
chaotic as a Computational Chaos (CC) phenomenon, shown in Fig. 7. From the
aspect of prediction, that means traditional methods cannot predict the
falling modes of trajectories near bifurcation point. The CC addresses a loose
similarity between heteroclinic bifurcation region and chaotic transition
region[21]. However, noting that when $I^{*}$ is close to $I^{*}_{c}$, the
periods are elongated. It is easy to find that the “chaotic” trajectories in
Fig. 7 have the characteristics of periodicity, compared with the chaos
discussed in the previous section in Fig. 3. The obvious differences confirm
that the heteroclinic bifurcation region and the chaotic transition region are
different.

Figure 7: The computational chaos resulted from numerical noises. The black
line is the trajectory computed by the CNS, which falls in the fluttering
modes. On the other hand, the trajectory computed by the fourth order Runge-
Kutta method in double precision is neither fluttering nor tumbling. Instead
they are in a special type of chaos resulted.

Considering the special characteristics of the CC here, it is worth for us to
distinguish it from the classic CC proposed by Lorenz[47] as a new type. They
have two main differences. First, this type of CC behaves characteristics of
periodicity as mentioned. The uncertainty only happens near the heteroclinic
bifurcation point that decides the disk would either flutter or tumble in the
next period. Therefore, we refer such chaos as heteroclinic bifurcation chaos.
Second, the CC proposed by Lorenz are only caused by too large time steps and
can be eliminated when traditional methods with smaller time steps are used.
On the other hand, the CC found in the Andersen-Pesavento-Wang model is
unavoidable when computed by traditional methods. No matter what time steps
are picked, the CC always happens. That can be explained by that such CC can
be caused both by truncation errors and round-off errors, so traditional
methods is powerless to control these errors small enough to ensure the
simulation “clean”.

Given the characteristics of heteroclinic bifurcation chaos, experimental
noises, as pointed out by Andersen et al.[36], could also trigger that
heteroclinic bifurcation chaos in that region. We can only conclude the
heteroclinic bifurcation of Andersen-Pesavento-Wang model is actually
different from the chaotic transition region. However, we cannot preclude the
existence of heteroclinic bifurcation chaos discovered here though it is
strange that it was never reported before in any experiments of disks falling
in the water. Whether this type of falling modes exists or the Andersen-
Pesavento-Wang model is not correct as describing the transition between
fluttering and tumbling is still an open question worth further study.

##  4 Discussions and concluding remarks

In this paper, we show that the uncertainty of the numerical simulation of
free fall motion in the water can be overcome by a novel but powerful tool,
the CNS. In the chaotic cases, numerical noises grow exponentially. Hence,
even though they are as small as $1\mathrm{E}{-16}$, they are amplified
exponentially. It is impossible to obtain long-term reliable prediction of
chaotic free fall motion in the water by means of traditional numerical
algorithms in double precision unless the CNS is applied. Besides, in the
heteroclinic bifurcation region, numerical noises can accumulate during the
long period and destroy the bifurcation structure, behaving like a new type of
computational chaos. Under this circumstance, traditional methods cannot even
predict the disk’s falling modes. Although numerical noises can ruin the
validity of this problem, a good way to control them at a small enough level
can still give us a trustful and reliable result in a long enough interval of
time.

In the past, the role of numerical noises was usually neglected. It is
magnificent that we implement physical models to simulate and predict
versatile phenomena in the nature, but it is not recommended to draw
conclusions without awareness of the fidelity of the simulation. Take the re-
entry problem as an example. In this kind of big engineering project, no doubt
a better prediction of the landing spot will bring tremendous economic and
safety benefits. Hence, this work shows people a correct direction to
radically solve chaotic problems by numerical methods, and it is of practical
importance.

## Acknowledgments

This work is supported by the National Natural Science Foundation of China
(No. 91752104) and the project “Development of large-scale spacecraft flight
and reentry surveillance and prediction system” of manned space engineering
technology (2018-14).

## Data Availability Statement

The data that support the findings of this study are available from the
corresponding author upon reasonable request.

## References

  * [1] J. Tang, L. Liu, H. Cheng, S. Hu, J. Duan, Long-term orbit prediction for tiangong-1 spacecraft using the mean atmosphere model, Advances in Space Research 55 (5) (2015) 1432–1444.
  * [2] R. W. Powell, Six-degree-of-freedom guidance and control-entry analysis of the hl-20, Journal of Spacecraft and Rockets 30 (5) (1993) 537–542.
  * [3] C. Zimmerman, G. Dukeman, J. Hanson, Automated method to compute orbital reentry trajectories with heating constraints, Journal of Guidance, Control, and Dynamics 26 (4) (2003) 523–529.
  * [4] V. S. Aslanov, A. S. Ledkov, Chaotic motion of a reentry capsule during descent into the atmosphere, Journal of Guidance, Control, and Dynamics 39 (8) (2016) 1834–1843.
  * [5] H. Poincaré, Sur le problème des trois corps et les équations de la dynamique, Acta mathematica 13 (1) (1890) A3–A270.
  * [6] E. N. Lorenz, Deterministic nonperiodic flow, Journal of the Atmospheric Sciences 20 (2) (1963) 130–141.
  * [7] J. C. Sprott, Elegant chaos: algebraically simple chaotic flows, World Scientific, 2010.
  * [8] J. C. Sprott, J. C. Sprott, Chaos and time-series analysis, Vol. 69, Citeseer, 2003\.
  * [9] E. N. Lorenz, Computational periodicity as observed in a simple system, Tellus A: Dynamic Meteorology and Oceanography 58 (5) (2006) 549–557.
  * [10] J. Teixeira, C. A. Reynolds, K. Judd, Time step sensitivity of nonlinear atmospheric models: numerical convergence, truncation error growth, and ensemble design, Journal of the Atmospheric Sciences 64 (1) (2007) 175–189.
  * [11] J. C. Maxwell, On a particular case of the descent of a heavy body in a resisting medium, Camb. Dublin Math. J 9 (1854) 145–148.
  * [12] F. Auguste, J. Magnaudet, D. Fabre, Falling styles of disks, Journal of Fluid Mechanics 719 (2013) 388–405.
  * [13] M. Chrust, G. Bouchet, J. Dušek, Numerical simulation of the dynamics of freely falling discs, Physics of Fluids 25 (4) (2013) 044102.
  * [14] C. Toupoint, P. Ern, V. Roig, Kinematics and wake of freely falling cylinders at moderate reynolds numbers, Journal of Fluid Mechanics 866 (2019) 82–111.
  * [15] K. Amin, J. Mac Huang, K. J. Hu, J. Zhang, L. Ristroph, The role of shape-dependent flight stability in the origin of oriented meteorites, Proceedings of the National Academy of Sciences 116 (33) (2019) 16180–16185.
  * [16] J.-T. Kim, Y. Jin, S. Shen, A. Dash, L. P. Chamorro, Free fall of homogeneous and heterogeneous cones, Physical Review Fluids 5 (9) (2020) 093801.
  * [17] L. B. Esteban, J. Shrimpton, B. Ganapathisubramani, Edge effects on the fluttering characteristics of freely falling planar particles, Physical Review Fluids 3 (6) (2018) 064302.
  * [18] L. B. Esteban, J. Shrimpton, B. Ganapathisubramani, Three dimensional wakes of freely falling planar polygons, Experiments in Fluids 60 (7) (2019) 114.
  * [19] S. P. Kuznetsov, Plate falling in a fluid: Regular and chaotic dynamics of finite-dimensional models, Regular and Chaotic Dynamics 20 (3) (2015) 345–382.
  * [20] W. W. Willmarth, N. E. Hawk, R. L. Harvey, Steady and unsteady motions and wakes of freely falling disks, The physics of Fluids 7 (2) (1964) 197–208.
  * [21] S. B. Field, M. Klaus, M. Moore, F. Nori, Chaotic dynamics of falling disks, Nature 388 (6639) (1997) 252–254.
  * [22] H. Zhong, C. Lee, Z. Su, S. Chen, M. Zhou, J. Wu, Experimental investigation of freely falling thin disks. part 1. the flow structures and reynolds number effects on the zigzag motion, Journal of Fluid Mechanics 716 (2013) 228.
  * [23] H. Zhong, S. Chen, C. Lee, Experimental study of freely falling thin disks: Transition from planar zigzag to spiral, Physics of Fluids 23 (1) (2011) 011702\.
  * [24] C. Lee, Z. Su, H. Zhong, S. Chen, M. Zhou, J. Wu, Experimental investigation of freely falling thin disks. part 2. transition of three-dimensional motion from zigzag to spiral, Journal of Fluid Mechanics 732 (2013) 77–104.
  * [25] Y. Wang, C. Shu, C. Teo, L. Yang, An efficient immersed boundary-lattice boltzmann flux solver for simulation of 3d incompressible flows with complex geometry, Computers & Fluids 124 (2016) 54–66.
  * [26] Y. Wang, C. Shu, C. Teo, L. Yang, Numerical study on the freely falling plate: Effects of density ratio and thickness-to-length ratio, Physics of Fluids 28 (10) (2016) 103603.
  * [27] T. Kim, J. Chang, D. Kim, Free-fall dynamics of a pair of rigidly linked disks, Physics of Fluids 30 (3) (2018) 034104.
  * [28] M. Lee, S. H. Lee, D. Kim, Stabilized motion of a freely falling bristled disk, Physics of Fluids 32 (11) (2020) 113604.
  * [29] G. Kirchhoff, Ueber die bewegung eines rotationskörpers in einer flüssigkeit., Journal für die reine und angewandte Mathematik 1870 (71) (1870) 237–262.
  * [30] A. V. Borisov, I. S. Mamaev, On the motion of a heavy rigid body in an ideal fluid with circulation, Chaos: An Interdisciplinary Journal of Nonlinear Science 16 (1) (2006) 013118.
  * [31] A. V. Borisov, V. V. Kozlov, I. S. Mamaev, Asymptotic stability and associated problems of dynamics of falling rigid body, Regular and Chaotic Dynamics 12 (5) (2007) 531–565.
  * [32] Y. Tanabe, K. Kaneko, Behavior of a falling paper, Physical Review Letters 73 (10) (1994) 1372.
  * [33] L. Mahadevan, H. Aref, S. Jones, Comment on “behavior of a falling paper”, Physical Review Letters 75 (7) (1995) 1420.
  * [34] Y. Tanabe, K. Kaneko, Tanabe and kaneko reply, Physical Review Letters 75 (7) (1995) 1421.
  * [35] A. Andersen, U. Pesavento, Z. J. Wang, Analysis of transitions between fluttering, tumbling and steady descent of falling cards, Journal of Fluid Mechanics 541 (1) (2005) 91–104.
  * [36] A. Andersen, U. Pesavento, Z. J. Wang, Unsteady aerodynamics of fluttering and tumbling plates, Journal of Fluid Mechanics 541 (1) (2005) 65.
  * [37] S. Liao, On the reliability of computed chaotic solutions of non-linear differential equations, Tellus A: Dynamic Meteorology and Oceanography 61 (4) (2009) 550–564.
  * [38] S. Liao, On the numerical simulation of propagation of micro-level inherent uncertainty for chaotic dynamic systems, Chaos, Solitons & Fractals 47 (2013) 1–12.
  * [39] S. Liao, Physical limit of prediction for chaotic motion of three-body problem, Communications in Nonlinear Science and Numerical Simulation 19 (3) (2014) 601–616.
  * [40] W. Tucker, Validated numerics: a short introduction to rigorous computations, Princeton University Press, 2011.
  * [41] P. Oyanarte, Mp-a multiple precision package, Computer Physics Communications 59 (1990) 345–358.
  * [42] X. Li, S. Liao, More than six hundred new families of newtonian periodic planar collisionless three-body orbits, SCIENCE CHINA Physics, Mechanics & Astronomy 60 (12) (2017) 129511.
  * [43] X. Li, Y. Jing, S. Liao, Over a thousand new periodic orbits of a planar three-body system with unequal masses, Publications of the Astronomical Society of Japan 70 (4) (2018) 64.
  * [44] X. Li, S. Liao, Collisionless periodic orbits in the free-fall three-body problem, New Astronomy 70 (2019) 22–26.
  * [45] Z. Lin, L. Wang, S. Liao, On the origin of intrinsic randomness of rayleigh-bénard turbulence, SCIENCE CHINA Physics, Mechanics & Astronomy 60 (1) (2017) 014712.
  * [46] T. Hu, S. Liao, On the risks of using double precision in numerical simulations of spatio-temporal chaos, Journal of Computational Physics (2020) 109629.
  * [47] E. N. Lorenz, Computational chaos-a prelude to computational instability, Physica D: Nonlinear Phenomena 35 (3) (1989) 299–317.


# Imitation and Supervised Learning of
Compliance for Robotic Assembly

Devesh K. Jha, Diego Romeres, William Yerazunis and Daniel Nikovski† †All
authors are with Mitsubishi Electric Research Laboratories (MERL), Cambridge,
MA, USA 02139<EMAIL_ADDRESS>
###### Abstract

We present the design of a learning-based compliance controller for assembly
operations for industrial robots. We propose a solution within the general
setting of learning from demonstration (LfD), where a nominal trajectory is
provided through demonstration by an expert teacher. This can be used to learn
a suitable representation of the skill that can be generalized to novel
positions of one of the parts involved in the assembly, for example the hole
in a peg-in-hole (PiH) insertion task. Under the expectation that this novel
position might not be entirely accurately estimated by a vision or other
sensing system, the robot will need to further modify the generated trajectory
in response to force readings measured by means of a force-torque (F/T) sensor
mounted at the wrist of the robot or another suitable location. Under the
assumption of constant velocity of traversing the reference trajectory during
assembly, we propose a novel accommodation force controller that allows the
robot to safely explore different contact configurations. The data collected
using this controller is used to train a Gaussian process model to predict the
misalignment in the position of the peg with respect to the target hole. We
show that the proposed learning-based approach can correct various contact
configurations caused by misalignment between the assembled parts in a PiH
task, achieving high success rate during insertion. We show results using an
industrial manipulator arm, and demonstrate that the proposed method can
perform adaptive insertion using force feedback from the trained machine
learning models.

##  I Introduction

Figure 1: Experimental setup with a Mitsubishi Electric Factory Automation
(MELFA) RV-$5$AS-D Assista $6$ DoF manipulator arm in a possible contact
configuration with the hole environment. The diameter of the peg is
approximately $20$ mm and the tolerance is approximately $1.5$ mm. The figure
also shows the industrial vision system that was used in the experiments for
detecting the hole.

Mating parts under tight tolerances (see Figure 1) is one of the most
important operations in robotic assembly. Over the last several decades,
robots have become very precise in performing pick and place operations.
However, complications arise when the positions of the parts involved in
assembly vary between repetitions of the operation. This can happen in cases
where the parts are deposited onto a surface by a feeder, and each time end up
in a different position; also, in a case where the parts arrive on a moving
conveyor belt. In such cases, industrial vision cameras can be used to
determine the parts’ position. However, the determined position of the parts
by the industrial vision cameras is usually fairly inaccurate, often by
several millimeters. This is significantly more than the required tolerances
during assembly. For this reason, standard compliance (e.g., stiffness)
controllers usually fail to perform the operation successfully. Besides, even
if the robotic device knows the exact position of the parts, its path (via
points) still needs to be modified in order to accommodate the variation in
the position of the parts. In practice, such a modification is performed by
means of dedicated software that takes as input the changed position of the
parts and outputs a new path for the robot. However, developing such dedicated
software is typically very difficult and laborious, and is currently one of
the main contributors to the high cost of deploying and re-tooling robotic
devices for new assembly operations.

Therefore, there is a strong industrial need for methods for creating robot
controllers that can accommodate varying positions of assembly parts with
minimal or no programming. That is, the robot controller should be able to
generalize over the positions of the parts. Arguably, generalization is one of
the hallmarks of intelligence, be it human or artificial, and endowing robots
with the ability to generalize over the starting conditions of assembly tasks
would amount to a qualitatively different type of intelligent assembly
functionality. Such a functionality would match the ability of human workers
to learn a new assembly skill and apply it to many new situations.

In this paper, we present a two-fold approach for assembly so that the agent
can generalize to a different initial condition and accommodate uncertainty in
the final goal position. A demonstration is used to learn the approach
movement of the robot from any initial position to a novel target position.
Then, a supervised learning method is used to learn a force-feedback policy to
accommodate and correct for the uncertainty in the estimated target position.
The purpose of the force-feedback policy is to interpret the contact
information and correct the possible misalignment. A commonly used
prototypical robot assembly operation is the insertion of a peg into a hole
(see Figure 1), and we will use this task as an example when describing and
evaluating the algorithm. We present a method that allows a robot equipped
with a camera and a force/torque sensor mounted on its wrist to learn this
operation, with minimal or no programming, even when the estimates of the
positions of the parts (returned by the camera) are not very precise. We
assume that a human instructor can provide at least one demonstration of a
successful operation, using the robot itself.

Contributions: Figure 2 provides an abstract idea of the proposed method for
peg-in-hole assembly. To enable safe sustained contacts with the environment
during insertion attempts, we present the design of an accommodation
controller. We use our novel accommodation controller as the low-level
controller during exploration as well as operation. The proposed method uses
this accommodation controller to maintain safe contact during the data
collection process. A machine learning model then learns a mapping from the
(quasi-)static contact wrench to the desired correction to perform insertion.
In particular, the proposed paper has the following contributions:

  * •

We present a design of a learning-based adaptive controller to perform
insertion with incorrect knowledge of the goal location, which can achieve
near-perfect insertion performance. This is demonstrated on a physical robotic
system with industrial-grade vision systems.

  * •

We present the design and analysis of a stable accommodation controller that
allows maintaining contacts with the external environment while ensuring
contact forces to always stay within desired bounds when following a dynamic
reference trajectory.

Figure 2: The proposed method for design of a force controller-based learning
method for part assembly. The proposed method uses an accommodation controller
as a low-level controller for safe interaction during data collection. The
data collected during exploration is used to train regression models to
correct for the misalignment between the peg and the hole. This machine
learning model-based controller is used at the higher level to correct the
misalignment, so as to perform insertion.

##  II Related Work

Research on robotic assembly has a long history, and a lot of this research
has used PiH insertion as a prototypical task. The main reason this task is
difficult to automate is that when the peg and hole are in contact, even a
small misalignment between their correct positions would cause a controller
that tracks a desired trajectory defined in terms of robot positions or
applied torques at the joints to result in very large contact forces that
might break the robot and/or the parts. What is needed is a method that
adjusts the motion of the parts in response to contact forces. Such techniques
are also known as adaptive assembly strategies (AAS) ([1, 2, 3, 4, 5, 6, 7, 8,
9, 10]).

It is generally desirable that a controller be able to make an active
interpretation of force/torque (F/T) signals to generate corrective actions.
Several ideas have been proposed in the literature to design controllers to
compensate for larger misalignment during insertion. The idea in many of these
controllers is to use a mapping from the contact wrench measured at the wrist
of the robot to a correction in the motion for the robot. For PiH problems, in
rare instances, it is possible to compute this mapping analytically — for
example, when the axes of a round peg and hole are aligned perfectly, the
point around which moments measured by the F/T sensor lies on axis of the peg,
and at least some overlap between the peg and hole exists, it is possible to
calculate the direction to the center of the hole using a closed-form
expression ([11]). However, this solution requires careful calibration of the
sensing system, as changing the location of the peg with respect to the end
tool of the robot and its F/T sensor will affect the sensed force signature in
a given contact configuration and render the analytical formulae incorrect.

Another approach for obtaining the mapping is to use machine learning
algorithms, such as supervised learning or reinforcement learning. An early
method for programmed compliance proposed in [12] used linear mappings in the
form of linear admittance or accommodation matrices, much like earlier
compliance (e.g. stiffness) controllers. Later, it was demonstrated that for
insertion tasks, the set of specified appropriate examples of corrective
action cannot be represented by a linear mapping, in general, and the use of
nonlinear regression models such as neural networks was proposed ([13, 14]).
This further significantly expanded the representation power of the learned
controller, but the question remained how to come up with this set of examples
of correct response to contact forces. To this end, a reinforcement learning
(RL) solution based on a process of trial and error was proposed in [4],
leading to a very impressive solution of a PiH problem with tight tolerances
that were orders of magnitude tighter than the accuracy (repeatability) of the
used robot. However, this method still assumed accurate knowledge of the goal
state, so it is not directly applicable to the problem of insertion with
uncertain goal state.

Following the seminal work in [4], there have been numerous works using
machine learning models for designing adaptive controllers. End-to-end
visuomotor training of robots by means of deep RL was demonstrated in [6, 8].
More recently, there has been some work on using tactile sensors for designing
feedback policies to perform insertion [15, 7]. The use of RL with recurrent
neural networks for insertion problems was demonstrated recently in [9].
However, although being a very powerful general method for learning control
policies, RL is also known for its unfavorable sample complexity, making it a
less attractive choice for learning on actual physical systems. Compared with
all these works, we present an approach that uses a combination of LfD and a
supervised learning-based method for local adjustment to the LfD trajectory
that is very sample-efficient in its interaction with the physical equipment.
Furthermore, we also present the use of a special-purpose accommodation
controller to safely interact with the environment during exploration for data
collection and actual insertion operations.

##  III Problem Statement

Peg-in-hole insertion is an entire class of problems whose difficulty can vary
widely depending on the type of misalignment present (positional only or
positional plus directional), its magnitude (whether the peg is already
partially in the hole, or if not, whether the peg and hole at least overlap
some), the tolerance of the insertion operation, whether the parts are
chamfered, what kind of sensors are used (e.g., visual, force, tactile), what
their relative accuracy is, etc. ([16]). An additional consideration is how
the robot is controlled (full joint torque control or position control only).

In this paper, we are considering a version of the problem that corresponds to
fairly typical circumstances often encountered in practice. We are designing
controllers to handle inaccuracies in the goal position that are much larger
than the tolerance of insertion, but still smaller than the size of the parts
being inserted. This situation arises, for example, when we use a camera to
estimate the position of the hole. The positional estimation error of modern
cameras for observed objects in the world frame tends to be on the order of
several millimeters, which is much larger than typical insertion tolerances in
assembly (often sub-millimeter), but still less than the radius of the parts
being inserted. This assumption entails that if the robot executes in open
loop a motion to where it believes the goal state is, even if the insertion
does not succeed and the peg collides with the rim of the hole, the resulting
contact configuration will have the peg and hole overlap to at least some
degree, allowing the measured force signature to be used to correct the
estimate of the goal position. (If this condition is not true, and there is no
overlap, a single F/T reading cannot disambiguate the contact configuration,
and the controller must resort to integrating multiple readings, or perform
search; this is a much more difficult version of the problem that we are not
exploring in this paper, and is usually addressed by very different kinds of
algorithms in the research community.) In addition, in this paper, we are
assuming that there is no (or very minimal) directional misalignment between
the peg and the hole, which corresponds to the case when the hole moves on a
plane (for example, a work surface in a factory), changing only its position,
but not the angle of its axis with respect to the plane.

In terms of contact sensors, we are considering a force/torque sensor mounted
at the wrist of the robot. Furthermore, we are using industrial manipulators
that are controlled only in (possibly compliant) position-control mode, that
is, we do not assume the availability of direct joint-torque-control mode.
This version of the problem matches almost exactly the one used in recent work
reported in [9], as well as earlier work on robotics assembly [11].

##  IV Proposed Method

The method proposed in this paper uses two different machine learning methods,
as described below. One of these methods is to learn a nominal trajectory for
insertion in the form of a dynamical movement primitive (DMP, [17]), and
adjust it appropriately for novel start and goal states. The second method is
based on a form of active supervised learning that collects training examples
relating various displacements of the peg from its correct position for
insertion and the F/T observations measured as a result, by means of active
experimentation with the peg and hole.

Learning starts by collecting one or more demonstrations of the operation by
an expert. The result of the demonstration(s) is one or more trajectories in
Cartesian space of the end tool of the robot of the form $\mathbf{x}(t)$,
$t\in[0,T]$, $\mathbf{x}(t)\in\mathbb{R}^{3}$. Pose of the parts can be
estimated by means of one or more cameras. Pose estimation by means of
computer vision is an active area of research ([18]), and various methods
provide different trade-offs between cost, speed, and accuracy. We assume that
there exists a pose estimation module which is properly calibrated with
respect to the start and end of the demonstration trajectory $\mathbf{x}(t)$,
that is, $\mathbf{x}_{0}=\mathbf{x}(0)$ and $g=\mathbf{x}(T)$ are indeed the
estimates returned by the camera.

After one or more trajectories $\mathbf{x}(t)$ have been recorded, a dynamic
movement primitive (DMP) learning algorithm is used to learn a separate DMP
for each of the components of $\mathbf{x}(t)$ ([17]). The resulting set of
DMPs can generate a new desired trajectory $\mathbf{x}_{r}(t)$, given a new
goal position $g_{d}$, by integrating the DMP forward in time. This kind of
generalization over the goal state is an essential function of DMPs. If the
new desired goal state $g_{d}$ was indeed the correct one, then following the
new desired trajectory $\mathbf{x}_{r}(t)$ would presumably reach the new goal
state. However, the estimate of the desired new goal state $g_{d}$ is usually
incorrect, because this estimate depends on the new position of the hole, as
described below. Any error in the estimate of the position of the hole
translates into an error in $g_{d}$. For this reason, the control law executed
by the robot needs to provide compliance. We use a control law of the form
$\mathbf{x}_{c}(t)=\mathbf{x}_{r}(t)+H(f)$, where $f$ is the vector of sensed
forces and torques by the FT sensor, and $H$ is a non-linear mapping that
produces corrections to the position of the robot. Here, $H$ is a general non-
linear mapping, which is learnt using a self-supervised learning method [19,
20, 12, 13].

The purpose of the next stage of the algorithm is to learn the mapping $H$
through self-experimentation. To this end, while the hole is still in its
original position, variations $d(t)$ are introduced to one of the original
trajectories $x(t)$ demonstrated by a human operator that successfully
completed the operation. To learn this mapping, we collect force signature
data of the form $(f_{i},d_{i})$, where $f_{i}$ is the force wrench
experienced by introducing a variation $d_{i}$. However, based on the
understanding that the mapping from displacements to forces is typically many-
to-one (multiple displacements would sometimes result in the same force),
following the analysis in [21], the inverse mapping would be one-to-many, that
is, not a function that can be learned by means of machine learning. However,
it can be realized that the exact magnitude of the displacement does not need
to be recovered for successful corrective action, and furthermore, multiple
displacements can generate the same forces only if the sign of these
displacements is the same, as long as the magnitude of the displacement does
not exceed the radius $R$ of the object being inserted. Based on this
realization, a supervised machine learning algorithm can be used to learn the
mapping $sign(d_{i})=H_{0}(f_{i})$, for all examples $i=1,N$ such that
$|d_{i}|\leq R$. Such a learned machine learning model is then used to a
design a corrective controller to correct the misalignment between the mating
parts.

A practical consideration about the exploration method is how to do it safely.
If the robot is commanded to follow the modified trajectory in standard
position control mode, the resulting collision will either damage the robot or
the parts, or result in a fault due to excessive force. This means that the
robot needs to be controlled in a compliant control mode. One accommodation
controller design for safe exploration and exploitation of the learned control
policy based on these principles and suitable for execution on industrial
robot arms is descried in the next section.

Figure 3: This plot shows the convergence of contact wrench obtained using the
proposed GAC during a contact configuration during an insertion attempt. The
reference trajectory advances linearly in the insertion direction at constant
velocity even after actual movement of the robot stops due to collision of the
peg with the rim of the hole. For the same reference trajectory, an admittance
controller would have produced ever increasing forces and moments.

##  V Force Controller Design for Safe Interaction

Figure 4: Block Diagram of the generalized accommodation controller that is
designed for the data collection process. $K_{a}$: Accommodation Matrix,
$K_{s}$: Stiffness Gains of the low-level stiffness controller.

While there has been significant work in the literature to perform active
control using the contact forces during insertion-type assembly operations,
design of controllers to ensure safe, sustained contacts has received less
attention. In this section, we present the design and analysis of a force
controller that ensures safe interaction by maintaining contact forces within
safe bounds. In particular, we try to answer and explain the following in this
section:

  1. 1.

Design of a Generalized Accommodation Controller (GAC) for safe exploration
during assembly attempts.

  2. 2.

Evaluate the convergence characteristics of the proposed controller.

It is important to note that the design of GAC is critical to the design of
the ACC, as this allows us to collect data in quasi-steady state by making the
contact forces converge to a value that is representative of the contact
configuration, and not of the applied effort by the robot, thus allowing
interpretation of the phenomena in quasi-static states (instead of in dynamic
states). It is also desirable to have such a controller to make sure that the
interaction force between the robot and its environment stays within some
desirable safe bounds. We first explain the design of this controller using
the stock stiffness controller provided in general industrial robots.

###  V-A Generalized Accommodation Controller Design

The idea of the proposed accommodation controller is presented in Figure 4. As
could be seen in the block diagram in Figure 4, the accommodation controller
manipulates the reference trajectory using force feedback. In particular, the
accommodation controller uses the following feedback law to manipulate the
reference trajectory of the robot. Let us denote the discrete-time reference
trajectory by $x_{r}[k]$, the trajectory commanded to the low-level position
controller by $x_{c}[k]$, the experienced forces by $f[k]$, the measured
position as $x[k]$, at any instant $t_{k}$. Let us denote the stiffness
constant of the compliant position controller by $K_{s}$ and the accommodation
matrix for the force feedback as $K_{a}$. For simplicity, we consider a
diagonal matrix $K_{a}$. The commanded trajectory sent to the robot is then
computed using the following update rule:

$x_{c}[k]=x_{c}[k-1]+\Delta x_{r}[k]+\sum_{i=0}^{k}\gamma^{k-i}K_{a}f[k]$ (1)

where $\gamma\in(0,1)$ is a discounting parameter for computing the integral
error, and $\Delta x_{r}[k]=x_{r}[k]-x_{r}[k-1]$ are desired position
increments computed from the reference trajectory (normally itself computed by
the DMP module based on its estimate of the goal position). Note that the
force experienced during interaction is governed by the stiffness constant of
the robot controller and is given by $f[k]=K_{s}(x[k]-x_{c}[k])$. We assume
that such a controller is available in the control software of the robot; if
not, it can be implemented by interpreting the contact forces accordingly
([22]).

###  V-B Convergence Behavior of the Accommodation Controller

In this section, we briefly describe the convergence behavior of the
accommodation controller that was described in the previous section. It has
the useful property that even when the reference trajectory advances
indefinitely, the contact force converges to a constant, as long as the
reference velocity is constant. Because the reference velocity (that is, the
velocity at which the desired reference trajectory $x_{r}[k]$ is traversed) is
entirely under the control of the designer of the controller, this effectively
means that the contact force is under their control, too, and can be set to a
desired safe value. This is in contrast to the behavior of a standard
compliance controller, such as a stiffness controller, where contact forces
would increase indefinitely if the reference trajectory keeps advancing for
any reason.

Convergence of the GAC can be understood by considering that by using suitable
values for the accommodation matrix $K_{a}$ and the discounting factor
$\gamma$, the contact forces can be regulated so that both they and the
commanded trajectory reach a constant value, when the position change in the
reference trajectory is kept constant, that is, $\Delta x_{r}[k]=v$ for some
constant velocity $v$ along the intended direction of insertion. To see this,
the commanded position will stop changing ($x_{c}[k]=x_{c}[k-1])$ when both
update terms to it cancel each other, that is, $v=\Delta
x_{r}[k]=-\sum\limits_{i=0}^{k}\gamma^{k-i}K_{a}f[k]$. Under this condition,
the commanded trajectory converges to $\tilde{x}_{c}=x_{c}[k]$. Then, if the
movement of the robot has also stopped at some equilibrium position
$\tilde{x}$, under the operation of the stiffness position controller of the
robot with a constant commanded position $\tilde{x}_{c}$, the steady-state
force experienced by the robot is given by
$\tilde{f}=K_{s}(\tilde{x}-\tilde{x}_{c})$. Then, at steady state,
$v=K_{a}/(1-\gamma)\tilde{f}$, or $\tilde{f}=v(1-\gamma)/K_{a}$. This
demonstrates that the contact force at equilibrium can be reduced by choosing
a lower speed of traversal of the reference trajectory, or a higher
discounting factor $\gamma$, or a lower value of the accommodation constant
$K_{a}$. The convergence behavior of the forces is also shown in Figure 3. As
could be seen in the plots, the force and moments experienced by the robot
converge to a steady-state value. The speed of convergence of the system is
governed by the discounting parameter $\gamma$. We observed that the
convergence is faster with higher values of $\gamma$. In practice, any
parameter $\gamma\in[0.3,0.65]$ could achieve satisfactory convergence. Values
of $\gamma$ higher than $0.65$ might sometimes lead to oscillations in the
system.

In the rest of the paper, we use this controller as the trajectory-following
controller which is used when following a reference trajectory obtained from a
higher level controller working towards achieving successful insertion for
assembly. This is described in detail in next section.

##  VI Experimental Results

Figure 5: Variation of contact wrench obtained from contact formations with
different misalignment w.r.t. the hole along the $y$ axis, for a constant
value of $x$. This plot shows that, in practice, we obtain a non-linear
dependence between the misalignment and the wrench experienced during the
contact formation.

In this section, we present results of learning the corrective action in
different contact configurations. In particular, we try to answer the
following questions:

  1. 1.

Can the proposed method resolve the direction of the error as well as its
magnitude?

  2. 2.

What is the range of misalignment that can be corrected using such a learned
model?

  3. 3.

How well does the proposed controller perform on insertion tasks with
inaccurate knowledge of goal the state?

In the rest of this section, we try to answer these questions.

###  VI-A Experiment Details

We use a Mitsubishi Electric Factory Automation (MELFA) RV-$5$AS-D Assista
$6$-DoF arm (see Figure 1) for the experiments. The robot has a pose
repeatability of $\pm 0.03$mm. The robot is equipped with Mitsubishi Electric
F/T sensor $1$F-FS$001$-W$200$ (see Figure 1). The peg has a diameter of $20$
mm and the hole has a diameter of $22$ mm. We use the underlying stiffness
controller of the robot to design the corresponding force controller. We use
the stiffness controller with $K_{z}=2$ N/mm and $K_{x}=K_{y}=10$ N/mm as the
stiffness constant during all the experiments. We use the GAC only in the axis
of approach during insertion operation but one can implement it for all the
six axes.

During the data collection procedure, error is introduced in the $x$ and $y$
position of the robot end-effector from the correct target state. This error
is sampled from a uniform distribution in the range $[-R,R]$ around the
circumference of the hole, where $R$ is the radius of the peg. Note that this
is sufficiently larger than the error introduced in previous work, e.g., [7,
9]. For example, in [9] the baseline success rate is already $40\%$ whereas in
the present study, it is $0\%$. Furthermore, whereas the method in [9] studied
insertion in a slot-like environment (i.e., misalignment is only along one
dimension), we study a more complex problem of peg insertion where error is
introduced simultaneously in two dimensions. The robot is moved along the
known approach direction with this error and using the force controller
described in Section V (see Figure 2). The controller maintains contact with
the surface after the initial contact, and the force data for the entire
interaction is stored. The total length of the trajectory is $30$ seconds, and
the data is collected at $280$Hz, which is the default control rate of the
robot.

###  VI-B Dependence of Contact Wrench on Misalignment

In this section, we analyze the behavior of the contact wrench experienced by
the peg during contact formation during an insertion attempt with varying
amount of misalignment. The variation of the contact wrench with varying
amount of misalignment along the $y$ axis, at a constant value of $x$, is
shown in Figure 5. These plots show the (quasi-)steady state values of the
contact wrench for a contact formation, that is, measured at the end of the
insertion attempt. There are several things to infer from these plots shown in
Figure 5. As seen in the plot of $F_{z}$, the force experienced in the
vertical direction is constant for all values of misalignment. This force
creates a moment about the x-axis which can be observed in the plot for
$M_{x}$. This is also explained by a planar analytical model for the static
equilibrium of the peg in the contact configuration. We also observe a moment
about the y-axis of the peg which can be observed in the plot for $M_{y}$. We
also observe that there is a non-linear dependence between the moments
experienced by the peg and the amount of misalignment. While a planar
analytical model would predict that the force and moments change linearly
during these sustained contacts, we see here that there is significant non-
linearity that is observed, and thus it provides the intuition for using a
non-linear machine learning model for adjusting the misalignment based on the
force readings obtained from the F/T sensor.

###  VI-C Predictive Model for Correcting Misalignment

In this section, we describe learning a model to predict the actual
misalignment between the hole and the peg. We train a models to predict the
misalignment and try to understand the accuracy and degree to which the
misalignment can be corrected using a supervised machine learning model. The
experimental setup that we use for data collection is shown in Figure 1 and
was described in the previous section. We use Gaussian processes for learning
the predictive model for misalignment. Gaussian processes are data efficient,
and have been shown to perform well in a wide range of problems [23, 24, 25].
The Gaussian process regression is trained using the data collected using the
accommodation controller described in Section V. As was shown earlier in
Figure 3, the system reaches (quasi) steady state using the force controller.
The Gaussian process regression is trained using the (quasi-)steady-state
force signature obtained by the force controller. Thus the input to the
regression model is a $6$ dimensional vector which is used to predict the
misalignment along the $x$ and $y$ axis. We train two GPR models, one for each
axis which are then used to predict the errors in each of the axis. We collect
data from a total of $1,200$ trajectories with randomly sampled misalignment
around the circumference of the hole. We train two different regression models
to predict error along the x-axis and the y-axis. We use a sum-kernel with the
RBF kernel and the white kernel to prevent overfitting of data. Figure 6 shows
the prediction of the regression models for misalignment along the two axis.
These plots show the mean RMSE predicted from $50$ different GP models trained
and tested over different subsets of data from the total collected data. We
also observe that the models have consistently better predictive behavior in
the $y$ axis compared to the $x$ axis (this is an artifact of the contact
configurations and not always true in general). We also train Gaussian process
classifiers to predict the direction of error using a similar sum-kernel to
analyze the separation of data in predicting the direction. The classification
results are listed in Table I which reflects almost perfect direction
recovery.

TABLE I: Classification accuracy in prediction of direction of misalignment along X and Y axis using Gaussian process classifiers. Axis | Classification Accuracy (higher is better)
---|---
X | 0.97
Y | 1.0
Figure 6: Mean RMSE error (lower is better) for misalignment prediction along
the $x$ and $y$ axis with increasing values of misalignment. As we can see
from the plots, the RMSE monotonically increases with increasing misalignment.
Figure 7: The hole surface is moved in the plane in an approximate rectangular
region of $12$ cm $\times$ $10$ cm. This figure shows the approximate region
in the plane where the hole could be detected by the external camera (see
Figure 1). Note that the images show the standard output of an industrial-
grade camera. The frame shown is the robot frame of reference (see also Figure
1).

###  VI-D Performance of the ML-based Controller

In this section, we evaluate performance of the learning-based controller in
correcting the misalignment during insertion experiments. We use the Gaussian
process classification model described in the previous section to correct for
the positional error using force feedback. In particular, we design a
controller that uses a reference trajectory and the learned correction model
to successfully perform insertion using force feedback. As shown in Figure 2,
we use the accommodation controller at the lower level, and then use the
predictions of the GP model to correct the reference trajectory using the
prediction. We perform two sets of experiments. In the first set of
experiments, the hole base is manually moved around in the plane (see Figure
7) inside a rectangular region of approximately $12$ cm $\times$ $10$ cm to
generate $50$ novel goal locations. The hole location is detected by the
external camera. This experiment tries to test the generalization of the
controller to novel contact formations as the hole base is moved around. We
observe that the camera error is biased in a particular direction, and thus
the robot always ends up approaching the correction from a particular
direction. Thus, we also design a second set of experiments where artificial
errors are added to known location by sampling from the uniform distribution.
This distribution is the same as during training time. During the force
feedback correction, we use the predictions made by the Gaussian process
classifier to move the robot in the direction predicted by the learned model.
A unit step of $2$mm is used to move the robot along each axis. It can be
adaptively changed based on the successive steps predicted by the model. The
robot uses the contact wrench at the end of the contact trajectory as input to
the learned Gaussian process model, uses the predictions from the model to
make corrections in the $x$ and $y$ axis as predicted by the classifier model,
makes contact with the hole surface again, and repeats the movement till
either it achieves successful insertion or declares failure. A failure is
declared based on either the total number of attempts exceeds the maximum
number of attempts ($10$) or the controller diverging from the desired goal
state more than $15$ mm. We implement two different controllers during test –
the first being the DMP trajectory (LfD) computed using the initial state and
the novel goal location. Note that this is computed using a single successful
demonstration directly provided on the robot using a joystick. The second
controller is the proposed DMP trajectory augmented with the additional force-
feedback controller (LfD+SL). The test results are listed in Table II. As can
be seen, the proposed method achieves very high success rate during test
times, while the DMP-based controller fails in all test condition in the final
insertion.

TABLE II: Success rate of different controllers Experiment type | LfD+ SL | LfD
---|---|---
Moving the hole | 48/50 | 0/50
Error added at known location | 79/80 | 0/80

##  VII Conclusions and Future Work

Insertion is the most important operation of robotic assembly systems.
Consequently, it has been widely studied to create robotic systems which can
perform autonomous insertion using sensor-based feedback. However, the problem
still remains open in several aspects, and there is no single technique to
perform generalized insertion. In this paper, we presented a learning method
for solving the problem of peg insertion during robotic assembly. Most of
robotic assembly problems involve sustained contact phenomena, which makes it
harder to ensure reliable and safe operation. We presented the design and
analysis of an adaptive force controller that allows us to maintain sustained
contacts with limited contact force in the presence of dynamic reference
trajectory. We show that our proposed controller results in (quasi-)steady-
state behavior by modifying the reference trajectory as the robot comes in
contact with the external environment. We use this controller for safe
exploration to collect data during the insertion process. This is used to
learn a predictive model that can correct the misalignment during insertion.
In our experiments, we show that we achieve approximately $98.75\%$ success
rate with error sampled in the range of $[-R,R]$ in both the $x$ and $y$ axis,
where $R$ is the radius of the inserted part. In contrast to the use of
closed-form expressions used for solving similar problems in [11], our method
does not need careful calibration, and works even when the axis of the peg is
not aligned with the F/T sensor. Similarly, in comparison to the method
recently reported in [9], our method uses supervised learning instead of
reinforcement learning, thus avoiding the need for prolonged random
exploration.

Our method uses the quasi steady-state force data for learning the corrective
compliance. However, converging to this quasi-steady state does take some
time, which probably makes the insertion operation slower than it could be. In
the future, we would like to use contact force readings from times before full
convergence, hopefully speeding up the insertion process. In addition, we plan
to study and analyze contact configurations with directional (rotational)
misalignment between the peg and the hole.

## References

  * [1] J. L. Nevins and D. E. Whitney, “Research on Advanced Assembly Automation,” _IEEE Computer_ , vol. 10, no. 12, pp. 24–38, 1977.
  * [2] M. T. Mason, “Compliance and Force Control for Computer Controlled Manipulators,” _IEEE Transactions on Systems, Man and Cybernetics_ , vol. 11, no. 6, pp. 418–432, 1981.
  * [3] H. Inoue, “Force Feedback in Precise Assembly Tasks,” 8 1974.
  * [4] V. Gullapalli, J. A. Franklin, and H. Benbrahim, “Acquiring Robot Skills via Reinforcement Learning,” _IEEE Control Systems_ , vol. 14, no. 1, pp. 13–24, 1994.
  * [5] F. J. Abu-Dakka, B. Nemec, J. A. Jørgensen, T. R. Savarimuthu, N. Krüger, and A. Ude, “Adaptation of manipulation skills in physical contact with the environment to reference force profiles,” _Autonomous Robots_ , vol. 39, no. 2, pp. 199–217, 2015.
  * [6] S. Levine, C. Finn, T. Darrell, and P. Abbeel, “End-to-end training of deep visuomotor policies,” _Journal of Machine Learning Research_ , vol. 17, pp. 1–40, 2016.
  * [7] S. Dong, D. Jha, D. Romeres, S. Kim, D. Nikovski, and A. Rodriguez, “Tactile-RL for insertion: Generalization to objects of unknown geometry,” in _2021 IEEE International Conference on Robotics and Automation (ICRA)_ , 2021.
  * [8] M. A. Lee, Y. Zhu, P. Zachares, M. Tan, K. Srinivasan, S. Savarese, L. Fei-Fei, A. Garg, and J. Bohg, “Making sense of vision and touch: Learning multimodal representations for contact-rich tasks,” _IEEE Transactions on Robotics_ , vol. 36, no. 3, pp. 582–596, 2020.
  * [9] P. Kulkarni, J. Kober, R. Babuška, and C. D. Santina, “Learning Assembly Tasks in a Few Minutes by Combining Impedance Control and Residual Recurrent Reinforcement Learning,” _Advanced Intelligent Systems_ , p. 2100095, 9 2021.
  * [10] Y. Liu, D. Romeres, D. K. Jha, and D. Nikovski, “Understanding multi-modal perception using behavioral cloning for peg-in-a-hole insertion tasks,” _CoRR_ , vol. abs/2007.11646, 2020.
  * [11] S. Gottschlich and Kak AC, “A dynamic approach to high-precision parts mating,” _IEEE Transactions on Systems, Man, and Cybernetics_ , vol. 19, pp. 797–810, 1989.
  * [12] M. A. Peshkin, “Programmed Compliance for Error Corrective Assembly,” _IEEE Transactions on Robotics and Automation_ , vol. 6, no. 4, pp. 473–482, 1990.
  * [13] H. Asada, “Teaching and learning of compliance using neural nets: Representation and generation of nonlinear compliance,” in _International Conference on Robotics and Automation_. IEEE, 1990, pp. 1237–1244.
  * [14] ——, “Representation and Learning of Nonlinear Compliance Using Neural Nets,” _IEEE Transactions on Robotics and Automation_ , vol. 9, no. 6, pp. 863–867, 1993.
  * [15] S. Dong and A. Rodriguez, “Tactile-Based Insertion for Dense Box-Packing,” in _IEEE International Conference on Intelligent Robots and Systems_. Institute of Electrical and Electronics Engineers Inc., 11 2019, pp. 7953–7960.
  * [16] D. E. Whitney, “Historical Perspective and State of the Art in Robot Force Control,” _The International Journal of Robotics Research_ , vol. 6, no. 1, pp. 3–14, 1987.
  * [17] A. J. Ijspeert, J. Nakanishi, H. Hoffmann, P. Pastor, and S. Schaal, “Dynamical movement primitives: Learning attractor models formotor behaviors,” _Neural Computation_ , vol. 25, no. 2, pp. 328–373, 1 2013\.
  * [18] T. Hodan, R. Kouskouridas, T.-K. Kim, F. Tombari, K. Bekris, B. Drost, T. Groueix, K. Walas, V. Lepetit, A. Leonardis, C. Steger, F. Michel, C. Sahin, C. Rother, and J. Matas, “A summary of the 4th International Workshop on Recovering 6D Object Pose,” _Lecture Notes in Computer Science (including subseries Lecture Notes in Artificial Intelligence and Lecture Notes in Bioinformatics)_ , vol. 11129 LNCS, pp. 589–600, 10 2018.
  * [19] R. Groome, _Force feedback steering of a teleoperator system (MSc Thesis)_ , 1972.
  * [20] D. Whitney, “Force feedback control of manipulator fine motions,” 1977.
  * [21] W. S. Newman, Y. Zhao, and Y. H. Pao, “Interpretation of force and moment signals for compliant peg-in-hole assembly,” in _IEEE International Conference on Robotics and Automation_ , vol. 1, 2001, pp. 571–576.
  * [22] K. Lynch and F. Park, _Modern Robotics_. Cambridge University Press, 2017.
  * [23] K. Ota, D. K. Jha, D. Romeres, J. van Baar, K. A. Smith, T. Semitsu, T. Oiki, A. Sullivan, D. Nikovski, and J. B. Tenenbaum, “Data-efficient learning for complex and real-time physical problem solving using augmented simulation,” _IEEE Robotics and Automation Letters_ , vol. 6, no. 2, pp. 4241–4248, 2021\.
  * [24] D. Romeres, D. K. Jha, W. Yerazunis, D. Nikovski, and H. A. Dau, “Anomaly detection for insertion tasks in robotic assembly using gaussian process models,” in _2019 18th European Control Conference (ECC)_. IEEE, 2019, pp. 1017–1022.
  * [25] D. Romeres, D. K. Jha, A. DallaLibera, B. Yerazunis, and D. Nikovski, “Semiparametrical gaussian processes learning of forward dynamical models for navigating in a circular maze,” in _2019 International Conference on Robotics and Automation (ICRA)_. IEEE, 2019, pp. 3195–3202.