# The complete inverse Kohn-Sham problem: from the density to the energy

A. Liardi  Dipartimento di Fisica “Aldo Pontremoli”, Università degli Studi di
Milano, 20133 Milano, Italy  F. Marino<EMAIL_ADDRESS>Dipartimento
di Fisica “Aldo Pontremoli”, Università degli Studi di Milano, 20133 Milano,
Italy  INFN, Sezione di Milano, 20133 Milano, Italy  G. Colò  Dipartimento di
Fisica “Aldo Pontremoli”, Università degli Studi di Milano, 20133 Milano,
Italy  INFN, Sezione di Milano, 20133 Milano, Italy  X. Roca-Maza
Dipartimento di Fisica “Aldo Pontremoli”, Università degli Studi di Milano,
20133 Milano, Italy  INFN, Sezione di Milano, 20133 Milano, Italy  E. Vigezzi
INFN, Sezione di Milano, 20133 Milano, Italy

###### Abstract

A complete solution to the inverse problem of Kohn-Sham (KS) density
functional theory is proposed. Our method consists of two steps. First, the
effective KS potential is determined from the ground state density of a given
system. Then, the knowledge of the potentials along a path in the space of
densities is exploited in a line integration formula to determine numerically
the KS energy of that system. A possible choice for the density path is
proposed. A benchmark in the case of a simplified yet realistic nuclear system
is shown to be successful, so that the method seems promising for future
applications.

###### pacs:

21.60.Jz, 21.60.-n

##  I Introduction

The Inverse Kohn-Sham (IKS) problem Jensen and Wasserman (2018); Shi and
Wasserman (2021) aims at reverse engineering the two steps that characterise
Density Functional Theory in the Kohn-Sham direct scheme (KS-DFT) Parr and
Yang (1994); Colò (2020). In the direct DFT problem, one usually assumes an
Energy Density Functional, or EDF, $E[\rho]$. Minimisation of this functional,
in the form $\delta E[\rho]=0$, provides a formulation of the ground state
(g.s.) problem, while the Runge-Gross theorem sets the formal basis for the
study of excited states within a time dependent DFT framework Runge and Gross
(1984).

Specifically, in the KS-DFT scheme a single-particle (s.p.) representation is
employed. Indeed, an auxiliary system of independent particles with the same
g.s. density as the true interacting system is introduced. The KS density,
then, is a function of the s.p. orbitals $\phi_{\alpha}(\mathbf{r})$ and reads

$\rho(\mathbf{r})=\sum_{\alpha}\phi^{*}_{\alpha}(\mathbf{r})\phi_{\alpha}(\mathbf{r}).$
(1)

The EDF is customarily written as

$E[\rho]=T+F+V_{\rm ext},$ (2)

where $T$ is the kinetic energy of a system of independent particles, $F$ is
the internal potential energy of the auxiliary KS system, and $V_{\rm ext}$ is
the contribution to the energy due to an external potential
$v_{ext}(\mathbf{r})$. A potential $v[\rho]$ is associated to $F$ by means of

$v[\rho]\equiv\frac{\delta F}{\delta\rho}.$ (3)

Then, the variational equation $\delta E[\rho]=0$, with the constraint that
the orbitals form an orthonormal set, gives rise to the so-called KS
equations, that are a set of one-particle Schrödinger-like equations. Hence,
the first step of the direct KS problem consists in calculating the effective
potential $v[\rho]$ from an assumed form for $F$. The second step consists in
solving the KS equations in order to determine the s.p. orbitals and the
corresponding density (1). In short,

$E[\rho]\ \ \ \rightarrow\ \ \ v[\rho]\ \ \ \rightarrow\ \ \ \rho$ (4)

is a view of the direct KS problem.

The IKS problem, that is

$\rho\ \ \ \rightarrow\ \ \ v[\rho]\ \ \ \rightarrow\ \ \ E[\rho],$ (5)

is far less simple. Inverse problems are in general difficult to tackle from a
mathematical and computational viewpoint. First, we should stress that there
are many more attempts and results concerning the first branch (density-to-
potential, hereafter D2P) than for the second one (potential-to-energy,
hereafter denoted by P2E). As for the D2P step, the first techniques to solve
the IKS problem have been proposed in the context of electronic DFT in Refs.
Görling (1992); Wang and Parr (1993); van Leeuwen and Baerends (1994); Zhao
_et al._ (1994). This topic has recently encountered a renewed interest both
in quantum chemistry Kanungo  _et al._ (2019); Naito  _et al._ (2019); Nam
_et al._ (2020, 2021); Pérez-Jordá (2020) and in nuclear physics Accorto  _et
al._ (2020, 2021) (our group had undertaken preliminary steps before Ref.
Accorto  _et al._ (2020), cf. e.g. Macocco (2015)). Several inversion
techniques that are relevant in this context are reviewed in the useful Refs.
Jensen and Wasserman (2018); Shi and Wasserman (2021); Kumar  _et al._ (2019),
and inversion codes either have been developed ex novo Nam  _et al._ (2021) or
have been integrated into existing libraries, such as Octopus Andrade  _et
al._ (2015). Complementary to these practical implementations, we should also
mention a recent attempt to delve more into the mathematical aspects of the
IKS problem Garrigue (2021).

However, in spite of the indubitable technical progress Kanungo  _et al._
(2019), the first branch of the IKS problem has found per se only a limited
number of applications. So far, indeed, the reverse-engineered KS potentials
have been mostly employed to benchmark existing models Umrigar and Gonze
(1994); Nam  _et al._ (2020). In our view, the topic of utmost interest is the
determination of the EDF itself. To this purpose, the P2E step is necessarily
involved. Among the few attempts in this direction, we mention that machine
learning approaches to DFT may benefit from the knowledge of the exact KS
potential Kanungo  _et al._ (2019). For example, in Ref. Schmidt  _et al._
(2019) it has been shown that training a neural network EDF not only on
energies, but also on potentials, greatly improves its performances, at least
in a simplified one-dimensional framework. A different line of research
combines density functional perturbation theory (DFPT), where an ansatz for
the EDF functional form must be assumed, and IKS to improve the initial model
towards the exact EDF Naito  _et al._ (2019); Accorto  _et al._ (2021).
Illustrative calculations have been performed in the case of covariant nuclear
DFT Accorto  _et al._ (2021).

A more direct approach to find the EDF from a solution of the whole IKS
problem, as sketched in (5), is our main interest in this work. The second
branch, P2E, has been the subject of some theoretical investigation, although
not of a full-fledged solution. A line integration formula has been proposed
by R. Van Leeuwen and E.J. Baerends in van Leeuwen and Baerends (1995) as a
mean to perform the functional integration and obtain the total energy $E$ of
a given system with density $\rho$ from $v$. Such formula requires the
knowledge of the system densities not only in the ground state but also along
a given path, together with the corresponding potentials. The relation between
the potential $v$ and the EDF has been thoroughly investigated by A.P. Gaiduk
et al. in Gaiduk  _et al._ (2009); Gaiduk and Staroverov (2010); Gaiduk
(2013). The line integration formula has been also discussed by the authors of
Ref. Karasiev and López-Boada (1998), as we shall mention below.

The novel contribution of this paper consists in outlining a possible path
towards the full solution of the IKS problem, and of its implementation for a
simple and yet realistic case, in nuclear physics. In particular, we have
merged the solution of the first branch of the IKS problem (as already
obtained in Accorto  _et al._ (2020)) with the line integration formula
proposed in Ref. van Leeuwen and Baerends (1995). We will demonstrate that, by
choosing a path of conveniently scaled densities $\rho_{t}({\bf r})$,
parametrized by a real number $t$, and by inserting densities and associated
potentials in the line integration formula, we can reconstruct the total
energy of the system, relative to the g.s. energy, as a function of the
scaling parameter $t$ with very good accuracy. This will allow to gain some
insights, as we shall discuss, about the functional form of the underlying
EDF.

Our method shows some resemblance with the scheme proposed by Dobaczewski et
al. in Refs. Dobaczewski (2016); Salvioni  _et al._ (2020). Indeed, both
approaches share the idea that information about nuclei driven outside their
ground state should be used in order to probe the exact EDF. To do that, some
external perturbation must be introduced. In Ref. Salvioni  _et al._ (2020)
extensive ab initio calculations have been performed, where both one- and two-
body operators were added to the Hamiltonian. Our approach differs in that we
generate a set of benchmark densities by means of a perturbing one-particle
external potential $V_{\rm ext}(r)$, simplifying in that way the required
benchmark calculations for the IKS method.

In order to appreciate the relevance of our work, a few words about the status
of nuclear DFT are in order. Nuclear DFT is a more challenging problem than
electronic DFT. In the case of electronic systems, the Coulomb interaction is
well-known and when building an EDF it is natural to single out the
contribution from the direct Hartree term. This, together with the kinetic
energy, provides the largest contribution to the EDF and leaves out only the
so-called exchange-correlation part to be determined. In contrast, the
effective Hamiltonian that describes the nucleon-nucleon interaction is not so
accurately known (an overview of ab initio nuclear physics can be found in
Refs. Hergert (2020); Hjorth-Jensen  _et al._ (2017); Leidemann and Orlandini
(2013)). A further issue is that three nucleon forces are relevant and cannot
be ignored.

In spite of all these difficulties, nuclear DFT has achieved a considerable
number of successes. These are described in the recent textbook Schunck
(2019). A shorter review can be found in Ref. Colò (2020). Throughout the
isotope chart, state-of-the-art nuclear EDFs can describe the overall trend of
masses (with about 1 MeV average accuracy) and nuclear radii (with about few
times 10-2 fm average accuracy), as well as some other properties like the
electric quadrupole moments or the main features of giant resonances. Hundreds
of EDFs exist because many possible functional forms can be assumed, and many
protocols to determine the EDF parameters can be designed. It is far from
obvious that new attempts to extend the EDF forms or to employ more refined
fitting protocols may lead to substantial improvements Kortelainen  _et al._
(2014). Current EDFs are frequently based on nuclear phenomenology and suffer
from large extrapolation errors Erler  _et al._ (2012), which are very
difficult to quantify. In other words, a well established theoretical scheme
for a systematic improvement does not seem to exist at present. For those
reasons, some practitioners believe that a re-thinking of the nuclear DFT
strategy would be timely Kortelainen  _et al._ (2014); Furnstahl (2020);
Dobaczewski (2016) (see also the work of some of us in Ref. Marino  _et al._
(2021)).

In this respect, non-empirical EDFs, or even strategies to better constrain
some terms of the EDFs, would be extremely welcome. Deriving the nuclear EDF
from ab initio is a possible avenue and is the subject of intense
investigations Salvioni  _et al._ (2020); Marino  _et al._ (2021); Zurek  _et
al._ (2021). Solving the IKS problem with experimental or ab initio densities
as an input is a complementary possibility. Our work has the goal of showing
that, at least in a simple situation, the IKS problem can indeed be solved.

In Sec. II, the KS framework is reviewed and the line integration formula is
presented. Then, in Sec. III, a few key aspects related to the use of this
formula are discussed in detail. Sec. IV provides a description of the method
we employ to perform the D2P inversion, that is, the constrained variational
(CV) method. In Sec. V, our method is applied to a case study. Lastly,
conclusions and perspectives of our work are outlined in Sec. VI.

##  II The line integration formula

We first remind the essential notions of DFT and clarify our notation. Useful
references are Parr and Yang (1994); Colò (2020); Schunck (2019).

In the KS scheme, the EDF is written as in Eq. (2) as a functional of the
auxiliary s.p. orbitals. $T$ is the non-interacting kinetic energy.

From the universal part $F$, that ultimately stems from all the inter-particle
interactions, we define the potential $v$ (also named self-consistent
potential in what follows),

$v([\rho];{\bf r})=\frac{\delta F}{\delta\rho({\bf r})}.$ (6)

On the other hand, the external potential $v_{\rm ext}({\bf r})$ and the
external contribution to the total energy are related by

$\displaystyle V_{\rm ext}=\int d^{3}r\ v_{\rm ext}({\bf r})\rho({\bf r})$ (7)
$\displaystyle v_{\rm ext}({\bf r})=\frac{\delta V_{\rm ext}}{\delta\rho({\bf
r})}.$ (8)

Applying the variational principle to Eq. (2) leads to the following KS
equations for the s.p. orbitals:

$\displaystyle\left(-\frac{\hbar^{2}}{2m}\nabla^{2}+v_{\rm KS}\right)\phi_{i}$
$\displaystyle=\epsilon_{i}\phi_{i},$ (9) $\displaystyle v_{KS}$
$\displaystyle=v[\rho]+v_{ext}.$ (10)

The KS potential $v_{KS}$ is the sum of the self-consistent potential and the
external term. The latter is a function of the position only, while the former
also depends on the density.

After setting the framework, we can now focus on the relation between the
self-consistent potential and the EDF, namely the aforementioned line
integration formula. In Ref. van Leeuwen and Baerends (1995) (cf. also Gaiduk
_et al._ (2009)), it is shown that if one knows the effective potential
$v[\rho]$ along a path of densities, then the corresponding change in the
energy functional can be reconstructed. In particular, a one-parameter family
of densities is considered. Accordingly, we shall write the densities as
$\rho_{t}(\bf r)$, with $A\leq t\leq B$. The reconstruction formula has been
discussed in van Leeuwen and Baerends (1995); Gaiduk  _et al._ (2009), which
focus on electronic DFT, for the exchange-correlation part of the functional
only. On the same grounds, we can write the following formula for $F$:

$F[\rho_{B}]-F[\rho_{A}]=\int_{A}^{B}dt\ \int d^{3}r\ v\left([\rho_{t}({\bf
r})],\bf r\right)\derivative{\rho_{t}({\bf r})}{t}\ .$ (11)

This formula is actually holds for any functional of the density and of its
gradients. We also remind that EDFs are generally written as the integral of
an energy density $f$

$F[\rho]=\int d^{3}r\ f(\bf r,\rho,\nabla\rho,\ldots).$

which is unique up to a gauge transformation, i.e. a term $\theta(\mathbf{r})$
whose volume integral vanishes.

For later convenience, we also define $I_{t}(R)$ by means of

$I_{t}(R)=\int_{0}^{R}drr^{2}\int d\Omega\,v\left([\rho_{t}({\bf r})],\bf
r\right)\derivative{\rho_{t}({\bf r})}{t},$ (12)

where $\Omega$ is the solid angle and $R$ is the radial upper bound for the
numerical integration. We expect $I_{t}(R)$ to be a convergent function as $R$
goes to infinity. Then, (11) can also be written as

$F[\rho_{B}]-F[\rho_{A}]=\int_{A}^{B}dt\ I_{t}(R\longrightarrow+\infty).$ (13)

In principle, any reasonable density path could be chosen in order to perform
the line integral. Three specific paths were discussed in Ref. Gaiduk  _et
al._ (2009). Our actual choice of the density path shall be discussed below.

The line integration formula, so far, has been mostly employed in cases where
an analytical dependence of the effective potential with respect to the
density is given, see Gaiduk  _et al._ (2009). The same perspective is taken
by the Levy-Perdew virial relation Levy and Perdew (1985); Ghosh and Parr
(1985) between exchange potential and exchange energy, where the energy can be
deduced from the potential evaluated for just one density. In passing, we
mention that line integration techniques have found applications also in
constructing approximation to the non-local kinetic energy in orbital free DFT
Xu  _et al._ (2019); Mi  _et al._ (2018).

Our idea is to move one step forward: the formula (11) shall be applied to
cases where such prior analytic knowledge $v[\rho]$ is not available.
Specifically, the relation between densities and potentials shall be defined
by a numerical inversion procedure, as detailed in Sec. IV.

##  III Conceptual remarks

We present here some remarks that concern the validity of the line integration
formula.

###  III.1 v-representable solutions

The first question concerns the choice of the density path. While it is simple
to design a path $\rho_{t}(\mathbf{r})$ that includes only $N$-representable
densities, it is not possible to guarantee, in general, that the path avoids
densities that are not $v$-representable Kryachko and Ludeña (2014). From a
physical viewpoint, a reasonable and somehow conservative choice may be that
of spanning a path of densities that are close to the actual ground state
density of a given system or to a realistic approximation thereof, e.g., as it
may be determined by an ab initio calculation. A theorem by Kohn Kohn (1983)
gives theoretical support to our argument.

As an example, in Sec. V results will be presented for densities originating
from monopole deformations of a nucleus Garg and Colò (2018). In this case,
the existence of a wide literature on the effects of small variations of the
shape of a nucleus around its equilibrium state can be exploited to design
clever and physically motivated density paths. To sum up, we are confident
that a sound physical intuition of the systems at hand allows to avoid the
$v$-representability problem.

###  III.2 Self-consistent potential

The second question is a key issue for our discussion. The D2P solution allows
us to determine the KS potential of (10), that is,

$v_{KS}=v[\rho]+v_{ext}.$

However, Eq. (11) specifically requires the knowledge of $v[\rho]$ alone.
While we cannot suggest a general solution, we propose a way to circumvent
this difficulty in specific cases. Indeed, the freedom in the choice of the
density path has to be exploited. Theoretical calculations of the systems
under study, subject to a known external potential $v_{ext,t}$, can be
performed and the corresponding densities $\rho_{t}(\mathbf{r})$ can be taken
as input. In Sec. V, we will show a demonstration of the feasibility of this
approach where constrained Hartree-Fock (CHF) calculations of an atomic
nucleus Ring and Schuck (1980) are employed to produce the benchmark
$\rho_{t}(\mathbf{r})$ from a given $v_{ext,t}$. In the longer term, it would
be interesting to employ ab initio simulations as a benchmark.

###  III.3 Energy conservation

In this context, an interesting identity that has been first introduced in
Ref. Karasiev and López-Boada (1998), is worth discussing, because it is an
useful test of our numerical procedure. Let us assume that
$\rho_{t}(\mathbf{r})$ is a number-preserving path, so that all densities are
normalized to the same number of particles. Let us consider the fundamental
DFT Euler equation Parr and Yang (1994)

$\functionalderivative{T[\rho]}{\rho(\mathbf{r})}+v_{KS}([\rho],\mathbf{r})=\mu[\rho],$
(14)

where $\mu$ is the chemical potential. Then, we can multiply Eq. (14) by
$\derivative{\rho_{t}}{t}$ and integrate over both $t$ and $\mathbf{r}$. Since
$\mu$ is independent from $\mathbf{r}$ and $\int
d\mathbf{r}\rho_{t}(\mathbf{r})$ is a constant with respect to $t$, the r.h.s.
vanishes. It follows immediately that

$\Delta T=-\int_{A}^{B}dt\int d\mathbf{r}\,v_{KS}\left(\left[\rho_{t}(\bf
r)\right],\bf r\right)\derivative{\rho_{t}(\mathbf{r})}{t}.$ (15)

Note that the right hand side of this equation is different from that of Eq.
(11) and it is equal to $-\Delta U_{KS}=-[\Delta F+\Delta U_{ext}]$. The total
variation of the energy $\Delta T+\Delta U_{KS}=0$ is conserved as it should
since the chosen path comes from a minimization of Eq. (2) for each value of
$t$. Therefore, this relation is a conservation law that holds on density and
potential paths $\rho_{t}$ and $v_{KS}\left[\rho_{t}\right]$. It is
potentially well suited to test the correctness of the implementation of both
the D2P step and Eq. (11). Indeed, as detailed in Sec. IV below, from an IKS
calculation both the KS potential and the kinetic energy are obtained. One can
then compare $\Delta T=T[\rho_{B}]-T[\rho_{A}]$ to the integral on the r.h.s.
of Eq. (15), thus gaining information on the relative numerical accuracy of
the two procedures (Sec. V).

##  IV Constrained-variational method

In this Section, the formalism and the implementation of the CV method Accorto
_et al._ (2020); Jensen and Wasserman (2018) are reviewed. In the KS scheme
Parr and Yang (1994), a system of $N$ independent particles is first
considered, which is described by a set of single-particle (s.p.) orbitals
$\phi_{\alpha}(\mathbf{r})$. The independent-particle kinetic energy is then a
functional of the orbitals and reads

$\displaystyle
T\left[\\{\phi_{\alpha}(\mathbf{r})\\}\right]=-\frac{\hbar^{2}}{2m}\sum_{\alpha=1}^{N}\int
d\mathbf{r}\
\phi_{\alpha}^{*}(\mathbf{r})\nabla^{2}\phi_{\alpha}(\mathbf{r}).$ (16)

Suppose that the ground state (g.s.) density of the system $\rho(\mathbf{r})$
(1) be equal to a given target density $\tilde{\rho}(\mathbf{r})$. Such
density may be known from experiment or an ab initio calculation. Under these
assumptions, the g.s. is characterized by the condition that the kinetic
energy (as a functional of the orbitals) is minimized under the constraint
$\rho(\mathbf{r})=\tilde{\rho}(\mathbf{r})$. Moreover, the additional
orthonormality constraints $\int d\mathbf{r}\
\phi_{\alpha}^{*}(\mathbf{r})\phi_{\beta}(\mathbf{r})=\delta_{\alpha\beta}$
are imposed on the orbitals. The problem has the form of a constrained
variation of an integral functional. The Lagrange multipliers method allows to
convert it into the free optimization of the objective functional

$\displaystyle
J\left[\\{\phi_{\alpha}(\mathbf{r})\\};v_{KS}(\mathbf{r}),\\{\epsilon_{\alpha\beta}\\}\right]=T\left[\\{\phi_{\alpha}(\mathbf{r})\\}\right]+$
$\displaystyle\hskip 85.35826pt+\int d\mathbf{r}\
v_{KS}(\mathbf{r})\left[\rho(\mathbf{r})-\tilde{\rho}(\mathbf{r})\right]$
$\displaystyle\hskip
85.35826pt-\sum_{\alpha<\beta}\epsilon_{\alpha\beta}\left(\int d\mathbf{r}\
\phi_{\alpha}^{*}(\mathbf{r})\phi_{\beta}(\mathbf{r})-\delta_{\alpha\beta}\right)$
. (17)

We stress that $J$ is a functional of the additional variables
$v_{KS}(\mathbf{r})$ and of $\epsilon_{\alpha\beta}$. We have readily
identified the Lagrange multiplier associated to the density constraint with
the KS potential $v_{KS}$. These two quantities coincide only if the minimum
of $J$ has been obtained.

By solving $\delta J=0$ plus the subsidiary conditions, the value of the
orbitals and of the multipliers can be determined. The Euler-Lagrange
equations for $J$ read

$\displaystyle-\frac{\hbar^{2}}{2m}\nabla^{2}\phi_{\alpha}(\mathbf{r})+v_{KS}(\mathbf{r})\phi_{\alpha}(\mathbf{r})=\sum_{\beta}\epsilon_{\alpha\beta}\phi_{\beta}(\mathbf{r}).$
(18)

These are s.p. Schrödinger equations in a non-canonical form Karasiev and
López-Boada (1998). After these have been solved, a diagonalization of the
$\epsilon$ matrix allows to determine the canonical eigenfunctions and the
s.p. energies, which satisfy

$-\frac{\hbar^{2}}{2m}\nabla^{2}\phi_{\alpha}(\mathbf{r})+v_{KS}(\mathbf{r})\phi_{\alpha}(\mathbf{r})=\lambda_{\alpha}\phi_{\alpha}(\mathbf{r}).$
(19)

Thus, Eq. (19) confirms the physical interpretation of the Lagarnge
multipliers.

In practice, works on the IKS problem have focused on spherically symmetric
systems Accorto  _et al._ (2020, 2021). As a consequence, the problem
simplifies to a one-dimensional radial equation, and the labels $\alpha,\beta$
will correspond to the usual quantum numbers $(nlj)$ with the standard shell-
model ordering. The angular and spin parts of the wave functions are then
known, and only the radial functions $u_{\alpha}(r)$ have to be determined.
The orthonormality conditions have to be imposed only if
$l_{\alpha}=l_{\beta}$ and $j_{\alpha}=j_{\beta}$ and read
$\int_{0}^{+\infty}dr\
u_{\alpha}(r)u_{\beta}(r)=\delta_{n_{\alpha}n_{\beta}}$. We note here that the
inclusion of a spin-orbit potential is left for future developments. The KS
potential is then orbital-independent.

The numerical implementation is based on a discretization scheme, where the
radial functions are evaluated on a radial mesh and the finite difference
method is used to approximate the derivatives Jensen and Wasserman (2018). The
solution proceeds as follows. First, the constrained optimization library
IPOPT Waechter (2009); Wächter and Biegler (2005) is used to find the orbitals
$u_{\alpha}(r)$ that minimize the kinetic energy and at the same time satisfy
the aforementioned constraints within a given tolerance (see also Ref. Accorto
_et al._ (2020)). Next, we aim at determining the KS potential and the energy
eigenvalues. The previously determined radial functions are plugged into the
non-canonical Schrödinger equations (18) that now turn into a linear system,
whose unknowns are $v_{KS}(r)$ on each point of the mesh and the
$\epsilon_{\alpha\beta}$’s Jensen and Wasserman (2018). This is an
overdetermined algebraic system, with one equation for each point of the mesh
and each different radial function. A least-squares algorithm is employed to
solve it. Lastly, the $\epsilon$ matrix is diagonalized and the canonical
solutions of (19) are thus found. The outcome of the CV method, therefore,
consists in the orbitals $u_{\alpha}(r)$, the non-interacting kinetic energy
$T$, the KS potential $v_{KS}(r)$ and the eigenvalues $\lambda_{\alpha}$ of
the effective s.p. description of the system.

A new code has been developed in Python which allows to perform the D2P
inversion once a density is given as input. A technical note on the numerical
inputs of the CV calculations is here in order. We have found that the
convergence of the method is mostly affected by the following IPOPT
parameters: the relative tolerance on the value of the objective function and
the absolute tolerance on the fulfillment of the constraints Waechter (2009).
The choice of the upper bound $R$ of the radial mesh is also important.

To be physically meaningful, a radial potential should vanish as
$r\longrightarrow+\infty$. It is well-known that the asymptotic behaviour of
the target density directly impacts that of the corresponding potential. For
example, Gaussian tails are notoriously problematic, cf. Accorto  _et al._
(2020); Kanungo  _et al._ (2019). The densities employed in this work,
however, are all well-behaved, as far as their asymptotic properties are
concerned. A further caveat relates to the fact that a potential is defined
only up to a constant. Thus both $v(r)$ and $\lambda_{\alpha}$ must be
shifted, so as to ensure that $v(r)$ goes to zero at large $r$. Actually, the
correct asymptotic behaviour of the potential is ensured by a correct
reproduction of the ionization energy in atomic physics or neutron/proton
separation energy in nuclear physics Perdew  _et al._ (1982); Perdew and Levy
(1997). This is a further check or recipe that should be applied whenever this
information is available.

Lastly, we mention that some inaccuracies may occasionally show up near the
border of the radial mesh. For example, the potential may take unnatural
values in a limited number of points in such region. However, we have verified
this has no impact on the energy integrals due to the exponentially decreasing
behaviour of the densities.

##  V Results

###  V.1 Physical case and numerical details

In this Section, we will implement the strategy outlined in the previous
Section and solve the two-step IKS problem. Namely, we will try to reconstruct
the functional starting from a set of densities $\rho_{\mu}$, with
$\mu_{A}<\mu<\mu_{B}$. The densities will be generated by means of a set of KS
calculations carried out with a Skyrme EDF, to which a harmonic external
potential $v_{ext,\mu}=\mu r^{2}$ is added. Note that in the nuclear physics
context these are often called self-consistent (constrained) Hartree-Fock
calculations.

For each value of $\mu$, we will use the procedure of Sec. IV to obtain the KS
potential $v([\rho_{\mu}],r)$ by inversion (first step, D2P). We will then
insert these potentials in Eq. (11) to reconstruct the energy difference
$\Delta F=F[\rho_{\mu_{B}}]-F[\rho_{\mu_{A}}]$ (second step, P2E). At the end,
we will assess to what accuracy the procedure reproduces the benchmark energy
difference associated with the original Skyrme EDF.

We highlight that the same procedure can be followed starting from a set of
densities generated by ab initio calculations, without making reference to any
pre-existing EDF. The two-step solution of the IKS problem can then lead to
values of $\Delta F$ that will provide precise information on the underlying
EDF.

For this first application of the method, we will present results obtained
with a simplified EDF. In particular, we benchmark our method to a
$t_{0}-t_{3}$ Skyrme EDF (see e.g. Ref. Chabanat  _et al._ (1997)) with which
we generate the target densities. This EDF is characterized by the energy
density

$F=\int d^{3}r\,\,\mathcal{F}=\int d^{3}r\,(\mathcal{F}_{0}+\mathcal{F}_{1}),$
(20)

where

$\displaystyle\mathcal{F}_{0}=\frac{1}{4}t_{0}(2\rho^{2}-(\rho_{p}^{2}+\rho_{n}^{2})),$
(21)
$\displaystyle\mathcal{F}_{1}=\frac{1}{24}t_{3}\,\rho^{\alpha}(2\rho^{2}-(\rho_{p}^{2}+\rho_{n}^{2})),$
(22) $\displaystyle\rho=\rho_{p}+\rho_{n}.$ (23)

Since $F$ depends on the densities, but not on their gradients, the potential
takes the form

$v_{q}([\rho({\bf r})],{\bf r})=\frac{\delta
F}{\delta\rho_{q}}=\frac{\partial\mathcal{F}}{\partial\rho_{q}},$ (24)

where $q=n,p$. The EDF parameters are reported in Tab. 1 for a few choices of
the parameter $\alpha$. Once $\alpha$ has been chosen, the coefficients
$t_{0}$ and $t_{3}$ are determined in such a way that the equation of state of
symmetric nuclear matter predicted by the EDF has a minimum at the empirical
saturation point, $\rho=$0.16 fm-3 and $E=-$16 MeV. In Sec. V.2 and V.3, we
test our method by employing the EDF with $\alpha=$0.32. To further simplify
the problem, we neglect Coulomb and spin-orbit and we pick up the
representative case of the 16O densities ($\rho_{n}=\rho_{p}=\rho/2$).

$\alpha$ |  $t_{0}$ (MeV fm3) |  $t_{3}$ (MeV fm3+3α)
---|---|---
0.16 |  -3068.63 |  19578.08
0.2 |  -2581.88 |  16853.70
0.32 |  -1851.75 |  13124.44
1 |  -1024.27 |  14602.57
Table 1: Parameters of the $t_{0}-t_{3}$ EDFs.

A family of scaled densities are generated by means of constrained Hartree-
Fock (CHF) calculations Ring and Schuck (1980); Colò  _et al._ (2013). As a
perturbation, a harmonic external potential $v_{ext,\mu}(r)=\mu r^{2}$, with
$\mu$ in the range [-0.25,0.25] MeV$\cdot$ fm-2 is employed. The physical
meaning is that of driving a scaling of the radius (and of the whole nuclear
density) with respect to the unperturbed case. This approach has its original
motivation in the study of monopole deformations Garg and Colò (2018). The
true ground state of the system is given by $\mu=0$ and is an absolute minimum
of the energy $T+F$.

The D2P inversion is performed by means of the CV method described in Sec. IV,
which provides both the effective KS potential $v_{KS}$ and the kinetic energy
associated to the neutron density.

###  V.2 Potentials

To exemplify the subtle point concerning the difference between $v_{KS}$ and
$v[\rho]$ (Sec. III), in Fig. 1 and 2 the CHF potential for either $\mu=-0.2$
or $\mu=0.2$ Mev fm-2 is compared to the IKS potential (KS). The two functions
are clearly different, while the third curve
$v[\rho_{\mu}]=v_{KS,\mu}-v_{ext,\mu}=v_{KS,\mu}-\mu r^{2}$, defined by
subtracting the harmonic term from the IKS potential, matches rather well the
CHF potential. This proves the accuracy of the CV method.

From now on, only the self-consistent part of the potential, $v[\rho_{\mu}]$,
shall be displayed. In Fig. 3, three densities (top) and the corresponding
potentials (bottom), determined by means of the CV inversion, are compared for
three different values of $\mu$. A positive value $\mu>0$ acts as a confining
potential and as a consequence leads to a more compact density profile than in
the unperturbed case. Conversely, a repulsive external potential ($\mu<0$)
leads to systems which are more spread out. Consequently, in the latter case
the density peak in the interior of the nucleus is less pronounced.

Figure 1:  The potential derived from $\rho_{\mu}$ by the IKS procedure
($v_{KS}$), the CHF potential $v_{CHF}$ and $v=v_{KS,\mu}-\mu r^{2}$ are
compared for $\mu=-0.2$ MeV$\cdot$ fm-2.  Figure 2: Same as Fig. 1, but for
$\mu=0.2$ MeV$\cdot$ fm-2.  Figure 3: 16O neutron densities (top) and
corresponding self-consistent potentials (bottom) for three values of the
perturbation strength $\mu$ (in MeV fm-2).

###  V.3 Line integration

We can now move on to the study of the P2E and the line integration formula. A
preliminary test is presented in Fig. 4. There, the function $I_{\mu}(R)$ (12)
is plotted as a function of the radius $R$ in the cases $\mu=-$0.2, 0, 0.2
MeV$\cdot$ fm-2. Our concern here is that of verifying the asymptotic
convergence of $I_{\mu}$ to a constant for large $R$. Indeed, the convergence
is quite fast and a stable result is reached already for $R=4$ fm. We remind
here that 16O has a root mean square radius of about 3 fm and, thus, $I_{\mu}$
is expected to converge for $R$ values that scale with $A^{1/3}$.

Next, we can turn to the study of the kinetic energy and the potential energy
contributions separately. In Fig. 5 different behaviours of the kinetic energy
differences, $\Delta T(\mu)=T(\mu)-T(\mu=0)$, are shown as functions of $\mu$.
We consider two independent ways of evaluating $\Delta T(\mu)$. Indeed, one
can simply subtract the kinetic energies provided by the CV method when
performing inversion at $\mu=0$ and at finite $\mu$. Otherwise, the formula by
Karasiev et al. Karasiev and López-Boada (1998), Eq. (15), can be exploited,
where the kinetic energy is computed from the knowledge of the potential
$v_{KS,\mu}$ obtained in the CV framework on a grid of intermediate values of
$\mu$. The two estimates of the kinetic energy difference, labelled as IKS and
Eq. (15) in the figure, are compared to the CHF results. The top panel
highlights that there is a very good agreement between the three curves for
all values of $\mu$. The deviations from the CHF results are plotted in the
bottom panel. It can be appreciated that the discrepancies are quite small,
being of the order of 10-2 MeV. In comparison, the neutron kinetic energy of
the unperturbed system amounts to $T_{\mu=0}=148.37$ MeV (including the center
of mass correction). We observe, however, that Eq. (15) guarantees an overall
better accuracy. Therefore, it will be used in the following to estimate
$\Delta T$.

Figure 4:  $I_{\mu}(R)$ (defined in (12)) as a function of the integration
limit $R$ for three representitative values of $\mu$. The integral function
converges asymptotically to a constant for large enough $R$.  Figure 5: Top:
kinetic energy differences $\Delta T(\mu)=T(\mu)-T(\mu=0)$ obtained with CHF
calculations, by inversion of the densities $\rho_{\mu}(r)$ (IKS), and by
employing Eq. (15). Bottom: discrepancy between the CHF and both IKS and Eq.
(15) results.

We now compute the universal term $F$ or, more precisely, $\Delta
F(\mu)=F(\mu)-F(\mu=0)$, by the line integration method introduced in Sec. II
according to Eq. (11). Fig. 6 shows that the reconstructed IKS curve is in
excellent agreement with the CHF values obtained with the original Skyrme EDF,
with discrepancies of about $10^{-3}$ MeV. This is rather encouraging, since
it implies that the line integration has been implemented to a high numerical
precision. At the same time, it is a strong support to the reliability of the
IKS machinery. This is confirmed by Fig. 7, where we display the total energy
difference $\Delta F(\mu)+\Delta T(\mu)$, accurate up to a fraction of keV.

Figure 6: Top: universal energy differences $\Delta F(\mu)=F(\mu)-F(\mu=0)$
obtained with CHF calculations, by inversion of the densities $\rho_{\mu}(r)$
(IKS). Bottom: discrepancy between the IKS and CHF results.  Figure 7: Top:
total energy $\Delta F(\mu)+\Delta T(\mu)$ obtained with IKS and CHF
calculations. Bottom: discrepancy between the IKS and CHF results.

A question now arises: how much information have we actually gained about the
original EDF, by solving the two-step inverse IKS problem? We will limit
ourselves to consider the family of $t_{0}-t_{3}$ Skyrme EDFs, obtained by
varying the value of the exponent $\alpha$ (cf. Tab. 1) and by assuming that
they are determined in symmetric nuclear matter, as discussed above (cf. Tab.
1). In other words, we will not discuss how to determine the absolute
strengths ($t_{0}$ and $t_{3}$) of the different terms in the underling
$F[\rho]$, but only if we can be sensitive to the exponent of the density
dependence. To this aim, we compare in Fig. 8 the values obtained for $\Delta
T+\Delta F$ as a function of $\mu$ for a few different values of $\alpha$. The
resolving power of the method will depend on the range of $\mu$ taken into
account in the calculation. Nonetheless, the figure already shows that, for a
modest change in $\mu$, our benchmark value $\alpha=0.32$ can be rather
clearly distinguished from $\alpha$=0.16 and $\alpha$=1. This type of
information may become instrumental in building new EDFs based on ab inito
calculations.

Figure 8:  The total energy $\Delta F(\mu)+\Delta T(\mu)$ reconstructed with
the two-step IKS method and associated with the value $\alpha=0.32$ is
compared to the values obtained with $t_{0}-t_{3}$ EDFs characterized by
different values of $\alpha$.

##  VI Conclusions and perspectives

A complete solution to the inverse Kohn-Sham problem of Density Functional
Theory has been proposed. First, the effective KS potential associated to a
given ground state density is determined. Second, a path of perturbed
densities is chosen, and the knowledge of the associated KS potentials is
exploited to compute the difference between the energies of the perturbed and
unperturbed states.

Benchmark calculations have been performed in a case relevant for nuclear
physics. In particular, we have shown that, for a simple $t_{0}-t_{3}$ Skyrme
EDF, and perturbing the system with an external harmonic potential, this
method allows to reconstruct the energies with a rather good numerical
accuracy.

These promising results open up a number of perspectives. First, one should
apply our method using ab initio calculations as an input. The feasibility of
these calculations when the nuclear systems are perturbed by external
potentials, and the assessment of the associated theoretical uncertainties,
are issues that are mandatory to clarify.

In principle, accurate ab initio calculations of nuclear densities perturbed
by a variety of external potentials contain a lot of useful information to
improve existing EDFs. In practice, one could apply the IKS methodology
outlined in the present work to produce a set of ab initio metadata on how
energies change under the action of perturbations.

In our work, we have not tackled the problem of how to fully determine the EDF
from these metadata. When ab initio inputs are available, we envisage to use
also the resulting absolute values of the energy. Implementing appropriate
fitting algorithms on these absolute energies may allow obtaining a set of
novel constraints on the different terms of the EDF. As the convergence of
this procedure is not guaranteed, because the choice of the different terms is
the result of an ansatz, and the EDF which is compatible with the ab initio
input may not correspond to that ansatz, one could devise as an alternative a
machine learning algorithm (see, e.g., Schmidt  _et al._ (2019)). This should
be able to select the terms in the EDF that are more likely to appear, and
quantitatively estimate their relevance.

In a later stage, further studies could be envisaged by considering a set of
different external potentials that are coupled to spin densities, isospin
densities, or other types of densities. This investigation could shed light on
the terms of the nuclear EDF that are not merely sensitive to the total
density.

## References

  * Jensen and Wasserman (2018) D. S. Jensen and A. Wasserman, International journal of quantum chemistry 118 (2018).
  * Shi and Wasserman (2021) Y. Shi and A. Wasserman, The journal of physical chemistry letters 12, 5308 (2021).
  * Parr and Yang (1994) R. G. Parr and W. Yang,  _Density-Functional Theory of Atoms and Molecules (International Series of Monographs on Chemistry)_ (Oxford University Press, USA, 1994).
  * Colò (2020) G. Colò, Advances in Physics: X 5, 1740061 (2020).
  * Runge and Gross (1984) E. Runge and E. K. U. Gross, Phys. Rev. Lett. 52, 997 (1984).
  * Görling (1992) A. Görling, Phys. Rev. A 46, 3753 (1992).
  * Wang and Parr (1993) Y. Wang and R. G. Parr, Phys. Rev. A 47, R1591 (1993).
  * van Leeuwen and Baerends (1994) R. van Leeuwen and E. J. Baerends, Phys. Rev. A 49, 2421 (1994).
  * Zhao  _et al._ (1994) Q. Zhao, R. C. Morrison, and R. G. Parr, Phys. Rev. A 50, 2138 (1994).
  * Kanungo  _et al._ (2019) B. Kanungo, P. M. Zimmerman, and V. Gavini, Nature communications 10, 4497 (2019).
  * Naito  _et al._ (2019) T. Naito, D. Ohashi, and H. Liang, Journal of Physics B: Atomic, Molecular and Optical Physics 52, 245003 (2019).
  * Nam  _et al._ (2020) S. Nam, S. Song, E. Sim, and K. Burke, Journal of chemical theory and computation 16, 5014 (2020).
  * Nam  _et al._ (2021) S. Nam, R. J. McCarty, H. Park, and E. Sim, The Journal of Chemical Physics 154, 124122 (2021).
  * Pérez-Jordá (2020) J. M. Pérez-Jordá, International Journal of Quantum Chemistry 120, e26052 (2020).
  * Accorto  _et al._ (2020) G. Accorto, P. Brandolini, F. Marino, A. Porro, A. Scalesi, G. Colò, X. Roca-Maza, and E. Vigezzi, Phys. Rev. C 101, 024315 (2020).
  * Accorto  _et al._ (2021) G. Accorto, T. Naito, H. Liang, T. Nikšić, and D. Vretenar, Phys. Rev. C 103, 044304 (2021).
  * Macocco (2015) I. Macocco, “Costruzione del potenziale nucleare efficace di Kohn-Sham dalle densità empiriche,”  (2015), Bachelor Thesis in Italian.
  * Kumar  _et al._ (2019) A. Kumar, R. Singh, and M. K. Harbola, Journal of Physics B: Atomic, Molecular and Optical Physics 52, 075007 (2019).
  * Andrade  _et al._ (2015) X. Andrade, D. Strubbe, U. De Giovannini, A. H. Larsen, M. J. T. Oliveira, J. Alberdi-Rodriguez, A. Varas, I. Theophilou, N. Helbig, M. J. Verstraete, L. Stella, F. Nogueira, A. Aspuru-Guzik, A. Castro, M. A. L. Marques, and A. Rubio, Physical chemistry chemical physics : PCCP 17, 31371 (2015).
  * Garrigue (2021) L. Garrigue, “Building Kohn-Sham potentials for ground and excited states,” (2021), arXiv:2101.01127 [math-ph] .
  * Umrigar and Gonze (1994) C. Umrigar and X. Gonze,  _High Performance Computing and Its Applications in The Physical Sciences_ (World Scientific, 1994).
  * Schmidt  _et al._ (2019) J. Schmidt, C. L. Benavides-Riveros, and M. A. L. Marques, The Journal of Physical Chemistry Letters 10, 6425 (2019).
  * van Leeuwen and Baerends (1995) R. van Leeuwen and E. J. Baerends, Phys. Rev. A 51, 170 (1995).
  * Gaiduk  _et al._ (2009) A. P. Gaiduk, S. K. Chulkov, and V. N. Staroverov, Journal of Chemical Theory and Computation 5, 699 (2009).
  * Gaiduk and Staroverov (2010) A. P. Gaiduk and V. N. Staroverov, The Journal of Chemical Physics 133, 101104 (2010).
  * Gaiduk (2013) A. P. Gaiduk,  _Theory of Model Kohn-Sham Potentials and its Applications_ , Ph.D. thesis, University of Western Ontario, Canada (2013).
  * Karasiev and López-Boada (1998) V. Karasiev and R. López-Boada, Phys. Rev. A 58, 1954 (1998).
  * Dobaczewski (2016) J. Dobaczewski, Journal of Physics G: Nuclear and Particle Physics 43, 04LT01 (2016).
  * Salvioni  _et al._ (2020) G. Salvioni, J. Dobaczewski, C. Barbieri, G. Carlsson, A. Idini, and A. Pastore, Journal of Physics G: Nuclear and Particle Physics 47, 085107 (2020).
  * Hergert (2020) H. Hergert, Frontiers in Physics 8 (2020).
  * Hjorth-Jensen  _et al._ (2017) M. Hjorth-Jensen, M. Lombardo, and U. van Kolck,  _An Advanced Course in Computational Nuclear Physics: Bridging the Scales from Quarks to Neutron Stars_ (2017).
  * Leidemann and Orlandini (2013) W. Leidemann and G. Orlandini, Progress in Particle and Nuclear Physics 68, 158 (2013).
  * Schunck (2019) N. Schunck, ed.,  _Energy Density Functional Methods for Atomic Nuclei_, 2053-2563 (IOP Publishing, 2019).
  * Kortelainen  _et al._ (2014) M. Kortelainen, J. McDonnell, W. Nazarewicz, E. Olsen, P.-G. Reinhard, J. Sarich, N. Schunck, S. M. Wild, D. Davesne, J. Erler, and A. Pastore, Phys. Rev. C 89, 054314 (2014).
  * Erler  _et al._ (2012) J. Erler, N. Birge, M. Kortelainen, W. Nazarewicz, E. Olsen, A. M. Perhac, and M. Stoitsov, Nature (London) 486, 509 (2012).
  * Furnstahl (2020) R. J. Furnstahl, The European Physical Journal A 56 (2020).
  * Marino  _et al._ (2021) F. Marino, C. Barbieri, A. Carbone, G. Colò, A. Lovato, F. Pederiva, X. Roca-Maza, and E. Vigezzi, Phys. Rev. C 104, 024315 (2021).
  * Zurek  _et al._ (2021) L. Zurek, E. A. Coello Pérez, S. K. Bogner, R. J. Furnstahl, and A. Schwenk, Phys. Rev. C 103, 014325 (2021).
  * Levy and Perdew (1985) M. Levy and J. P. Perdew, Phys. Rev. A 32, 2010 (1985).
  * Ghosh and Parr (1985) S. K. Ghosh and R. G. Parr, The Journal of Chemical Physics 82, 3307 (1985).
  * Xu  _et al._ (2019) Q. Xu, Y. Wang, and Y. Ma, Phys. Rev. B 100, 205132 (2019).
  * Mi  _et al._ (2018) W. Mi, A. Genova, and M. Pavanello, The Journal of Chemical Physics 148, 184107 (2018).
  * Kryachko and Ludeña (2014) E. S. Kryachko and E. V. Ludeña, Physics reports 544, 123 (2014).
  * Kohn (1983) W. Kohn, Phys. Rev. Lett. 51, 1596 (1983).
  * Garg and Colò (2018) U. Garg and G. Colò, Progress in Particle and Nuclear Physics 101, 55 (2018).
  * Ring and Schuck (1980) P. Ring and P. Schuck,  _The nuclear many-body problem_ (Springer-Verlag, New York, 1980).
  * Waechter (2009) A. Waechter, in  _Combinatorial Scientific Computing_ (2009).
  * Wächter and Biegler (2005) A. Wächter and L. T. Biegler, SIAM J. on Optimization 16, 1 (2005).
  * Perdew  _et al._ (1982) J. P. Perdew, R. G. Parr, M. Levy, and J. L. Balduz, Phys. Rev. Lett. 49, 1691 (1982).
  * Perdew and Levy (1997) J. P. Perdew and M. Levy, Phys. Rev. B 56, 16021 (1997).
  * Chabanat  _et al._ (1997) E. Chabanat, P. Bonche, P. Haensel, J. Meyer, and R. Schaeffer, Nuclear Physics A 627, 710 (1997).
  * Colò  _et al._ (2013) G. Colò, L. Cao, N. V. Giai, and L. Capelli, Computer Physics Communications 184, 142 (2013).


# Imitation learning for sim-to-real transfer of robotic cutting policies
based on residual Gaussian process disturbance force model

Jamie Hathaway1,2,†, Rustam Stolkin1,2, Alireza Rastegarpanah1,2,† 1
Department of Metallurgy & Materials Science, University of Birmingham,
Birmingham, UK, B15 2TT.2 The Faraday Institution, Quad One, Harwell Science
and Innovation Campus, Didcot, UK, OX11 0RA.$\dagger$ These authors
contributed equally to this work.This work was supported in part by the
project called “Research and Development of a Highly Automated and Safe
Streamlined Process for Increase Lithium-ion Battery Repurposing and
Recycling” (REBELION) under Grant 101104241 and in part by the UK Research and
Innovation (UKRI) project “Reuse and Recycling of Lithium-Ion Batteries”
(RELIB) under RELiB2 Grant FIRG005 and RELIB3 Grant FIRG057.The authors would
like to acknowledge Abdelaziz Wasfy Shaarawy, Cesar Contreras, Carl Meggs and
Christopher Gell respectively for assistance with experimental validation,
design of material holder and cutter tool for experiments herein.This work has
been submitted to the IEEE for possible publication. Copyright may be
transferred without notice, after which this version may no longer be
accessible.

###### Abstract

Robotic cutting, or milling, plays a significant role in applications such as
disassembly, decommissioning, and demolition. Planning and control of cutting
in real-world scenarios in uncertain environments is a complex task, with the
potential to benefit from simulated training environments. This letter focuses
on sim-to-real transfer for robotic cutting policies, addressing the need for
effective policy transfer from simulation to practical implementation. We
extend our previous domain generalisation approach to learning cutting tasks
based on a mechanistic model-based simulation framework, by proposing a hybrid
approach for sim-to-real transfer based on a milling process force model and
residual Gaussian process (GP) force model, learned from either single or
multiple real-world cutting force examples. We demonstrate successful sim-to-
real transfer of a robotic cutting policy without the need for fine-tuning on
the real robot setup. The proposed approach autonomously adapts to materials
with differing structural and mechanical properties. Furthermore, we
demonstrate the proposed method outperforms fine-tuning or re-training alone.

##  I Introduction

Contact cutting processes feature extensively in a range of applications, most
notably manufacturing. However, disassembly applications face particular
challenges due to limited availability of Computer-Aided Design (CAD) models,
material diversity, and product condition uncertainties. These imply differing
challenges from manufacturing as lack of models and prior knowledge about the
task introduces challenges for traditional path and process planning
approaches. Recently, learning-based approaches have proven successful for
online adaptation to uncertainties in a range of application environments,
albeit mostly non-destructive [1, 2, 3, 4]. For destructive tasks, simulation
environments present a compelling option to avoid the overhead of costly data
collection. For example, previous work [5] successfully demonstrated learning
and adaptation to a range of uncertain milling environments; however,
challenges remain in sim-to-real adaptation.

While high precision dynamometer setups are extensively used in research
applications, challenges are presented for disassembly applications due to
geometric limitations and relatively high cost [6, 7]. Force-torque sensors,
either mounted at the wrist, or directly integrated into the robot, offer
comparatively greater flexibility and lower cost. However, sensor limitations
result in considerable differences between simulation and real world from
factors such as tool vibration, motor cogging torque and residual unmodelled
dynamic effects such as chattering, which are difficult to model without
laborious identification. Learning-based approaches suffer poor generalisation
when dealing with distributional mismatch between source and target domain
observations. Furthermore, direct re-training is problematic due to the
problem of catastrophic forgetting, resulting in large reductions in
performance [8]. Hence, direct sim-to-real adaptation (i.e. without re-
training / direct supervision on the target domain) is desirable due to safety
issues and limited sample availability.

This letter proposes an imitation learning based sim-to-real adaptation
approach for learning milling tasks based on residual modelling of
disturbances directly from offline real world cutting demonstrations. We
demonstrate successful sim-to-real transfer of a cutting policy without
requirement for supervised domain adaptation on the target domain, based on a
series of real world cutting trials on unknown materials. We furthermore
demonstrate the proposed framework outperforms fine-tuning or direct re-
training alone. The approach is prior-knowledge independent in the sense that
only assumption of a periodic disturbance force is made.

The letter is structured as follows: we survey related work in the area of
robotic milling, transfer learning and domain adaptation approaches in Section
II. Section III introduces the proposed method based on residual Gaussian
process modelling of real world cutting examples and knowledge transfer via
imitation learning. Section IV details evaluation of the method for simulated
and real world cutting trials. Finally, Section V concludes the letter.

##  II Related Work

Previous approaches to sim to real adaptation of policies can be broadly
categorised into domain adaptive and transfer learning approaches [9]. In the
context of reinforcement learning (RL), the aim of the former is to obtain a
policy representation that can adapt to the real world via the alignment of
conditional or marginal distributions of observations between source and
target domains. Previously, domain adaptation approaches have been proposed
based on unpaired translation models [10]. In [11], a 3-stage framework for
imitation learning from visual demonstrations was proposed. Another approach
for unsupervised domain adaptation was proposed in [12] based on a deep
spherical manifold Gaussian kernel to align source and target domain
distributions. Zhao et al. [13] introduces a reinforcement learning-based
approach for deformation control in machining across different samples. A
joint invariant feature representation is learned for pairs of samples which
is used to train a meta-model which adapts to new tasks by fine-tuning. A
notable drawback is that for each sample a new policy must be trained to
initially generate the meta-model; however, once the meta model is generated,
the approach can be fine-tuned with minimal data. A related concept for domain
adaptation is proposed in [14] by learning policies within a unified latent
space representation of the source and target domain observations.

Another family of approaches focuses on domain adaptation for classification
based on a combination of classifier and discriminator models. These are
reminiscent of image classification works by leveraging a backbone trained
across multiple tasks to generalise across different domains [15]. In [16],
this is employed with a weighted combination of global and local
discriminators for generalisation of tool wear classification across different
milling tools. However, these approaches are still costive in terms of data
collection. Consider that in addition to labelled source domain data, a large
amount of _unlabelled_ target domain data are necessary to train the
discriminators. To address the problem of incomplete target domain data, [17]
proposed a domain adaptation method for condition monitoring of different
milling tools based on a generative CNN. Besides data collection, domain
adaptive approaches typically make assumptions about the differences between
domains. A common assumption is that the conditional distributions of outputs,
such as classifiers, are domain invariant, while differences between the
domains arise from differences in the marginal distributions over
observations, which is not always the case [18]. In this work, we consider the
case that the conditional distributions differ between domains, which
represents the “conditional shift” case from [18].

Related to the domain adaptation view on sim-to-real transfer, another
approach is to compensate directly for measured disturbances on the real
setup. For example, an online disturbance rejection approach is presented in
[19] based on a combination of a disturbance observer with a GP regression
model of position-based disturbances, while more specifically to milling, [7]
proposes a 4-inertial disturbance observer approach for identification of
milling force models. However, the authors suggest that some level of prior
knowledge may be required to obtain initial estimates for accurate parameter
identification.

Transfer learning approaches contrast domain adaptive approaches in that they
aim to improve performance on a target domain via transfer of knowledge
between domains. This transfer can be unidirectional (e.g. direct fine-tuning
of an existing policy on real world), or bi-directional, such as incorporating
knowledge obtained from real world experiments to inform simulation design.
Recently, sim-to-real transfer approaches have been proposed based on
augmentation of simulations with models learned from trajectories collected
from the real environment. In [20] a neural augmented simulation approach was
proposed based on a deep neural network dynamic model learned from real world
trajectories. Similarly, [21] proposed learning of an inverse dynamics model
as a deep neural network for sim-to-real adaptation. However, these approaches
require a large number of real world samples to train. On the other hand,
Gaussian processes (GPs) have been demonstrated to be efficient at learning
from a small number of samples. This has resulted in sim-to-real and domain
adaptation approaches in other application areas [22]. Another GP-based
approach based on synthesis of data from multiple experts is proposed in [23]
for identification of facial expressions, where the target domains are
different faces. Alternative sim-to-real transfer approaches have employed
system identification methods to optimise physical parameters of a simulation
to maximise target domain performance [24]. Therefore, these methods are
capable of only modelling parametric differences between source and target
domains. Related to this concept and this present work, [25] leverages an
imitation-learning-based approach for domain adaptation based on real world
demonstrations on a misspecified simulator.

##  III Methodology

Figure 1: Overview of the proposed framework.

###  III-A Dataset Preparation

We consider an initially unstructured dataset comprising time series
$\mathcal{D}_{u}=\\{\mathbf{X}_{0}\ldots,\mathbf{X}_{n}\\}$. Each times series
contains force measurements $\mathbf{X}_{i}=\\{\boldsymbol{f_{e}}\\}$. Many
periodic disturbances are parameterised as a function of position instead of
time [19]. In the absence of position measurements, the data must first be
aligned in the time domain. Each time series is first normalised as:

$\mathbf{X}_{i}=\Sigma_{\mathbf{X}_{i}}^{-1}\left(\mathbf{X}_{i}-\bar{\mathbf{X}}_{i}\right)$
(1)

where $\bar{\mathbf{X}}_{i}$, $\Sigma_{\mathbf{X}_{i}}$ are the mean and
standard deviations of each example time series $\mathbf{X}_{i}$. The
alignment in time domain consists of an initial coarse alignment stage,
followed by a fine alignment stage. The coarse alignment is performed by
maximising the cross-correlation:

$t_{\mathrm{delay}}=\underset{t}{\mathrm{argmax}}\left((\mathbf{X}_{i}*\mathbf{X}_{j}))(t)\right)$
(2)

Subsequently, the fine alignment is performed with pairwise dynamic time
warping (DTW) between each example trajectory pair $\mathbf{X}_{i}$,
$\mathbf{X}_{j}$. DTW finds an optimal warping path:

$(l^{*}_{n},m^{*}_{n})=\underset{(l_{n},m_{n})}{\mathrm{argmin}}\sum_{n}\frac{d(\boldsymbol{x}_{l_{n}},\boldsymbol{x}_{m_{n}})w_{n}}{\sum_{n^{\prime}}{w_{n^{\prime}}}}\quad$
(3)

with respect to a cumulative weighted distance between features
$d(\boldsymbol{x}_{l},\boldsymbol{x}_{m})$ (here Euclidean distance), where
$w_{n}$ are weighting terms. Adopting the convention from [26], the
“symmetric2” step pattern was used. To account for variable length examples,
we employ an open-ended DTW approach, which allows for more flexible alignment
by permitting variable-length warping paths, which further computes the
minimum (3) over all truncated sequences of $\mathbf{X}_{j}$. To re-index a
time series according to the warping path, each element of the original time
series is assigned to a new time point based on the optimal mapping. This re-
indexed time series $\hat{\mathbf{X}}$ effectively accounts for time
distortions, allowing synchronization and alignment of the data across
different examples. Hence, the final dataset comprises aligned time series
$\mathcal{D}=\\{\hat{\mathbf{X}}_{0}\ldots,\hat{\mathbf{X}}_{n}\\}$.

###  III-B GP Force Model

We consider the regression problem:

$\boldsymbol{f}_{e}=\boldsymbol{f}(t)+\boldsymbol{d}(t)+\epsilon$ (4)

where $f(t)$ represents the force computed from the underlying process force
model, while $\boldsymbol{d}(t)$ represents a disturbance force. The
disturbance $\boldsymbol{d}(t)$ is assumed to be periodic and with independent
and identically distributed (i.i.d) noise
$\epsilon\sim\mathcal{N}(0,\sigma^{2})$. We model the disturbance force
$\boldsymbol{d}(t)$ as a Gaussian Process (GP) with zero mean and covariance
function $k(\boldsymbol{x},\boldsymbol{x}^{\prime})$

$\boldsymbol{d}(\boldsymbol{t})\sim\mathcal{GP}(\boldsymbol{0},k(\boldsymbol{t},\boldsymbol{t}^{\prime}))$
(5)

where the distribution of observed and unobserved data is modeled as a joint
multivariate Gaussian distribution:

$\begin{bmatrix}\boldsymbol{d}\\\
\boldsymbol{d}_{*}\end{bmatrix}=\mathcal{N}\left(\mathbf{0},\begin{bmatrix}\mathbf{K}&\mathbf{K}_{*}\\\
\mathbf{K}_{*}^{\mathsf{T}}&\mathbf{K}_{**}\end{bmatrix}\right)$ (6)

where $\mathbf{K}$, $\mathbf{K}_{*}$, $\mathbf{K}_{**}$ are the training,
train-test and test covariance matrices respectively. The posterior predictive
distribution for the disturbance force $\boldsymbol{d}_{*}$ given test points
$\mathbf{X}_{*}$ is then

$\displaystyle
p(\boldsymbol{d}_{*}|\mathcal{D},\boldsymbol{d},\mathbf{X}_{*})=$
$\displaystyle\mathcal{N}\left(\boldsymbol{\mu},\mathbf{\Sigma}\right)$ (7)
$\displaystyle\boldsymbol{\mu}=$
$\displaystyle\mathbf{K}_{*}^{\mathsf{T}}\left[\mathbf{K}+\sigma^{2}\mathbf{I}\right]^{-1}\boldsymbol{d}$
(8) $\displaystyle\mathbf{\Sigma}=$
$\displaystyle\mathbf{K}_{**}-\mathbf{K}_{*}^{\mathsf{T}}\left[\mathbf{K}+\sigma^{2}\mathbf{I}\right]^{-1}\mathbf{K}_{*}$
(9)

To capture the periodic nature of the disturbance force, we use the
exponential sine covariance function

$k(\boldsymbol{x},\boldsymbol{x}^{\prime})=\exp\left(\frac{2\sin^{2}\left(\frac{\pi
d(\boldsymbol{t},\boldsymbol{t}^{\prime})}{p}\right)}{l^{2}}\right)$ (10)

with periodicity $p$, length scale $l$ and Euclidean distance
$d(\boldsymbol{t},\boldsymbol{t}^{\prime})$.

Based on the assumptions of the mechanistic force modelling approach [27], the
cutting force $\boldsymbol{f}$ can be related to each cutting edge
$p\in{1\ldots N_{p}}$ of a fluted cutting tool via mechanistic constants
$\boldsymbol{k}_{c}$, $\boldsymbol{k}_{e}$:

$\boldsymbol{f}^{\mathcal{P}}=b_{p}\boldsymbol{k}_{e}+b_{p}\boldsymbol{k}_{c}h_{p}$
(11)

where $b_{p}$ is the thickness of the cutting edge, $h_{p}$ is the uncut chip
thickness, computed from the cutting edge rotation angle $\theta_{p}$ in the
tool model reference frame, tool feed rate $v$ and spindle speed $\omega$:

$h_{p}=\sin\theta_{p}\frac{v}{N_{p}\omega}$ (12)

The total cutting force is computed as the sum over cutting elements weighted
by the Boolean vector $\mathbf{G}\in\mathbb{B}^{N_{p}}$ specifying the
engagement of each element $k$ with the workpiece.

$\boldsymbol{f}=\sum_{p}^{N_{p}}G_{p}\mathbf{R}_{\mathcal{P}}^{\mathcal{}}\boldsymbol{f}^{\mathcal{P}}$
(13)

where $\mathbf{R}_{\mathcal{P}}^{\mathcal{}}$ is the rotation of cutting edge
$p$ about the tool axis to the tool model frame.

###  III-C Imitation Learning Framework

The cutting task is formulated in the Mujoco simulation environment, into
which the model of cutting mechanics from (13) is embedded. The expert reward
function $r$ expresses a weighted sum of task-specific performance metrics
$\mathrm{MRV}$ (material removed volume) and cutting time $t_{\mathrm{cut}}$,
and feasibility reward shaping comprising path error $\boldsymbol{e}$, process
force $\boldsymbol{f}$

$r=Q_{\textrm{MRV}}\cdot\mathrm{MRV}-Q_{\mathrm{cut}}t_{\mathrm{cut}}-\boldsymbol{e}\mathbf{Q}_{d}\boldsymbol{e}^{\mathsf{T}}-\boldsymbol{f}\mathbf{Q}_{f}\boldsymbol{f}^{\mathsf{T}}$
(14)

with observations
$\boldsymbol{o}=\begin{bmatrix}\dot{\boldsymbol{c}}(t)\dot{\boldsymbol{v}}_{EE}^{\mathsf{T}}&\boldsymbol{e}&\boldsymbol{v}_{EE}&\boldsymbol{f}_{e}&\boldsymbol{a}\end{bmatrix}^{\mathsf{T}}$,
where $\boldsymbol{c}(t)$ is the reference cutting path, end-effector velocity
$\boldsymbol{v}_{EE}$, and agent actions $\boldsymbol{a}$. Due to the problems
of catastrophic forgetting with fine-tuning and high cost of data collection
in the real environment, we propose an imitation-learning based approach to
train a target policy that can adapt to the target domain via a surrogate
target domain. We establish a test case comprising “offline” and “online”
imitation learning algorithms; behavioural cloning (BC) and DAgger, which we
contrast with the case of the expert policy directly transferred with no fine-
tuning, and the expert with fine-tuning directly on the surrogate target
domain.

During training, we employ the source domain expert policy $\pi_{e}$ to train
a surrogate target domain policy $\hat{\pi}_{t}$. To resemble the fine-tuning
case, we initialise the target domain policy as $\hat{\pi}_{t}=\pi_{e}$. At
each step, the tuple $(\boldsymbol{o}_{e},\pi_{e}(\boldsymbol{o}_{e}))$ is
sampled from the base environment using the data collection policy

$\pi_{d}=\beta\pi_{e}+(1-\beta)\pi_{t}$ (15)

where $\beta$ is the non-expert action probability, which is non-zero for
DAgger and zero otherwise. Then, the expert observations $\boldsymbol{o}_{e}$,
actions $\boldsymbol{a}$ $(\boldsymbol{o}_{e},\boldsymbol{a})$ are modified by
sampling from the posterior distribution (7) for each point in the trajectory
to generate new experiences $(\boldsymbol{o}_{t},\boldsymbol{a})$.
Subsequently, the surrogate target domain experiences are used to update the
target policy as the standard behavioural cloning procedure. However, the
method is applicable to other imitation learning algorithms such as GAIL or
AIRL. This procedure is summarised in Algorithm 1. We employ the proximal
policy optimisation (PPO) algorithm for learning and fine-tuning of policies,
although the principles are independent of learning algorithm. For each
approach, we conduct training with all strategies for 50 episodes, a learning
rate of $1\times 10^{-3}$ and batch size of 64. For DAgger, $\beta$ is varied
according to a 0–1 linear schedule for 45 episodes.

Algorithm 1 Imitation-learning sim-to-real transfer

Expert policy, $\pi_{e}$, target policy $\hat{\pi}_{t}$

Source domain (expert) environment $\mathcal{E}$

Disturbance model
$p(\boldsymbol{d}^{\prime}|\mathcal{D},\boldsymbol{x},\boldsymbol{d},\boldsymbol{x}^{\prime})$

$\hat{\pi}_{t}\leftarrow\pi_{e}$

for $i=0$ to $N$ do

while episode not done do

Sample $(\boldsymbol{o}_{e},\pi_{e}(\boldsymbol{o}_{e}))$ from $\mathcal{E}$

Sample $\boldsymbol{d}^{\prime}\sim
p(\boldsymbol{d}^{\prime}|\mathcal{D},\boldsymbol{x},\boldsymbol{d},\boldsymbol{x}^{\prime})$

$\boldsymbol{o}_{t}\leftarrow\boldsymbol{o}_{e}+\boldsymbol{d}^{\prime}$

Append $\mathcal{D}_{e}$ with
$(\boldsymbol{o}_{t},\pi_{e}(\boldsymbol{o}_{e}))$

Update $\mathcal{E}$ with $a\sim\hat{\pi}_{t}(\boldsymbol{o}_{e})$

end while

end for

Train $\hat{\pi}_{t}$ with $\mathcal{D}_{e}$ as
$\mathsf{BC},\mathsf{DAgger},\ldots$

###  III-D Experimental Setup

The real world setup consists of a KUKA LBR iiwa R820 14kg collaborative robot
equipped with a wrist mounted motorised slitting saw tool. Although the iiwa
has built-in torque sensing capabilities, we consider the case where a force-
torque sensor (FT-AXIA 80) is mounted between the cutter tool and the robot
flange which measures the process force. The policy controls the robot via an
operational space computed torque tracking controller, with control law:

$\boldsymbol{\tau}=\mathbf{J}^{\mathsf{T}}\left[\mathbf{\Lambda}(\boldsymbol{q})\left[\mathbf{K}_{d}(t)\dot{\boldsymbol{e}}+\mathbf{K}_{p}(t)\boldsymbol{e}\right]+\mathbf{\Gamma}+\boldsymbol{\mu}\right]$
(16)

where $\boldsymbol{\tau}$ are the commanded joint torques, $\mathbf{J}$ the
robot Jacobian, and $\mathbf{\Lambda}$, $\boldsymbol{\Gamma}$,
$\boldsymbol{\mu}$ the robot dynamic parameters corresponding to inertia,
Coriolis / centrifugal forces and gravitational forces respectively. The
policy outputs correspond to the controller translational stiffness
$\mathbf{K}_{p}$ and translational setpoint adjustment, consisting of a feed
rate modification from a nominal feed rate and depth of cut (DoC) offset from
the reference trajectory. The damping gain $\mathbf{K}_{d}$ is computed
according to give a critically damped behaviour for the selected stiffness.
For each strategy, we consider a single cutting pass of a planar material by
conventional milling. The reference tool path is defined with respect to the
material surface position, with a reference DoC of $0$mm – hence, the actual
DoC is directly selected via the policy DoC offset.

Figure 2: Overview of the experimental setup used for real world cutting
experiments. Figure 3: Plot of measured disturbance force in feed (Y)
direction from dataset of examples after temporal alignment. The Gaussian
process model fit is shown; the shaded areas show 1-$\sigma$ and 2-$\sigma$
standard deviations from the mean respectively. The dashed line shows the
transition from training data (left) to extrapolation (right).

##  IV Results & Discussion

In this section, the performance of the proposed method is evaluated in
simulation in the context of the performance of the expert in the source
domain (i.e. simulation without GP augmentation) and surrogate target domain
(simulation augmented with GP). The performance and behaviour is subsequently
compared for the true target domain for a series of real world cutting tasks.

###  IV-A Simulation studies

(a) Actions

(b) States

Figure 4: Comparison of policy actions and states between source domain expert
in source and surrogate target domains (Expert, Expert GP), and surrogate
target domain policies using fine-tuning, behavioural cloning (BC) and DAgger
imitation learning strategies. States include the path error transverse
$e_{x}$ and normal $e_{z}$ to the planned path and forces in the feed $F_{y}$
and normal direction $F_{z}$. Forces are shown with a 50-point (1s) moving
average filter.

We first establish a case study for cutting of a planar material from a set of
50 trials with randomly chosen mechanistic constants for the simulation
augmented with the learned GP model. We compare the performance of the
“expert” agent, which is trained in the unaugmented simulation for 32000
episodes and examine the behaviour with respect to the policies trained with
the proposed GP + imitation strategy, with the performance in the base case,
i.e. unaugmented simulation as a benchmark. Figure 4 shows the actions adopted
by each strategy, consisting of the relative feed rate versus the nominal
(0.75m/min), DoC and controller stiffness $\boldsymbol{K}_{p}$.
Correspondingly, the path deviations in the transverse ($e_{x}$) and normal
($e_{z}$) directions are shown, along with forces in the feed and normal
directions ($F_{y}$, $F_{z}$ respectively).

From Figure 4(a) the behaviour of the expert in the source domain can be seen.
When transferred directly into the surrogate target domain (“Expert GP”), the
policy reacts aggressively to the added disturbance, exhibiting sporadic
variations in feed rate and stiffness selection. In comparison, the surrogate
target domain policies trained with both imitation learning strategies closely
tracks the source domain expert behaviour. Conversely, while the actions
adopted by the original expert policy after fine-tuning differ, particularly
for the stiffness Z component, however, the fine-tuned policy remains
similarly sensitive to the disturbances. Comparing the states in Figure 4(b)
further corroborates the similar performance achieved by both imitation
learned policies, with DAgger tracking the expert behaviour more closely than
BC. The fine-tuned policy adopts a pattern to the expert as directly applied
to the surrogate target domain, with a delay of $\sim$0.1s. As an aside, note
that with low-pass filtering of the forces, the underlying trend, as seen in
the “expert” evaluation, is not recovered, even with aggressive cutoff. This
also has the effect of introducing delay into the measured signal; for
example, a 10th-order Butterworth filter with a cutoff defined at 5Hz has a
maximal group delay of $\sim$0.5s, or 25 policy evaluations.

Broadly, for all strategies, similar patterns in actions can be observed
during each phase of the cutting task. For example, the policy adopts a high
stiffness and feed rate broadly around the nominal during approach, followed
by transitioning to a high feed rate and low stiffness in the normal direction
(Z) after impact. Furthermore, the DoC remains consistent across all trials.
Therefore, the comparison of each strategy implies predominantly _behavioural_
benefits of the proposed approach, rather than fundamentally altering the
“decision-making” of the original expert strategy.

Figure 5: Comparison of training curves between base environment and
environment with GP residual force model. Figure 6: Violin plot of reward
distribution between source domain expert policy and target domain strategies
– fixed process parameters for all materials (baseline), unmodified and fine-
tuned expert policies, and target policies trained with BC and DAgger
imitiation learning approaches.

Figure 5 shows the training curves for expert policies trained from scratch in
source and surrogate target domains respectively. The source policy converges
rapidly in the first 3M training steps, before a phase of gradual improvement
over the remaining 13M steps to a final reward of -0.766. For the surrogate
target policy, an initial reduction in performance is recorded, followed by
gradual improvement from steps 5.5–12M until convergence to a final reward of
-2.24. Besides the advantage of training time reduction (50 episodes vs.
32000, $\sim$0.16% of training time), the final rewards from the surrogate
target expert policy are notably reduced. This follows from the inclusion of
the GP model representing a more challenging task for the agent; the force
observations in particular are directly related to the reward and encode
important information about the interaction which relates to the selection of
milling process parameters.

We consider the performance of the original expert policy, the target domain
fine tuned policy and policies with each imitation learning strategy (BC,
DAgger) in both the original simulation environment and simulation augmented
with the learned GP model from real-world trials. We additionally compare the
performance with fixed process parameters (Baseline) at the nominal (feed rate
0.75m/min and 1mm DoC) for all materials. In each instance, we compare the
episodic rewards over 50 simulation rollouts. From Figure 6 it is clear that
the overall performance of the expert and all strategies is robust to
distributional mismatch between source and target domains. For the fine-tuned,
BC and DAgger strategies, the performance is slightly more consistent in the
extreme case as indicated by the higher minimum rewards obtained. The greatest
deviation is observed for the DAgger trained target policy from the expert for
both the original and GP-augmented environments. As DAgger allows for the
trajectories to deviate from those obtained using the expert demonstrations
alone, it enables the policy to explore and adapt further to the surrogate
target domain. It is therefore unsurprising that DAgger demonstrates superior
performance in this case.

###  IV-B Real world cutting trials

(a) Polyurethane foam

(b) Mica sheet

Figure 7: Comparison of policy actions between source domain expert and
surrogate target domain policies using behavioural cloning (BC) and DAgger
imitation learning strategies. Actions include the relative feed rate
adjustment vs. nominal (0.75m/min), depth of cut (DoC) and controller
stiffness $\mathbf{K}_{p}$.

(a) Polyurethane foam

(b) Mica sheet

Figure 8: Comparison of path error transverse $e_{x}$ and normal $e_{z}$ to
the planned path and forces in the feed $F_{y}$ and normal direction $F_{z}$
between source domain expert, BC and DAgger surrogate target domain policies.
Forces are shown with a 50-point (1s) moving average filter.

We evaluate the cutting policy for conventional milling of 5 different
materials – mica, cardboard, corrugated plastic, high-density polyurethane
(PU) foam and aluminium – which exhibit broadly differing mechanical and
structural properties. Figure 7 and Figure 8 show the actions and key states
for the expert, BC and DAgger trained policies for foam and mica. Emblematic
of the source domain expert as applied directly to the real world is a high
degree of variability in the agent actions over time, particularly of the
controller stiffness $\mathbf{K}_{p}$, corroborating the behaviour observed in
simulation. Selection of feed rate appears sporadic and largely uncorrelated
to any particular evolution of the cutting process; for example, an increase
in loading of the tool. This can be contrasted with the feed rate selection of
the DAgger trained policy, where reductions in feed rate are associated with
increases in the DoC due to increased loading on the tool. The behaviour of
the expert policy is inconsistent, with path error and force progressively
increasing during the cutting task, whereas for the BC and DAggger policies,
these are relatively consistent. However, it is important to note that since
the focus of this work is on “rough” cutting for disassembly or
decommissioning applications, and not on achieving high dimensional
tolerances, the expected errors are higher than would be expected for milling,
e.g. for a manufacturing application.

Figure 9: Breakdown of average reward and reward components by material and
strategy. The performance of each strategy, labelled “Ex”, “BC”, “DA” and “Ba”
respectively are grouped by the material; here abbreviated Mica, Cardboard,
Plastic, Foam and Aluminium respectively. The total reward components for each
strategy (i.e. total reward) is marked.

Some behavioural characteristics are not preserved from the simulation case
study; for example, the selected feed rate during the approach stage of
cutting remains nearly at the maximum for all strategies, contrasting the with
simulation case. Similarly, the stiffness along the normal direction (Z) is
increased during the approach phase, however to a lesser extent than in
simulation. However, during the cutting phase, these are broadly more similar
to the simulation case. As previously established, the force observations
encode key information about the interaction, including the discrimination of
contact and non-contact states.

A comparison of the average rewards over all materials for the expert policy
and BC and DAgger trained target policies is shown in Figure 9. Across all
trials, consistently lower forces are encountered with the policy-based
approach, with higher total rewards observed for all materials with the
exception of foam. Similarly, this trend is observed between the expert and
imitation-trained policies, while (excluding foam) retaining similar or
improved MRV across the range of materials, implying an improved capability to
regulate the cutting force for the real world case studies. For overall
rewards, a significant difference in between the BC and expert was found for
mica (unpaired two-tailed T-test, $p=0.015$; 95% significance level), with the
observed effect being improvement, but not for DAgger ($p=0.33$); the higher
MRV and force implying a more aggressive behaviour. In cardboard, significant
differences were found in both cases (BC $p=0.021$, DAgger $p=0.020$), with
the observed effect being improvement in both cases. For plastic, mixed
results are seen (BC $p=0.14$, DAgger $p=0.038$) consisting of inconclusive
improvement for BC, but observed improvement in the DAgger case. Similarly for
foam, BC performed worse than DAgger, differing markedly from the expert with
a reduction in reward, ($p=0.00014$), with an inconclusive difference in
performance for DAgger ($p=0.084$). Finally for aluminium, both strategies
showed a self-consistent improvement over the expert, but individually were
inconclusive ($p=0.094$, $p=0.090$). Overall, it can generally be said that
DAgger performed better on the softer materials, whereas BC performed well on
stronger materials such as mica. For DAgger, the performance on the range of
case studies was either improved, or – in the least case – did not
significantly differ from the expert. This can similarly be said for BC with
the exception of foam.

##  V Conclusion

An imitation-learning based approach for sim-to-real transfer of a robotic
cutting policy was proposed. We demonstrate how the residual process dynamics
and disturbances can be modelled from a small number of real world trials. We
validate the proposed method on a real robot setup, demonstrating the policies
transferred to the real world in many cases to have significantly improved
performance over the expert as directly transferred from simulation.
Furthermore, the behavioural characteristics of the target domain policies
(sensitivity to disturbances) were improved. We also demonstrate the proposed
method based on imitation learning outperforms fine-tuning or re-training
alone, and that improvements can be obtained even with relatively simple
offline imitation learning approaches.

A notable limitation of this work is that while the proposed method can
incorporate data from multiple examples, the disturbances are assumed to be
sampled from the same underlying process. If the disturbances differ between
demonstrations, the alignment between demonstrations will be poor. In this
case, a separate GP model would be required for each variation. The assumption
of a periodic disturbance force furthermore remains a notable limitation. To
address these limitations, future work will explore direct synthesis of
surrogate real-world data from a minimal unstructured dataset of offline
demonstrations.

## References

  * [1] C. C. Beltran-Hernandez, D. Petit, I. G. Ramirez-Alpizar, and K. Harada, “Variable compliance control for robotic peg-in-hole assembly: A deep-reinforcement-learning approach,” _Applied Sciences_ , vol. 10, no. 19, 2020.
  * [2] Y. Wang, C. C. Beltran-Hernandez, W. Wan, and K. Harada, “Hybrid trajectory and force learning of complex assembly tasks: A combined learning framework,” _IEEE Access_ , vol. 9, pp. 60 175–60 186, 2021.
  * [3] X. Li, J. Xiao, W. Zhao, H. Liu, and G. Wang, “Multiple peg-in-hole compliant assembly based on a learning-accelerated deep deterministic policy gradient strategy,” _Industrial Robot: the international journal of robotics research and application_ , vol. 49, no. 1, pp. 54–64, Jan 2022.
  * [4] X. B. Peng, M. Andrychowicz, W. Zaremba, and P. Abbeel, “Sim-to-real transfer of robotic control with dynamics randomization,” in _2018 IEEE International Conference on Robotics and Automation (ICRA)_ , 2018, pp. 3803–3810.
  * [5] J. Hathaway, A. Rastegarpanah, and R. Stolkin, “Learning robotic milling strategies based on passive variable operational space interaction control,” _IEEE Transactions on Automation Science and Engineering_ , pp. 1–14, 2023\.
  * [6] B. Gibson, C. Cox, M. Aguilar, A. Strauss, and G. Cook, “Low-cost wireless force sensor design with applications in friction stir welding,” 06 2012, pp. 70–75.
  * [7] K. Takahei, N. Suzuki, and E. Shamoto, “Identification of the model parameter for milling process simulation with sensor-integrated disturbance observer,” _Precision Engineering_ , vol. 78, pp. 146–162, 2022.
  * [8] S.-I. Ao and H. Fayek, “Continual deep learning for time series modeling,” _Sensors_ , vol. 23, no. 16, 2023.
  * [9] E. Salvato, G. Fenu, E. Medvet, and F. A. Pellegrino, “Crossing the reality gap: A survey on sim-to-real transferability of robot controllers in reinforcement learning,” _IEEE Access_ , vol. 9, pp. 153 171–153 187, 2021\.
  * [10] P. M. Scheikl, E. Tagliabue, B. Gyenes, M. Wagner, D. Dall’Alba, P. Fiorini, and F. Mathis-Ullrich, “Sim-to-real transfer for visual reinforcement learning of deformable object manipulation for robot-assisted surgery,” _IEEE Robotics and Automation Letters_ , vol. 8, no. 2, pp. 560–567, 2023\.
  * [11] D. Zhang, W. Fan, J. Lloyd, C. Yang, and N. F. Lepora, “One-shot domain-adaptive imitation learning via progressive learning applied to robotic pouring,” _IEEE Transactions on Automation Science and Engineering_ , pp. 1–14, 2022.
  * [12] Y. Zhang and B. D. Davison, “Deep spherical manifold gaussian kernel for unsupervised domain adaptation,” in _2021 IEEE/CVF Conference on Computer Vision and Pattern Recognition Workshops (CVPRW)_. Los Alamitos, CA, USA: IEEE Computer Society, Jun 2021, pp. 4438–4447.
  * [13] Y. Zhao, C. Liu, Z. Zhiwei, K. Tang, and D. He, “Reinforcement learning method for machining deformation control based on meta-invariant feature space,” _Visual computing for industry, biomedicine, and art_ , vol. 5, p. 27, 11 2022\.
  * [14] J. Xing, T. Nagata, K. Chen, X. Zou, E. Neftci, and J. L. Krichmar, “Domain adaptation in reinforcement learning via latent unified state representation,” _CoRR_ , vol. abs/2102.05714, 2021.
  * [15] M. Ragab, E. Eldele, W. L. Tan, C.-S. Foo, Z. Chen, M. Wu, C.-K. Kwoh, and X. Li, “Adatime: A benchmarking suite for domain adaptation on time series data,” _ACM Trans. Knowl. Discov. Data_ , vol. 17, no. 8, May 2023.
  * [16] K. Li, M. Chen, Y. Lin, Z. Li, X. Jia, and B. Li, “A novel adversarial domain adaptation transfer learning method for tool wear state prediction,” _Knowledge-Based Systems_ , vol. 254, p. 109537, 2022.
  * [17] C.-B. Chou and C.-H. Lee, “Generative neural network-based online domain adaptation (GNN-ODA) approach for incomplete target domain data,” _IEEE Transactions on Instrumentation and Measurement_ , vol. 72, pp. 1–10, 2023.
  * [18] K. Zhang, B. Scholkopf, K. Muandet, and Z. Wang, “Domain adaptation under target and conditional shift,” in _International Conference on Machine Learning_ , 2013.
  * [19] H. Jung and S. Oh, “Gaussian process and disturbance observer based control for disturbance rejection,” in _2022 IEEE 17th International Conference on Advanced Motion Control (AMC)_ , 2022, pp. 94–99.
  * [20] F. Golemo, A. A. Taiga, A. Courville, and P.-Y. Oudeyer, “Sim-to-real transfer with neural-augmented robot simulation,” in _Proceedings of The 2nd Conference on Robot Learning_ , ser. Proceedings of Machine Learning Research, A. Billard, A. Dragan, J. Peters, and J. Morimoto, Eds., vol. 87. PMLR, 29–31 Oct 2018, pp. 817–828.
  * [21] P. F. Christiano, Z. Shah, I. Mordatch, J. Schneider, T. Blackwell, J. Tobin, P. Abbeel, and W. Zaremba, “Transfer from simulation to real world through learning deep inverse dynamics model,” _CoRR_ , vol. abs/1610.03518, 2016\.
  * [22] K. Wang, J. Ma, K. L. Man, K. Huang, and X. Huang, “Sim-to-real transfer with domain randomization for maximum power point estimation of photovoltaic systems,” in _2021 IEEE International Conference on Environment and Electrical Engineering and 2021 IEEE Industrial and Commercial Power Systems Europe (EEEIC / I &CPS Europe)_, 2021, pp. 1–4.
  * [23] S. Eleftheriadis, O. Rudovic, M. P. Deisenroth, and M. Pantic, “Gaussian process domain experts for model adaptation in facial behavior analysis,” in _2016 IEEE Conference on Computer Vision and Pattern Recognition Workshops (CVPRW)_. Los Alamitos, CA, USA: IEEE Computer Society, Jul 2016, pp. 1469–1477.
  * [24] M. Kaspar, J. D. Muñoz Osorio, and J. Bock, “Sim2real transfer for reinforcement learning without dynamics randomization,” in _2020 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_ , 2020, pp. 4383–4388.
  * [25] S. Jiang, J.-C. Pang, and Y. Yu, “Offline imitation learning with a misspecified simulator,” in _Proceedings of the 34th International Conference on Neural Information Processing Systems_ , ser. NIPS’20, H. Larochelle, M. Ranzato, R. Hadsell, M. Balcan, and H. Lin, Eds., vol. 33, Red Hook, NY, USA, 2020, pp. 8510–8520.
  * [26] T. Giorgino, “Computing and visualizing dynamic time warping alignments in R: the dtw package,” _Journal of statistical Software_ , vol. 31, pp. 1–24, 2009.
  * [27] E. Armarego and R. Brown, _The Machining of Metals_. Prentice-Hall, 1969.