# Forward Sensitivity Analysis of the FitzHugh-Nagumo System: Parameter
Estimation

Shady E. Ahmed
School of Mechanical & Aerospace Engineering,
Oklahoma State University,
Stillwater, OK 74078, USA.
<EMAIL_ADDRESS>&Omer San
School of Mechanical & Aerospace Engineering,
Oklahoma State University,
Stillwater, OK 74078, USA.
<EMAIL_ADDRESS>&Sivaramakrishnan Lakshmivarahan
School of Computer Science,
University of Oklahoma,
Norman, Oklahoma - 73019, USA.
<EMAIL_ADDRESS>

###### Abstract

The FitzHugh-Nagumo (FHN) model, from computational neuroscience, has
attracted attention in nonlinear dynamics studies as it describes the behavior
of excitable systems and exhibits interesting bifurcation properties. The
accurate estimation of the model parameters is vital to understand how the
solution trajectory evolves in time. To this end, we provide a forward
sensitivity method (FSM) approach to quantify the main model parameters using
sparse measurement data. FSM constitutes a variational data assimilation
technique which integrates model sensitivities into the process of fitting the
model to the observations. We analyse the applicability of FSM to update the
FHN model parameters and predict its dynamical characteristics. Furthermore,
we highlight a few guidelines for observations placement to control the shape
of the cost functional and improve the parameter inference iterations.

_K_ eywords Forward sensitivity, parameter estimation, FitzHugh-Nagumo model,
data assimilation.

##  1 Introduction

Dynamical systems are ubiquitous around us and in every scientific discipline.
Examples from physical sciences include atmospheric and oceanic flows, heat
and mass transfer, the behavior of moving objects (e.g., cars, ships,
airplanes, rockets, pendulums, etc.), chemical reactions, and signal
transmission. In social sciences, the population increase and distribution,
human interactions, and cultural developments over centuries have been
following interesting dynamical patterns. Researchers and practitioners in
life sciences have also found that the application of dynamical systems
theories to the bodies, organs, and cells yields significant advancement in
our understanding and treatment of the body. In neurosciences, the
understanding of brain performance and response to external stimulus has been
critical for epilepsy prevention and treatment. Several dynamical models have
been historically proposed and investigated to study and analyze the neuronal
activity (e.g., see [1]). The FitzHugh-Nagumo (FHN) equations [2, 3, 4]
represent one of the very popular and simple models in the study of neuro-
physiology. In addition to its utility for the modeling of biological
behavior, it is considered a prototypical model in the study of nonlinear
dynamics due to its interesting characteristics such as the bifurcation
properties [5].

The two-equation FHN model, describing neuronal spike discharges, can be
defined as

$\displaystyle\dfrac{\mathrm{d}v}{\mathrm{d}t}$
$\displaystyle=v-\dfrac{1}{3}v^{3}-w+I,$ (1)
$\displaystyle\tau\dfrac{\mathrm{d}w}{\mathrm{d}t}$ $\displaystyle=v+a-bw,$
(2)

where $v$ defines the membrane potential, $w$ stands for a recovery variable,
and $\tau$ is the time scale. $I$ represents the external input current, while
$a$ and $b$ are controlling parameters. The FHN model might appear in various
forms, which can be related to Eqs. 1–2 by a set of changes of variables and
coordinate transformations. It describes the dynamics of excitable systems
which can be observed in various natural systems such as neuronal dynamics,
electrocardiology, chemical reactions, and climate dynamics. However, the
parameters in the FHN model are difficult to be computed directly in a real-
world experimentation and the estimation of these parameters has gained the
interest of a lot of researchers in physiological sciences. We shall see in
the following discussions that the specification of the model parameters is
crucial for the prediction of the system’s behavior. For instance, the system
can either converge to a stable fixed point or exhibit a limit cycle. Thus,
the knowledge of such parameters can be very useful for diagnostic as well as
prediction purposes, and the objective of the current study is to estimate the
model’s parameters from a few (possibly noisy) measurements of the system’s
state.

The parameter estimation framework for FHN model can be generally formulated
via standard techniques such as simulated annealing, genetic algorithms,
differential evolution, and Kalman filtering extensions. Besides, the known
model’s structure and characteristics can be utilized to customize an
algorithm to estimate the parameters of the respective model. For example, the
time-scale separation in the FHN model has been exploited to infer the model’s
parameters [6]. Che et al. [7] solved the parameter estimation problem by
deriving a second order differential equation for the membrane potential,
being the observed quantity. A least-squares based regression was then applied
and equipped by a wavelet denoising technique to reduce the effect of noise
contamination. Geng et al. [8] applied an expectation maximization based
algorithm to identify generic FHN model parameters and estimate the variance
of the interfering Gaussian noise. Jensen et al. [9] applied a Markov chain
Monte Carlo method to infer the parameters in a stochastic FHN model,
constructed by adding a noise term governed by a Brownian motion. Melnykova
[10] proposed a contrast estimator technique to infer the model’s parameters
in the asymptotic setting.

In the present study, we utilize a variational data assimilation technique,
namely the forward sensitivity method (FSM) [11, 12], to identify the correct
parameter values. The inherent sensitivity analysis reveals the relative
dependence of the cost functional, defined by the discrepancy between the
identified model’s predictions and the actual observations, onto the
respective parameters. We also investigate the effect of observation placement
instants on the shape of the cost functional and the corresponding
sensitivities. We finally highlight measurement collection guidelines that
potentially improve the parameter inference iterations.

##  2 Parameter Estimation Framework

The FHN model can be described as

$\dot{\mathbf{x}}=f(\mathbf{x},\boldsymbol{\alpha}),$ (3)

where $\mathbf{x}=[v(t),w(t)]^{T}$ denotes the system’s state,
$\boldsymbol{\alpha}=[a,b]^{T}$ is the model’s parameters, and $f$ represents
the continuous-time dynamics of the FHN model (i.e.,
$f(\mathbf{x};\alpha)=[v-\dfrac{1}{3}v^{3}-w+I,(v+a-bw)/\tau]^{T}$). Assuming
the the model $f$ is continuously differentiable in its arguments (i.e.,
$\mathbf{x}$ and $\boldsymbol{\alpha}$), its Jacobians with respect to the
state $\mathbf{x}$ and the parameter $\boldsymbol{\alpha}$ can be defined as
below

$\displaystyle Df_{\mathbf{x}}=\begin{bmatrix}1-v^{2}&-1\\\
1/\tau&-b/\tau\end{bmatrix},\qquad
Df_{\boldsymbol{\alpha}}=\begin{bmatrix}0&0\\\ 1/\tau&-w/\tau\end{bmatrix},$
(4)

where $Df_{\mathbf{x}}$ and $Df_{\boldsymbol{\alpha}}$ define the model’s
sensitivity with respect to the state $\mathbf{x}$ and the parameters
$\boldsymbol{\alpha}$, respectively.

###  2.1 Forward sensitivities

Using a suitable temporal integration scheme, the FHN can be rewritten in a
discrete-time form as follows,

$\mathbf{x}({k+1})=\mathbf{M}(\mathbf{x}(k),\boldsymbol{\alpha}),$ (5)

where $\mathbf{x}(k)=[v(t_{k}),w(t_{k})]^{T}\in\mathbb{R}^{2}$ defines the
system’s state at time $t_{k}$, and
$\mathbf{M}:\mathbb{R}^{2}\times\mathbb{R}^{2}\to\mathbb{R}^{2}$ represents
the one-step state transition map. Thus, the following discrete-time Jacobians
can be computed,

$\displaystyle\mathbf{D}\mathbf{M}_{\mathbf{x}}(k)=\bigg{[}\dfrac{\partial
M_{i}(\mathbf{x},\boldsymbol{\alpha})}{\partial
x_{j}}\bigg{]}_{\mathbf{x}=\mathbf{x}(k)},\qquad\mathbf{D}\mathbf{M}_{\boldsymbol{\alpha}}(k)=\bigg{[}\dfrac{\partial
M_{i}(\mathbf{x},\boldsymbol{\alpha})}{\partial\alpha_{j}}\bigg{]}_{\mathbf{x}=\mathbf{x}(k)}.$
(6)

Furthermore, we define the sensitivity of the model forecast at any time
$t_{k}$ with respect to the model’s parameters as follows,

$\mathbf{V}(k)=\bigg{[}\dfrac{\partial
x_{i}(k)}{\partial\alpha_{j}}\bigg{]}\in\mathbb{R}^{2\times 2}.$ (7)

Equation 5 can be used to evaluate the forward sensitivity matrices at
different times in a recursive way as

$\mathbf{V}(k+1)=\mathbf{D}\mathbf{M}_{\mathbf{x}}(k)\mathbf{V}(k)+\mathbf{D}\mathbf{M}_{\boldsymbol{\alpha}}(k),$
(8)

with $\mathbf{V}(0)=\mathbf{0}$ since the initial condition $\mathbf{x}(0)$ is
independent of the model’s parameters $\boldsymbol{\alpha}$.

###  2.2 Forecast error

Assuming $\mathbf{z}(k)\in\mathbb{R}^{m}$ to be the vector of measurements at
time $t_{k}$, given by

$\mathbf{z}(k)=\mathbf{h}(\bar{\mathbf{x}}(k))+\boldsymbol{\xi}(k),$ (9)

where $\mathbf{h}:\mathbb{R}^{2}\to\mathbb{R}^{m}$ defining the observational
operator that relates the model space to the observation space, and
$\bar{\mathbf{x}}$ defines the _true_ system’s state while $\boldsymbol{\xi}$
denotes the measurement noise. For simplicity, we suppose that we directly
measure the system’s state (i.e., $\mathbf{h}(\mathbf{x}(k))=\mathbf{x}(k)$).
We also assume that $\boldsymbol{\xi}$ is a white Gaussian noise with zero
mean and a covariance matrix $\mathbf{R}$ (i.e.,
$\boldsymbol{\xi}(k)={\mathcal{N}}(\mathbf{0},\mathbf{R}(k))$).

We define the difference between the model forecast and measurements as
$\mathbf{e}(k)=\mathbf{z}(k)-\mathbf{h}(\mathbf{x}(k))$, which is called the
innovation or forecast error (computed in the observation space). With the
assumption that the dynamical model is perfect (i.e., correctly encapsulates
all the relevant processes) and the initial condition $\mathbf{x}(0)$ is
known, then the deterministic part of the forecast error can be attributed to
the inaccuracy of the model’s parameters values, defined as
$\delta\boldsymbol{\alpha}=\bar{\boldsymbol{\alpha}}-\boldsymbol{\alpha}$,
where $\bar{\boldsymbol{\alpha}}$ denotes the true values of the parameters.
Thus, we can define a cost functional $J:\mathbb{R}^{2}\to\mathbb{R}$ as

$J(\boldsymbol{\alpha})=\sum_{k=1}^{N}\dfrac{1}{2}\|\mathbf{e}(k)\|^{2}_{\mathbf{R}^{-1}(k)}=\sum_{k=1}^{N}\dfrac{1}{2}\mathbf{e}(k)^{T}\mathbf{R}^{-1}(k)\mathbf{e}(k),$
(10)

where $N$ is the number of measurement instants. The minimization of the cost
function $J$ can be solved as a strong constrained problem with the standard
Lagrangian multiplier method, resulting in the adjoint framework.
Alternatively, we utilize the forward sensitivity matrices to evaluate an
optimal estimate for the parameters $\boldsymbol{\alpha}$. Let
$\delta\mathbf{x}(k)=\bar{\mathbf{x}}(k)-\mathbf{x}(k)$ be the difference
between the model’s forecast and the true state, with $\delta\mathbf{x}(0)=0$
since the initial conditions are perfectly known. With first order Taylor
expansions of $\mathbf{e}(k)$ and $\delta\mathbf{x}(k)$, the following
expressions can be defined,

$\displaystyle\mathbf{e}(k)=\mathbf{D}\mathbf{h}(k)\delta\mathbf{x}(k),\qquad\delta\mathbf{x}(k)=\mathbf{V}(k)\delta\boldsymbol{\alpha},$
(11)

where $\mathbf{D}\mathbf{h}$ is the Jacobian of the observational operator
$\mathbf{h}$. Therefore, the forecast error can be related to the correction
to the model’s parameters as
$\mathbf{e}(k)=\mathbf{D}\mathbf{h}(k)\mathbf{V}(k)\delta\boldsymbol{\alpha}$.
Since we assume that $\mathbf{h}(\mathbf{x}(k))=\mathbf{x}(k)$, we deduce that
$\mathbf{D}\mathbf{h}$ reduces to the identity matrix. The previous forecast
error formulation can be written for all $N$ time instants at which
observations become available, and the following linear equation is obtained,

$\mathbf{H}\delta\boldsymbol{\alpha}=\mathbf{e}_{F},$ (12)

where the matrix $\mathbf{H}\in\mathbb{R}^{Nm\times 2}$ and the vector
$\mathbf{e}_{F}\in\mathbb{R}^{Nm}$ are defined as follows,

$\mathbf{H}=\begin{bmatrix}\mathbf{D}\mathbf{h}(1)\mathbf{V}(1)\\\
\mathbf{D}\mathbf{h}(2)\mathbf{V}(2)\\\ \vdots\\\
\mathbf{D}\mathbf{h}(N)\mathbf{V}(N)\end{bmatrix},\qquad\mathbf{e}_{F}=\begin{bmatrix}\mathbf{e}(1)\\\
\mathbf{e}(2)\\\ \vdots\\\ \mathbf{e}(N)\end{bmatrix}.$ (13)

The inverse problem can be solved in a weighted least squares sense to find an
optimal correction vector $\delta\boldsymbol{\alpha}$, with $\mathbf{R}^{-1}$
as a weighting matrix, where $\mathbf{R}$ is a an $Nm\times Nm$ block-diagonal
matrix with $\mathbf{R}(k)$ being its $k$-th diagonal block. We assume that
$\mathbf{R}$ is a diagonal matrix defined as
$\mathbf{R}=\sigma^{2}\mathbf{I}_{Nm}$, where $\mathbf{I}_{Nm}$ is the
$Nm\times Nm$ identity matrix. Then, the solution to Eq. 12 can be written as

$\delta\boldsymbol{\alpha}=\left(\mathbf{H}^{T}\mathbf{R}^{-1}\mathbf{H}\right)^{-1}\mathbf{H}^{T}\mathbf{R}^{-1}\mathbf{e}_{F}.$
(14)

###  2.3 Placement of observations using forward sensitivity

In order to select the time instants at which measurement data are collected,
we relate the cost functional given in Eq. 10 to the forward sensitivity
matrix $\mathbf{V}(k)$. This is based on the method proposed by Lakshmivarahan
et al. [13] to control the shape of the cost functional and keep its gradient
away from zero to accelerate the convergence. By substituting
$\mathbf{e}(k)=\mathbf{D}\mathbf{h}(k)\mathbf{V}(k)\delta\boldsymbol{\alpha}$
into Eq. 10, we get the following,

$J(\boldsymbol{\alpha})=\sum_{k=1}^{N}\dfrac{1}{2}\delta\boldsymbol{\alpha}^{T}\bigg{(}\mathbf{E}(k)^{T}\mathbf{R}^{-1}(k)\mathbf{E}(k)\bigg{)}\delta\boldsymbol{\alpha}=\sum_{k=1}^{N}\dfrac{1}{2}\delta\boldsymbol{\alpha}^{T}\mathbf{G}(k)\delta\boldsymbol{\alpha},$
(15)

where $\mathbf{E}(k)=\mathbf{D}\mathbf{h}(k)\mathbf{V}(k)$ and
$\mathbf{G}(k)=\mathbf{E}(k)^{T}\mathbf{R}^{-1}(k)\mathbf{E}(k)$. We note that
$\mathbf{G}(k)$ is called the observability Gramian. The gradient of the cost
functional with respect to the parameter vector $\boldsymbol{\alpha}$ can be
written as below

$\displaystyle\nabla_{\boldsymbol{\alpha}}J(\boldsymbol{\alpha})$
$\displaystyle=\sum_{k=1}^{N}-\mathbf{G}(k)\delta\boldsymbol{\alpha},$ (16)

which relates the gradient of the cost functional and the parameterization
error/correction. From Eq. 16, a necessary condition for the minimization of
the cost functional is that $\mathbf{G}(k)$ is positive definite. For the case
considered here, $\mathbf{D}\mathbf{h}(k)=\mathbf{I}$ and
$\mathbf{R}^{-1}(k)=\dfrac{1}{\sigma^{2}}\mathbf{I}$. Thus,
$\mathbf{G}(k)=\dfrac{1}{\sigma^{2}}\mathbf{V}(k)^{T}\mathbf{V}(k)$, where
$\mathbf{V}(k)=\begin{bmatrix}V_{11}&V_{12}\\\ V_{21}&V_{22}\end{bmatrix}$.
Therefore, one way to guarantee that the gradient of the cost functional does
not hit zero and improve the convergence is to select the measurement instants
in such a way that the diagonal entries (i.e., $V^{2}_{11}+V^{2}_{21}$ and
$V^{2}_{12}+V^{2}_{22}$) are as large as possible.

##  3 Results and Discussions

We analyze the capability of the forward sensitivity approach to identify the
FHN model’s parameters. In particular, we study an arbitrary case where the
true parameters values are $a=0.15$ and $b=0.35$. Initial conditions of
$(v(0),w(0))=(0.0,1.0)$ are considered and the fourth order Runge-Kutta scheme
is applied for time integration with a time step of $\Delta t=0.1$, time scale
$\tau=10$, and a maximum time of $t_{m}=100$. We assume that the measurements
are collected every $200$ time steps, corrupted by an additive Gaussian noise
with a zero mean and a standard deviation of $\sigma=0.1$.

###  3.1 Fixed input

As a first investigation, we study the case with zero input (i.e., $I=0$).
This corresponds to a fixed point of $(v^{*},w^{*})=(-0.229,-0.225)$ with a
model Jacobian of $\begin{bmatrix}0.948&-1\\\ 0.1&-0.035\end{bmatrix}$. The
eigenvalues of this matrix are $\lambda_{1}=0.832,\lambda_{2}=0.080$, implying
unsteady equilibrium points. However, a Lyapunov function analysis reveals
that the solution of this system is bounded and exhibits an attractive limit
cycle [14, 15, 16, 17, 18]. In Figure 1, we plot the time evolution of the
membrane potential, $v$, and the recovery variable, $w$, for the true system
compared to the case with the inferred parameters values. Starting from a
prior guess of $a=0.2$ and $b=0.2$ to initiate the FSM iterations, a
parameterization of $a=0.159$ and $b=0.364$ is identified, very close to the
true values. Thus, we can see that the adopted FSM approach is adequately
capable of assimilating these noisy data to estimate the model’s parameters
for this case.

Figure 1: Results for $a=0.15$, $b=0.35$, and $I=0$, with measurements every
200 time steps and $\sigma=0.1$. Estimated parameters are $a=0.159$ and
$b=0.364$.

A second testing situation is to apply a constant input of $I=5$, with the
same parameters values as before. We find that this case corresponds to a
stable fixed point. In other words, the solution trajectory converges to the
equilibrium point (which is $(v^{*},w^{*})=(1.652,5.149)$) and resides there.
We apply the same procedure to estimate the model’s parameters starting with
an initial guess of $a=0.2$ and $b=0.2$. The plots in Figure 2 show that the
iterative algorithm fails to correctly approximate the parameters values and
produces a periodic solution, instead. To understand this, we compute the
fixed point and the eigenvalues of the corresponding model’s Jacobian. We find
that with $I=0$ (the previous case), both the true values $(a,b)=(0.15,0.35)$
and the initial guess $(a,b)=(0.2,0.2)$ induce a periodic limit cycle. On the
other hand, for $I=5$, the true parameter values correspond to a stable fixed
point, while the initial guess still yields a cyclic behavior. Therefore, the
estimation process should cross the bifurcation points in order to predict the
correct parameterization, which is a common problem in parameter estimation
frameworks.

Figure 2: Results for $a=0.15$, $b=0.35$, and $I=5$, with measurements every
200 time steps and $\sigma=0.1$. Estimated parameters are $a=-0.930$ and
$b=-0.021$, starting from an initial guess of $(0.2,0.2)$.

In order to mitigate this issue, prior information about the regime of the
solution trajectory can be utilized to make an intelligent guess. For
instance, an initial guess of $(a,b)=(0.5,0.5)$ with $I=5$ yields a stable
fixed point, and hence can be chosen as an alternative starting point. Results
are presented in Figure 3, where we can see that both the true and predicted
trajectories converge to the equilibrium state. However, the estimated
parameters values ($a=0.0344$ and $b=0.334$) are slightly far from the true
ones.

Figure 3: Results for $a=0.15$, $b=0.35$, and $I=5$, with measurements every
200 time steps and $\sigma=0.1$. Estimated parameters are $a=0.0344$ and
$b=0.334$, starting from an initial guess of $(0.5,0.5)$.

In order to explore the effect of the measurements on the forward
sensitivities, we plot the variation of $V^{2}_{ij}$ for $i,j\in\\{1,2\\}$
with time in Figure 4. We observe that the initial period has the least
influence on the forward sensitivities, while the measurements around and
after $t=50$ have the largest effects. Therefore, we redistribute our
measurement instants based on the approach described in Section 2.3. In
particular, we collect data at $t\in\\{30,35,50,60,70\\}$ and apply the FSM
framework to estimate the model’s parameter. Starting from an initial guess of
$(a,b)=(0.5,0.5)$, a parameterization of $(a,b)=(0.128,0.355)$ is estimated,
showing significant improvement with respect to the case with equispaced
measurement signals. Results are shown in Figure 5 for the true and predicted
trajectories. We can also notice that the optimized measurements are
concentrated towards the equilibrium state.

Figure 4: Forward sensitivities for $a=0.15$, $b=0.35$, and $I=5$. Selected
observations instants are denoted with orange circles. Figure 5: Results for
$a=0.15$, $b=0.35$, and $I=5$, with measurements instants selected based on
the forward sensitivity criteria. Estimated parameters are $a=0.128$ and
$b=0.355$, starting from an initial guess of $(0.5,0.5)$.

###  3.2 Varying input

Finally, we vary the input excitation as $I=5t/t_{m}$ (i.e., linearly
increasing from $0$ to $5$). This corresponds to a moving fixed point,
beginning with a cyclic trajectory and followed by a convergence to the stable
equilibria. Parameter estimation results for equidistant measurement instants
are depicted in Figure 6 beginning from an initial guess of $(a,b)=(0.2,0.2)$.
We find that the predicted trajectory sufficiently match the true one, but the
estimated parameters are not very accurate.

Figure 6: Results for $a=0.15$, $b=0.35$, and varying $I$, with measurements
every 200 time steps and $\sigma=0.1$. Estimated parameters are $a=0.023$ and
$b=0.330$, starting from an initial guess of $(0.2,0.2)$.

We then investigate the effects of observation times on the forward
sensitivities of the model predictions. We find a spike in the sensitivity of
$v$ predictions with respect to the $v$ measurements around $t=12.5$. We also
see a relatively large dependence on the $w$ measurements about $t=37.5$. On
the other hand, the $w$ predictions show an increasing sensitivity on either
$v$ or $w$ measurements at final times. Therefore, we reallocate our
observation times to capture these trends as demonstrated in Figure 7. Results
based on this enhanced parameter estimation methodology are described in
Figure 8, where the approximated parameters values ($a=0.128$, and $b=0.356$)
are closer to the true values.

Figure 7: Sensitivities for $a=0.15$, $b=0.35$, and varying $I$. Selected
observations instants are denoted with orange circles. Figure 8: Results for
$a=0.15$, $b=0.35$, and varying $I$, with measurements instants selected based
on the forward sensitivity criteria. Estimated parameters are $a=0.128$ and
$b=0.356$ starting from an initial guess of $(0.2,0.2)$.

##  4 Concluding Remarks

We put forth a forward sensitivity analysis for the FitzHugh-Nagumo (FHN)
system to infer the model’s parameterization from sparse observations. The
approach relies on the investigation of the forward sensitivity matrices that
encapsulates the temporal dependence of model’s predictions onto its
parameters. The presented methodology shows substantial success in
assimilating noisy observational data to identify the unknown parameters. We
find that the convergence of the predicted parameters to the true values
relatively depends on the first guess used to initialize the algorithm. In
particular, the initial guess has to yield equilibrium points with similar
stability characteristics to the true one. We study three test cases,
including zero input, constant non-zero current, and time-dependent
excitation. We also formulate measurement collection guidelines based on the
relation between the cost functional and the forward sensitivity components.
We demonstrate that this approach provides more accurate estimates of unknown
parameters than those resulting with arbitrary measurement placements.

## Acknowledgments

This material is based upon work supported by the U.S. Department of Energy,
Office of Science, Office of Advanced Scientific Computing Research under
Award Number DE-SC0019290. O.S. gratefully acknowledges the U.S. DOE Early
Career Research Program support.

Disclaimer: This report was prepared as an account of work sponsored by an
agency of the United States Government. Neither the United States Government
nor any agency thereof, nor any of their employees, makes any warranty,
express or implied, or assumes any legal liability or responsibility for the
accuracy, completeness, or usefulness of any information, apparatus, product,
or process disclosed, or represents that its use would not infringe privately
owned rights. Reference herein to any specific commercial product, process, or
service by trade name, trademark, manufacturer, or otherwise does not
necessarily constitute or imply its endorsement, recommendation, or favoring
by the United States Government or any agency thereof. The views and opinions
of authors expressed herein do not necessarily state or reflect those of the
United States Government or any agency thereof.

## References

  * [1] Eugene M Izhikevich.  Dynamical systems in neuroscience: The geometry of excitability and bursting.  MIT Press, 2007.
  * [2] Richard FitzHugh.  Mathematical models of threshold phenomena in the nerve membrane.  The bulletin of Mathematical Biophysics, 17(4):257–278, 1955\.
  * [3] Richard FitzHugh.  Impulses and physiological states in theoretical models of nerve membrane.  Biophysical Journal, 1(6):445, 1961.
  * [4] Jinichi Nagumo, Suguru Arimoto, and Shuji Yoshizawa.  An active pulse transmission line simulating nerve axon.  Proceedings of the IRE, 50(10):2061–2070, 1962.
  * [5] S Sehgal and AJ Foulkes.  Numerical analysis of subcritical Hopf bifurcations in the two-dimensional FitzHugh-Nagumo model.  Physical Review E, 102(1):012212, 2020.
  * [6] Rose T Faghih, Ketan Savla, Munther A Dahleh, and Emery N Brown.  The Fitzhugh-Nagumo model: Firing modes with time-varying parameters & parameter estimation.  In 2010 Annual International Conference of the IEEE Engineering in Medicine and Biology, pages 4116–4119. IEEE, 2010.
  * [7] Yanqiu Che, Li-Hui Geng, Chunxiao Han, Shigang Cui, and Jiang Wang.  Parameter estimation of the FitzHugh-Nagumo model using noisy measurements for membrane potential.  Chaos: An Interdisciplinary Journal of Nonlinear Science, 22(2):023139, 2012.
  * [8] Li-Hui Geng, Terefe Bayisa Ayele, Jin-Cang Liu, and Brett Ninness.  Expectation maximization based FitzHugh-Nagumo model identification under unknown gaussian measurement noise.  In 2020 Chinese Control And Decision Conference (CCDC), pages 2167–2172. IEEE, 2020.
  * [9] Anders Chr Jensen, Susanne Ditlevsen, Mathieu Kessler, and Omiros Papaspiliopoulos.  Markov chain Monte Carlo approach to parameter estimation in the FitzHugh-Nagumo model.  Physical Review E, 86(4):041114, 2012.
  * [10] Anna Melnykova.  Parametric inference for hypoelliptic ergodic diffusions with full observations.  Statistical Inference for Stochastic Processes, 23(3):595–635, 2020.
  * [11] S Lakshmivarahan and John M Lewis.  Forward sensitivity approach to dynamic data assimilation.  Advances in Meteorology, 2010:1–12, 2010.
  * [12] Sivaramakrishnan Lakshmivarahan, John M Lewis, and Rafal Jabrzemski.  Forecast error correction using dynamic data assimilation.  Springer, Switzerland, 2017.
  * [13] S Lakshmivarahan, John M Lewis, and Junjun Hu.  On controlling the shape of the cost functional in dynamic data assimilation: Guidelines for placement of observations and application to saltzman’s model of convection.  Journal of the Atmospheric Sciences, 77(8):2969–2989, 2020.
  * [14] Matthias Ringkvist.  On dynamical behaviour of FitzHugh-Nagumo systems.  PhD thesis, Department of Mathematics, Stockholm University, 2006.
  * [15] E Kaumann and U Staude.  Uniqueness and nonexistence of limit cycles for the FitzHugh equation.  In Equadiff 82, pages 313–321. Springer, 1983.
  * [16] KP Hadeler, U An Der Heiden, and K Schumacher.  Generation of the nervous impulse and periodic oscillations.  Biological Cybernetics, 23(4):211–218, 1976.
  * [17] SA Treskov and EP Volokitin.  On existence of periodic orbits for the FitzHugh nerve system.  Quarterly of Applied Mathematics, 54(4):601–607, 1996.
  * [18] Mattias Ringkvist and Yishao Zhou.  On the dynamical behaviour of FitzHugh–Nagumo systems: revisited.  Nonlinear Analysis: Theory, Methods & Applications, 71(7-8):2667–2687, 2009.


# Bimanual rope manipulation skill synthesis through context dependent
correction policy learning from human demonstration

T. Baturhan Akbulut1,∗, G. Tuba C. Girgin1,∗, Arash Mehrabi2, Minoru Asada3,4,
Emre Ugur1, Erhan Oztop2,3 ∗ Equal contribution, alphabetical order1 T.
Baturhan Akbulut, G. Tuba C. Girgin and Emre Ugur are with Deparment of
Computer Engineering, Bogazici University, Istanbul, Turkey. 2 Arash Mehrabi
and Erhan Oztop are with Department of Computer Engineering, Ozyegin
University, Istanbul, Turkey. 3 Erhan Oztop and Minoru Asada are with
Symbiotic Intelligent Systems Research Center, Institute for Open and
Transdisciplinary, Research Initiatives, Osaka University, Japan. 4 Minoru
Asada is also affiliated with International Professional University of
Technology in Osaka, Osaka, Japan.

###### Abstract

Learning from demonstration (LfD) provides a convenient means to equip robots
with dexterous skills when demonstration can be obtained in robot intrinsic
coordinates. However, the problem of compounding errors in long and complex
skills reduces its wide deployment. Since most such complex skills are
composed of smaller movements that are combined, considering the target skill
as a sequence of compact motor primitives seems reasonable. Here the problem
that needs to be tackled is to ensure that a motor primitive ends in a state
that allows the successful execution of the subsequent primitive. In this
study, we focus on this problem by proposing to learn an explicit correction
policy when the expected transition state between primitives is not achieved.
The correction policy is itself learned via behavior cloning by the use of a
state-of-the-art movement primitive learning architecture, Conditional Neural
Motor Primitives (CNMPs). The learned correction policy is then able to
produce diverse movement trajectories in a context dependent way. The
advantage of the proposed system over learning the complete task as a single
action is shown with a table-top setup in simulation, where an object has to
be pushed through a corridor in two steps. Then, the applicability of the
proposed method to bi-manual knotting in the real world is shown by equipping
an upper-body humanoid robot with the skill of making knots over a bar in 3D
space. The experiments show that the robot can perform successful knotting
even when the faced correction cases are not part of the human demonstration
set.

##  I INTRODUCTION

Learning from demonstration (LfD) [1] is an effective robot learning
framework, where a robot is taught a desired skill through human
demonstration. The demonstration may be kinaesthetic, visual, or through
teleoperation [2]. The latter places the human operator in the control loop.
Thus human motor learning capacity enables teaching more complex tasks to the
robot [3]. The learning can be based on policy cloning [4], i.e. learning a
state-to-action mapping, or simple motor-tape approach [5], i.e. playing back
the stored motor commands in an open-loop fashion. For complex and long tasks,
both methods suffer from the so called covariate shift problem [6]. When the
system enters an off-learned distribution state, it can not recover. Even
worse, it can diverge to states that are even further away. Assume a robot
executing a rather long trajectory in order to make a knot using its two
grippers. Even small mistakes the are made during the initial parts of the
action such as holding the rope from a slightly different position or making
slight errors in the movement trajectory that cause the gripper to miss the
rope, would result in failure of task execution. A remedy is to consider the
actions as the composition of simpler action primitives that are less prone to
covariate shift, allowing corrective actions to be inserted between primitives
to ensure the successful execution of the subsequent primitive. Corrective
action policies can be learned by the robot itself or again taught by a human
demonstrator. In this study, we explore the latter approach, first in
simulation, then show that it can be applied to real hardware by synthesizing
a complex rope manipulation skill on an upper-body humanoid robot (Fig 1).

Figure 1:  Bi-manual knotting with Torobo in the real world.

In summary, we propose a method that enables the robot to recover from failure
states that might be encountered during critical points in complex and long
action executions. To do this, the action is segmented based on a set of
identified critical points. During these critical points, the human operator
demonstrates the movements that needs to be taken in order to make the task
execution succeed, and the demonstration is learned by our system as a
function of robot-environment context. The corrective actions that might
require very different trajectories are encoded in a single corrective
movement primitive. Therefore, when encountered with a failure context, a
corrective action is automatically generated to ensure successful task
execution.

Overall, the main contributions of our paper are (1) reducing the human effort
for data generation for LfD by focusing the learning on correction policies,
(2) showing that the approach is general by realizing it on real hardware, and
(3) providing dexterous rope manipulation in 3D space superseding earlier rope
knotting skills with real robots.

##  II Related Work

There has been a significant amount of research on robotic manipulation of
deformable objects. Nair et al. [7] combined self-supervised learning with
imitation for rope manipulation. First, they trained an inverse dynamics model
that predicts the action given the images of the current and next states.
Then, they used human demonstration as a high-level plan of what to do and the
inverse dynamics model to determine how to do it. The training data, which
consisted of 60K pairs of images, were used for training the inverse dynamics
model. The inverse dynamics model of [7] works well in predicting the actions
that make small deformations, however it lacks the ability to construct a
trajectory to reach to goal states. Wang et al. [8] extended the ideas of Nair
et al. [7] by replacing human demonstration with a Casual InfoGAN (CIGAN)
model. The CIGAN can generate a trajectory to reach a given goal state from
the initial state. However, high success rates could not be reached in
generating plausible plans, as explained in [8].

Yan et al. [9] developed a visual model-based learning framework that jointly
learns the underlying visual latent representations and the dynamics models
for deformable objects. The latent dynamics model was learned in a simulation
environment, then transferred to a real PR2 robot. In another work, Yan et al.
[10] proposed a model that learns the dynamics of linear deformable objects in
state space. However, bi-manual rope manipulation was not addressed in either
of the works, and they did not consider local corrective movements targeted at
critical points as in our work.

Chang et al. [6] proposed the network MILO where offline RL achieves learning
without needing online demonstrations. This method suffers from out-of-
distribution configurations that a robot can encounter in real-life
experiments. Reichlin et al. [11] suggested a network to combine it with the
offline RL setup to overcome the distribution-shift problem. The proposed
recovery network encodes observations and creates a latent representation. The
agent corrects its motion by using the encoded information and approximating
the training dataset’s density in latent space. However, the recovery network
affects the predicted motion in each step, while our proposed method only
needs correction in critical points, simplifying its logic.

Multiple other approaches have been tried for deformable object manipulation.
Notably [12] realized a model-free visual reinforcement learning algorithm for
fabric and rope manipulation both in simulated and real robot environments,
and a max-margin learning algorithm for rope untangling demonstrated in the
physical environment. Sundaresan et al. [13] used an interpretable deep visual
representations for rope. They learned point-pair correspondences between
initial and goal rope configurations in simulation from synthetic depth data,
enabling them to plan actions to reach the goal configuration, and applied
their method also on a real robot. In another study, Seita et al. [14] used
deep imitation learning for fabric smoothing. They trained policies in
simulation using dataset aggregation (DAgger) with an algorithmic supervisor,
given images of a fabric sample as input. The difference between our algorithm
and DAgger is that our algorithm only asks the expert for the correct action
when it reaches a critical state, while DAgger needs continuous access to the
supervisor’s policy.

Unlike previous works, our algorithm learns to recover from failure states. We
tackle the problem of learning an explicit correction policy using expert
demonstration when the robot encounters a state where executing subsequent
primitive is not possible.

Figure 2: Illustration of the proposed concept with a single critical point.
Demonstrated trajectory (dashed gray curve) passes through a critical region
C, which depends on the task. For the robot execution, a correction movement
is inserted to counteract the possible case of execution missing C, thereby
splitting the execution into three motor primitives (MP): pre-MP, cor-MP, and
post-MP. post-MP must be learned context-dependent to recover a wide range of
failure cases. This is achieved by obtaining multiple correction
demonstrations through teleoperation.

##  III Method

Figure 3: Conditional Neural Movement Primitive architecture (CNMP) with a
context input is shown. In the illustration, the context is provided by an
autoencoder. See Section III for more details.

Complex tasks, e.g., rope manipulation, are often linked through critical
points significant for successfully executing the task, where critical points
are considered to be certain parts of the execution trajectory that may lead
to failure. Furthermore, long and complex tasks are fragile to the
accumulation of errors. Therefore, we propose to segment the trajectory by
using the critical points as the segment boundaries and insert corrective
motor primitives to ensure the subsequent segment’s preconditions are met.
Fig. 2 illustrates a trajectory with two such segments, encoded with two
movement primitives, namely Pre-Movement Primitive (pre-MP) and Post-Movement
Primitive (post-MP), and a single critical point where a corrective movement
primitive (cor-MP) enables the robot to recover from different failure
situations. In this work, we assume the critical points are known and focus on
cor-MP learning through human demonstration. Since the correction would depend
on the robot and environment configuration, i.e. the context reached after the
execution of the preceding segment, we learn a context-dependent correction
policy. In plain terms, a context is a vector extracted from sensors available
to the robot that captures the task-relevant environment-robot relations. By
imposing such a task decomposition, rather than learning a monolithic task
policy from demonstration, we aim to improve performance and robustness, while
achieving a general correction policy for the critical point under
consideration. To achieve a similar performance from a monolithic behavior
cloning system, one would need to ask for multiple demonstrations for the full
action; whereas in our proposal, the corrections for the critical points have
to be demonstrated multiple times while the main task execution can be
demonstrated a few times (typically once suffices).

###  III-A General Framework

Our model is built on top of Conditional Neural Movement Primitives (CNMPs)
architecture [15]. CNMPs is a deep neural architecture capable of representing
sensorimotor data conditioned on the given set of observation points as a
latent vector, which then can be used to unfold the trajectory based on
queried time points. Thus, when coupled with a feedback controller, it can act
as a drive for the movements of a robot. In our work, for the critical segment
learning, we feed CNMPs with the context information so that the correction is
conditioned on the context. A CNMPs is composed of two networks: an encoder
and a decoder, as shown in Fig. 3. The encoder part learns to generate a
latent representation (i.e. latent vector) from the samples taken from an
input trajectory. By taking the mean of the representation for each sample, a
general representation is obtained that is sufficient to reconstruct the
trajectories learned. Since training and querying are done with a context
input appended to the latent representation (Fig. 3), the generated
trajectories become context-dependent. Thus, for the corrective motor
primitives (cor-MP), CNMPs learn to generate different motor patterns based on
the context. Priming the CNMPs with a context not only enforces the selection
of the right correction policy but also allows generalization to unseen
contexts, thereby bringing robustness to overall task execution.

In this study, we test our method on two tasks realized in simulated and real
robots. Both tasks involve a single critical region. Therefore, we can
unambiguously refer to pre-critical motor primitive (pre-MP), corrective motor
primitive (cor-MP), and post-critical motor primitive(post-MP) in our
descriptions (see Fig 2). Further details are explained comprehensively in
Section IV.

For the simulation task, human demonstration is emulated through a
deterministic algorithm. On the other hand, for the hardware task, human
demonstration is obtained through a teleoperation system that is developed by
using two haptic devices (GeoMagic Inc.).

###  III-B Learning the Context

The context vector provides information about the current robot state and
environment relation. As such, it can be a manually designed feature
extractor, or it can be learned in an end-to-end fashion through an auto-
encoder (AE). To test the latter approach, for the simulation experiment we
designed an AE with 2D convolutional layers of kernel sizes 3, strides [2, 1,
2, 1, 2], and filter sizes [32, 32, 64, 64, 64] in the encoder block. The
output of the encoder is flattened and fed into a fully connected layer to get
the latent representation of the image. The decoder mimics the encoders
operations in reverse and used only in training the AE. In the simulation
experiments, the encoder receives an image and the resulting latent vector is
used as the context. The AE is trained on a data set consisting of 400 images.
360 of them are used for training and the remaining 40 for testing. The images
are taken after the first push and correction actions (described in Section
IV-A). The AE is trained for 70 epochs using Adam optimizer [16] with a
learning rate of 3e-4. The size of the latent space representation of images
was set to 8 for the experiments reported in the current paper.

Figure 4:  Simulation environment prepared in Gazebo. Figure 5:  A typical
action execution to accomplish the goal.

##  IV Experiments and Results

###  IV-A Push object through corridor with AE-based context

We use the simulated upper body humanoid robot Torobo (Tokyo Robotics Inc.) in
Gazebo simulation environment (Fig. 4) in order to teach the robot to push a
target object through a corridor in two steps: First the robot needs to align
the object with the corridor, and then push it trough, which is a simple task
the robot can perform. However, to test our method, we, perturb the position
of the corridor stochastically in a 6 cm range, thereby deliberately creating
a critical region for the task execution that would require a corrective step.
Following the logic of the proposed model, the robot asks demonstrations for
the those critical points that are encountered after the first push. As stated
in Section III, the human demonstration is emulated by an algorithmic logic
that guides the robot to do the ’right’ corrective action. While the
correction is being carried out, the data collected as the corrective action
demonstration data, and used for training cor-MP neural network.

The characteristics of the corrective motor primitives (cor-MP) differ
according to the position of the corridor. When the corridor moves away from
the robot, the robot needs to push the object to align it with the corridor.
On the other hand, if the corridor comes closer, the robot needs to pull the
object so as to align it with the corridor. A typical task execution with a
pull type correction is shown in Fig. 5.

Two hundred trajectories are collected together with environment information
in total to create a data set. Each data unit consists of the Cartesian
position of the gripper, the corridor location, and the target object. In
addition, for each trajectory, two images are captured from the camera of the
Torobo Robot. The first image is captured after the execution of the pre-
critical movement primitive (pre-MP). The second image is captured after the
execution of cor-MP. These image are used to obtain autoencoder based context
for the experiments that involve learned contexts (Section IV.B.2). Typical
images captured from the camera of the robot head camera is shown in Fig. 6.

TABLE I: Simulation Training Parameters Network Type |  Learning Rates |  Encoder Layer Count |  Decoder Layer Count |  Layer Size
---|---|---|---|---
All CNMPs |  0.0001 |  3 |  4 |  128

As mentioned in Section III, the task execution is composed of the executions
of three sequential motor primitives. Thus our system is required to produce
three trajectories, namely pre-MP network, cor-MP network, post-MP network to
successfully complete the task. The trajectories are encoded in task space in
this experiment. Particularly, the context is used by the post-MP network to
make the final push since after cor-MP network correction, the location of the
object would be, in general, changed. Two simulation experiments are
conducted, one with ’perfect context’, and one with the ’AE-based’ context.
For the perfect context, the corridor location and the object location after
the correction are used. To be concrete, cor-MP network uses the position of
the corridor on the posterior-anterior axis as the context, and post-MP
network uses the position of the target object on the same axis. The context
values are normalized before being used in cor-MP network and post-MP network.

(a)

(b)

Figure 6:  Camera images of the scene, showing the configurations before
(left) and after (right) the correction action.

For the AE-based context based experiments, the captured images from the
camera are fed into an autoencoder (AE) to extract the information related to
the target object and the corridor positions. Then, the encoded images (i.e
latent vectors) are given as contexts to the cor-MP network and post-MP
network, respectively. Table I shows the learning parameters used for
training. To reiterate, the CNMP networks are trained by using (normalized)
time and x, y, and z positions of the gripper collected during demonstrations
for both perfect and AE-based context experiments. The CNMP weights are
trained by using Adam optimizer[16]. During the training, 100 of 200
trajectories are used for training our model for each context type. The
remaining ones are used for validation. Each model received 10M weight
updates.

###  IV-B Results of Push through corridor task Experiment (Simulation)
Experiments

TABLE II: The action decision and task success of our proposed architecture with perfect context (PC) and autoencoder features (AE) are provided. For comparison, the results obtained from monolithic networks (Mono) is also provided.  Expected Correction | Action Choice | Decision Match | Task Success
---|---|---|---
PC | AE | Mono | PC | AE | Mono | PC | AE | Mono
Push | 9 Push | 8 Push, 1 No Correction | 9 Push | $100\%$ | $89\%$ | $100\%$ | $100\%$ | $100\%$ | $44\%$
No Correction | 6 No Correction | 6 No Correction | 6 No Correction | $100\%$ | $100\%$ | $100\%$ | $83\%$ | $100\%$ | $67\%$
Pull | 10 Pull | 9 Pull, 1 No Correction | 10 Pull | $100\%$ | $90\%$ | $100\%$ | $100\%$ | $100\%$ | $90\%$
TABLE III: Prediction Errors (cm) of MP trajectories.  Network Type | cor-MP | post-MP
---|---|---
mean | max | mean | max
PC | 0.80 | 1.45 | 1.69 | 3.19
AE | 0.80 | 0.99 | 2.64 | 4.83
B | 1.30 | 1.92 | 5.75 | 7.24

The trained pre-MP, corr-MP and post-MP networks are integrated into the
simulation environment to assess the performance of the proposed system in
driving the Torobo Robot for successfully pushing the target object through
the corridor. For execution, the predicted pre-MP Cartesian trajectory is
mapped to a joint trajectory through inverse kinematics which is fed to the
robot system as desired joint angles at a nominal rate. The robot controller
then tracks the obtained joint trajectory and ends up in the critical region.
The robot observes the scene, obtains the latent vector from the trained auto-
encoder (in the case of AE based context) or logs the position of the corridor
(in the perfect-context case), and uses it as the context for the cor-MP
network. By taking the last gripper position as the initial observation for
cor-MP network, cor-MP trajectory is predicted. Since the cor-MP network has
learned corrections from demonstrations, the generated trajectories are
expected to correspond to push, pull, or no action (as these have been
strategies used in the demonstrations). After executing cor-MP, the post-MP
trajectory is predicted according to the observed context and the position of
the end effector. After the execution, the trial is counted as a success one
if the target center of the object passes the corridor’s midpoint without
touching the corridor’s entrance. To asses the performance of the system 25
trials (test data set) with unseen contexts are executed with the predicted
pre-MP, cor-MP, and post-MP for both perfect context and AE-based context
experiments.

####  IV-B1 Results of CNMPs with Perfect Context (PC-CNMPs)

In order to successfully complete the task, the robot was required to make 9
push and 10 pull correction actions in the 25 cases tested. The remaining 6
cases required no correction (as the object was aligned with the corridor to a
tolerable level). The PC-CNMPs results show that the network was able to
predict an appropriate correction trajectory when faced with previously unseen
contexts. This shows that the cor-MP network can learn to generate correct
decisions and correct movement amounts. As a result of these corrective
actions, the robot was able to move the object to a position to align it with
the corridor. The only task failure was detected during a push type of
correction generated by the cor-MP network, and the failure was due the object
slightly contacting the corridor at the entrance.

The mean and maximum error values for cor-MP and post-MP are shown in Table
III. The error for each motion is calculated by taking the average of the
Euclidean distances between the predicted action points and the ground truths.
The mean error is 0.80 cm for cor-MP and 1.69 cm for post-MP. post-MP error is
higher than the previous one because the desired motion for post-MP is more
complex than the previous one. These errors cause small noise in the end
position of the target object after each motion. Because these errors are
small enough, 24 trials out of 25 yielded successful task executions.
Predictions, robot executions, and the ground-truths for two different action
type for cor-MP are compared in Fig. 7. As seen in the Figure, the predictions
overlap the ground truths. Therefore, the simulated robot was able to execute
the predicted motions accordingly.

####  IV-B2 Results of CNMPs with AE (AE-CNMPs)

As shown in Table II, when AE features were used as context by our system, in
24 trials out of 25, the robot was observed to accomplish the task by pushing
the target object in front of the narrow corridor, correcting its position,
and pushing it through the narrow corridor to the goal position. When the
execution of each predicted motion is examined, the cor-MP network was mostly
able to choose the right action type given the context. When we examine the
predictions that perform differently from the test data set used to assess the
performance of the predictions it is observed that the predicted trajectories
by the cor-MP network slightly shifted in comparison to demonstrated
trajectories. Since the captured image at boundary cases between push and no-
correction as well as pull and no-correction are similar, the AE
representation for those cases are slightly smeared. Therefore, the cor-MP
network may choose to generate different action types around the decision
boundaries. Still, because the correction amount that is required around the
decision boundaries is small compared to the width of the corridor, the
different decisions at these boundaries have not affected the success of task
execution. As a result, the success of the system with AE context has been
found to be 100%.

Figure 7:  Comparison of ground-truth, prediction, and robot execution for
different cor-MPs.

The mean cor-MP error is similar to that of perfect context version of the
system. The mean error resulted from post-MP network prediction is 2.64 cm,
which is higher than the one found for the perfect context experiments. This
might be due to the fact that the AE based context version of the system made
different choices when doing the correction, leaving the target object in
different location experienced during the demonstrations, which is acted upon
by post-MP.

Figure 8:  Knot over bar execution illustrated with Torobo robot.

####  IV-B3 Comparison with Monolithic Baseline CNMP

To assess the benefits of using segment action representation with explicit
corrective motor primitives for the critical regions, we implement a
monolithic baseline CNMP network. Given a fixed context, the full trajectory
can be learned with the Vanilla CNMPs. However, since the position of the
corridor (i.e. context) changes trial-to-trial, without knowing any
information about the context, the learning system cannot decide which
trajectory to generate. Therefore, to be fair to the baseline case, we provide
the corridor location to the Vanilla CNMPs and ask it to learn the
demonstrated full trajectories (which include the proper corrections). The
number of parameters of Vanilla CNMPs is adjusted to be the sum of number of
parameters of the pre-MP, cor-MP and post-MP networks The comparison of the
Vanilla CNMPs network generated actions to those of PC-CNMPs and AE-CNMPs are
shown in Table II with the Mono label. As shown, the Vanilla CNMPs network can
generate the right correction action; however, when executed, only 17 trials
out of 25 succeed in completing the task. This is due to the errors stemming
from learning a longer trajectory, albeit the Vanilla CNMPs network uses an
equal number of parameters as of the sum of the PC-CNMPs and AE-CNMPs
networks. The mean and max error values for the Vanilla CNMPs network are also
observed to be higher than the PC-CNMPs and AE-CNMPs networks. The results
show that the proposed method outperforms the baseline with the total number
of learnable parameters equalized.

###  IV-C Knotting with the real robot

For the hardware implementation we chose bi-manual rope knotting over a bar,
as it demonstrates the applicability of the proposed method for a range of
complex tasks. With the implemented teleoperation setup (Section III), the
users can control the pose of the individual grippers attached to the arms of
the robot Torobo robot. With this, we obtain basic knotting demonstration on
the robot by an expert demonstrator. Learning from demonstration with the
Vanilla CNMPs yields knotting behavior; however small variations such as how
the gripper holds the rope ends, the length of the rope between the grippers,
and invisible twists in the rope cause the robot to miss the critical ‘pass
through the loop and pull the free end’ (see Fig. 8) part of the knotting
action. So, analogous to the simulation experiments, at the point just before
the gripper would enclose (see Fig. 8-C), we collected multiple corrective
demonstrations from the user to synthesize a corrective motor primitive to be
inserted at this critical point with the aim of increasing the success of
knotting. The actions before and after this point i.e., (pre-MP, post-MP) can
be learned easily by individual CNMPs. For the corrective part, we have
collected 28 corrective action demonstrations which involve posterior-to-
anterior, anterior-to-posterior, lateral corrections and their combinations. A
few demonstrations without any corrections were also included. All the
demonstrations were learned with a single corrective CNMP (cor-MP network) as
in the simulations giving us a corrective motor primitive (cor-MP). The
correction was given as relative Cartesian positions with respect to the
gripper position. Thus the correction could be applied in any place in the
workspace of the robot. For the Cartesian to joint space transformation we
used inverse kinematics with fixed gripper orientation (as attained after the
execution of pre-MP). As the context information, the cor-MP network received
the relative position of the rope end with respect to the gripper mid-point.
The rope end position is extracted with a straightforward color-segmentation
and clustering pipeline on the RGB-D image captured by the camera (Intel
SR300) mounted on the head of Torobo.

Figure 9:  Real robot knotting performance. Knot correction demonstrations
(gray) with predictions generated (green) and the corresponding contexts (left
panel).

The cor-MP network could learn the demonstrations well, thus near exact
context (i.e. rope-gripper configuration) at the critical point successfully
generated corrective movements. More critical test was to see the performance
on unseen contexts. For that, we asked the robot to execute pre-MP and perform
the correction by using the trajectory generated by the cor-MP network, which
was followed by the execution of the post-MP, that included the enclose
command of the gripper. Thus if the correction was successful, the knot
operation would be completed successfully. In those trials, if the reached
critical points did not require a correction, we moved the rope outside the
fingers to see whether the required correction can be made. The experimental
results are provided in Fig. 9. The gripper-robot configurations used in
training are shown with gray, and the test configurations with successful
correction and task execution are shown with green. In addition some
indicative failures are also shown (yellow and red). Overall, we can conclude
that the correction skill was acquired by the robot and generalization to
inexperienced contexts could be handled. Most of the failures were to due to
kinematic constraints (i.e. out of reach) and others are near misses (the rope
is contacted by one of the gripper fingers but the rope was not straddled).

##  V Discussion and Future Work

In this paper, we have proposed a novel LfD approach which can be applied to
deformable object manipulation. By segmenting a given task into motor
primitives and inserting explicit correction policies between the segments, we
show that the potential failures that would be faced with a monolithic
approach are eliminated, and accuracy and generalization of the task execution
can be increased, while at the same human effort for demonstration is reduced.

In the current paper, we have focused on learning correction policies. As such
the success of the proposed method relies on the knowledge on which regions
are likely to need correction, i.e. critical for task success. Finding these
regions automatically without excessive number of full task demonstrations
would make the proposed method much more potent and applicable to tasks where
domain knowledge is not readily available.

## References

  * [1] S. Schaal, “Learning from demonstration,” _Advances in neural information processing systems_ , vol. 9, 1996.
  * [2] N. Amirshirzad, A. Kumru, and E. Oztop, “Human adaptation to human–robot shared control,” _IEEE Transactions on Human-Machine Systems_ , vol. 49, no. 2, pp. 126–136, 2019.
  * [3] J. Babic, J. G. Hale, and E. Oztop, “Human sensorimotor learning for humanoid robot skill synthesis,” _Adaptive Behavior_ , vol. 19, no. 4, pp. 250–263, 2011. [Online]. Available: ¡Go to ISI¿://000294157400002
  * [4] D. A. Pomerleau, “Alvinn: An autonomous land vehicle in a neural network,” in _Advances in Neural Information Processing Systems_ , D. Touretzky, Ed., vol. 1. Morgan-Kaufmann, 1988.
  * [5] C. Atkeson, J. Hale, F. Pollick, M. Riley, S. Kotosaka, S. Schaul, T. Shibata, G. Tevatia, A. Ude, S. Vijayakumar, E. Kawato, and M. Kawato, “Using humanoid robots to study human behavior,” _IEEE Intelligent Systems and their Applications_ , vol. 15, no. 4, pp. 46–56, 2000.
  * [6] J. Chang, M. Uehara, D. Sreenivas, R. Kidambi, and W. Sun, “Mitigating covariate shift in imitation learning via offline data with partial coverage,” in _Advances in Neural Information Processing Systems_ , M. Ranzato, A. Beygelzimer, Y. Dauphin, P. Liang, and J. W. Vaughan, Eds., vol. 34. Curran Associates, Inc., 2021, pp. 965–979.
  * [7] A. Nair, D. Chen, P. Agrawal, P. Isola, P. Abbeel, J. Malik, and S. Levine, “Combining self-supervised learning and imitation for vision-based rope manipulation,” in _2017 IEEE International Conference on Robotics and Automation (ICRA)_ , 2017, pp. 2146–2153.
  * [8] A. Wang, T. Kurutach, K. Liu, P. Abbeel, and A. Tamar, “Learning robotic manipulation through visual planning and acting,” 2019. [Online]. Available: https://arxiv.org/abs/1905.04411
  * [9] W. Yan, A. Vangipuram, P. Abbeel, and L. Pinto, “Learning predictive representations for deformable objects using contrastive estimation,” 2020. [Online]. Available: https://arxiv.org/abs/2003.05436
  * [10] M. Yan, Y. Zhu, N. Jin, and J. Bohg, “Self-supervised learning of state estimation for manipulating deformable linear objects,” _IEEE Robotics and Automation Letters_ , vol. 5, no. 2, pp. 2372–2379, 2020.
  * [11] A. Reichlin, G. L. Marchetti, H. Yin, A. Ghadirzadeh, and D. Kragic, “Back to the manifold: Recovering from out-of-distribution states,” 2022. [Online]. Available: https://arxiv.org/abs/2207.08673
  * [12] W. H. Lui and A. Saxena, “Tangled: Learning to untangle ropes with rgb-d perception,” in _2013 IEEE/RSJ International Conference on Intelligent Robots and Systems_ , 2013, pp. 837–844.
  * [13] P. Sundaresan, J. Grannen, B. Thananjeyan, A. Balakrishna, M. Laskey, K. Stone, J. E. Gonzalez, and K. Goldberg, “Learning rope manipulation policies using dense object descriptors trained on synthetic depth data,” 2020. [Online]. Available: https://arxiv.org/abs/2003.01835
  * [14] D. Seita, A. Ganapathi, R. Hoque, M. Hwang, E. Cen, A. K. Tanwani, A. Balakrishna, B. Thananjeyan, J. Ichnowski, N. Jamali, K. Yamane, S. Iba, J. Canny, and K. Goldberg, “Deep Imitation Learning of Sequential Fabric Smoothing From an Algorithmic Supervisor,” in _IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_ , 2020.
  * [15] M. Y. Seker, M. Imre, J. Piater, and E. Ugur, “Conditional neural movement primitives,” in _Proceedings of Robotics: Science and Systems_ , FreiburgimBreisgau, Germany, June 2019.
  * [16] D. P. Kingma and J. Ba, “Adam: A method for stochastic optimization,” _arXiv preprint arXiv:1412.6980_ , 2014.