# Remarks on existence and uniqueness of the solution for stochastic partial
differential equations

Benny Avelin  Uppsala University, Department of Mathematics, Sweden
<EMAIL_ADDRESS>and  Lauri Viitasaari  Aalto University School of
Business, Department of Information and Service Management, Finland
<EMAIL_ADDRESS>
###### Abstract.

In this article we consider existence and uniqueness of the solutions to a
large class of stochastic partial differential of form
$\partial_{t}u=L_{x}u+b(t,u)+\sigma(t,u)\dot{W}$, driven by a Gaussian noise
$\dot{W}$, white in time and spatial correlations given by a generic
covariance $\gamma$. We provide natural conditions under which classical
Picard iteration procedure provides a unique solution. We illustrate the
applicability of our general result by providing several interesting
particular choices for the operator $L_{x}$ under which our existence and
uniqueness result hold. In particular, we show that Dalang condition given in
[5] is sufficient in the case of many parabolic and hypoelliptic operators
$L_{x}$.

Mathematics Subject Classifications (2010): 60H15, 60G15, 35C15, 35K58, 35S10.

Keywords: Stochastic partial differential equations, existence and uniqueness,
mild solution, semilinear parabolic equations, hypoelliptic equations.

##  1\. Introduction

In this article we consider the stochastic partial differential equation
(SPDE) of form

$\frac{\partial u}{\partial
t}(t,x)=(L_{x}u)(t,x)+b(t,u(t,x))+\sigma(t,u(t,x))\dot{W}(t,x),\hskip
14.22636ptt\geq 0,x\in\mathbb{R}^{d}$ (1.1)

with initial condition $u(0,x)=u_{0}(x)$. Here $b$ and $\sigma$ are assumed to
be Lipschitz continuous functions and bounded on compacts in the spatial
variable $x$, and uniformly in $t$. The stochastic force $\dot{W}$ is assumed
to be Gaussian with correlation structure that is white in the time variable
$t$ and given by a generic covariance $\gamma$ in the spatial variable $x$.

SPDEs have been a subject of active research in the literature for recent
years, and the basic theory is already rather well-established. Especially,
initiated by the seminal paper by Dalang [5], stochastic heat equations, in
which case $L_{x}=\Delta$ is the Laplace operator, with Gaussian noise that is
white in time have received a lot of attention. In this case one obtains the
existence and uniqueness of the solution provided that the following so-called
Dalang’s condition

$\int_{\mathbb{R}^{d}}\frac{1}{\beta+2|\xi|^{2}}\hat{\gamma}(d\xi)<\infty$
(1.2)

holds for some $\beta_{0}>0$ (in which case it holds for all
$\beta>\beta_{0}$). Here $\hat{\gamma}$ denotes the non-negative measure
arising as the Fourier transform of the spatial covariance $\gamma$ of the
noise $\dot{W}$. One also observes that here the term $|\xi|^{2}$ corresponds
to the Fourier multiplier arising from the Laplace operator $\Delta$. This
leads to a natural extension to the case where $L_{x}$ is the
$L^{2}$-generator of a Levy process. In this case 1.2 is replaced by

$\int_{\mathbb{R}^{d}}\frac{1}{\beta+2\text{Re}\Psi(\xi)}\hat{\gamma}(d\xi)<\infty,$
(1.3)

where $\Psi(\xi)$ arises from the characteristic exponent of the associated
Levy process. Here $\text{Re}\Psi(\xi)$ is again a non-negative function, and
if 1.3 holds for some $\beta>0$, then it holds for all $\beta>0$. As a
particular interesting example, this covers the case of the stochastic
fractional heat equation where $L_{x}$ is given by the fractional Laplace
operator $-(-\Delta)^{\alpha}$, $\alpha\in(0,1]$. Indeed, then
$\Psi(\xi)=|\xi|^{2\alpha}$ and one recovers the classical heat equation and
condition 1.2 by plugging in $\alpha=1$. This generalisation is studied, among
others, in [10] where the existence and uniqueness result was given under
suitable assumptions on the coefficients $b,\sigma$.

A standard technique to prove existence and uniqueness of the solution to 1.1
is based on the so-called fundamental solution $G_{t}(x)$ (or the Green
kernel) associated to the equation $\partial_{t}u=L_{x}u$. Then one obtains a
candidate for the solution through convolutions with the kernel $G_{t}(x)$. In
particular, in the settings mentioned above one obtains that $G_{t}$ is the
density of the underlying Levy process and forms a semigroup. After that,
using the positivity of $G_{t}$ as well, one obtains the solution through
Picard iteration.

In this article we provide a general existence and uniqueness result for the
solution to 1.1. Our condition, given by Equation 2.6 below, is similar to
1.3. However, in our results we do not require non-negativity of the
fundamental solutions $G_{t}$ as we only consider upper bounds for $|G_{t}|$.
Furthermore, we do not require the associated function, given by $\Psi(\xi)$
in 1.3, to be non-negative. This fact allows to consider a larger class of
operators $L_{x}$ in 1.1, in the case where we do not have a Gaussian upper
bound on the fundamental solution. Our main contribution is that, by a careful
analysis on the essential requirements in the classical arguments, we get a
remarkably strong generalization. As such, we are able to cover a host of
equations which was previously not known in the literature, including the
important Kolmogorov equations. Finally, we stress that similar considerations
can be applied to equations of form 1.1, where $\partial_{t}$ is replaced with
more general operator $L_{t}$, see Section 4.

The rest of the paper is organised as follows. In Section 2 we present our
assumptions and provide the existence and uniqueness result, Theorem 2.2. In
Section 2.1, we illustrate the applicability of our result by providing a
detailed discussion on several interesting examples. The proof of Theorem 2.2
is postponed to Section 3 where we also recall some basic facts on stochastic
calculus. We end the paper with some concluding remarks.

##  2\. General existence and uniqueness result

We consider the stochastic partial differential equation

$\frac{\partial u}{\partial
t}(t,x)=(L_{x}u)(t,x)+b(t,u(t,x))+\sigma(t,u(t,x))\dot{W}(t,x),\hskip
14.22636ptt\geq 0,x\in\mathbb{R}^{d}$ (2.1)

where $L_{x}$ is a suitable differential operator (acting on the variable
$x$), and coefficients $b$ and $\sigma$ are assumed to be spatially Lipschitz
continuous functions, uniformly in $s$ over compacts. That is, for any $T>0$
we have for all $t\in[0,T]$ that

$|b(t,x)-b(t,y)|\leq L_{b,T}|x-y|$

and

$|\sigma(t,x)-\sigma(t,y)|\leq L_{\sigma,T}|x-y|.$

Note that then the functions $b$ and $\sigma$ are spatially locally bounded,
uniformly in $s$ over compacts. That is,

$\sup_{s\in[0,T],y\in K}\max(|\sigma(s,y)|,b(s,y)|)<\infty$

for all compact sets $K\subset\mathbb{R}^{d}$ and all finite $T>0$. For the
(centered) Gaussian noise $\dot{W}$ we assume that the covariance is given by

$\mathbb{E}\left[\dot{W}(t,x)\dot{W}(s,y)\right]=\delta_{0}(t-s)\gamma(x-y),$

where $\delta_{0}$ denotes the Dirac delta and $\gamma$ are non-negative and
non-negative definite measures. That is, the noise can be described by a
centered Gaussian family with covariance

$\mathbb{E}[W(\phi)W(\psi)]=\int_{0}^{\infty}\int_{\mathbb{R}^{d}}\mathcal{F}\phi(\cdot,s)(\xi)\overline{\mathcal{F}(\psi)(\cdot,s)}(\xi)\hat{\gamma}(d\xi)ds,$
(2.2)

where $\hat{\gamma}$ is the spectral measure of $\gamma$, $\mathcal{F}$
denotes the Fourier transform, and $\phi,\psi$ are suitable functions. We
provide rigorous treatment on the construction of the Gaussian noise $\dot{W}$
in Section 3.1 that contains a brief introduction to Gaussian analysis and
stochastic integration required for our analysis.

We denote by $G_{t}(x)$ the fundamental solution (or the Green kernel)
associated to the operator $M=\partial_{t}-L_{x}$, in the sense that for any
"nice enough" intial data $u_{0}$ we have that

$\displaystyle u(x,t)=\int G_{t}(x-y)u_{0}(y)dy$

satisfies $Mu=0$ and $u\to u_{0}$ as $t\to 0$.

Consider now our original equation 2.1. For the initial condition
$u(0,x)=u_{0}(x)$, we assume that $u_{0}(x)$ is deterministic and satisfies,
for every $T>0$,

$\sup_{t\in(0,T],x\in\mathbb{R}^{d}}\int_{\mathbb{R}^{d}}|G_{t}(x-y)||u_{0}(y)|dy<\infty.$
(2.3)

We prove that under certain conditions, equation 2.1 admits a unique mild
solution in the following sense.

######  Definition 2.1.

We say that a random field $u(t,x)$, adapted to the filtration generated by
$\dot{W}$, is a mild solution to 2.1 if for all
$(t,x)\in\mathbb{R}_{+}\times\mathbb{R}^{d}$, the process $(s,y)\rightarrow
G_{t-s}(x-y)\sigma(u(s,y))\textbf{1}_{[0,t]}(s)$ is integrable and we have

$\displaystyle u(t,x)$ $\displaystyle=$
$\displaystyle\int_{\mathbb{R}^{d}}G_{t}(x-y)u_{0}(y)dy+\int_{0}^{t}\int_{\mathbb{R}^{d}}G_{t-s}(x-y)b(s,u(s,y))dyds.$
(2.4) $\displaystyle+$
$\displaystyle\int_{0}^{t}\int_{\mathbb{R}^{d}}G_{t-s}(x-y)\sigma(s,u(s,y))W(ds,dy).$

Here the stochastic integral exists in the sense of Dalang-Walsh, see Section
3.1.

The following existence and uniqueness result is the main result of this
paper. The proof follows standard arguments and is presented in Section 3.

######  Theorem 2.2.

Let $G_{s}$ be the fundamental solution to 2.1 and assume that there exists an
integrable function $g_{s}(x)=\mathcal{F}^{-1}(e^{-s\hat{f}})$,
$\hat{f}\geq-C$ such that $|G_{s}(x)|\leq g_{s}(x)$ and there exists
$\beta_{0}>2C$ such that,

$\int_{0}^{\infty}e^{-\beta_{0}s}\|g_{s}\|_{L^{1}(\mathbb{R}^{d})}ds<\infty$
(2.5)

and

$\int_{\mathbb{R}^{d}}\frac{1}{\beta_{0}+2\hat{f}(\xi)}\hat{\gamma}(d\xi)<\infty$
(2.6)

Then 2.1 admits a unique mild solution.

######  Remark 1.

Observe that the condition 2.5 is independent of the chosen covariance
$\gamma$, and is rather mild. Indeed, in many examples $g_{s}$ is given by a
density of some generating stochastic process at time $s$, and hence we have
$\|g_{s}\|_{L^{1}(\mathbb{R}^{d})}=1$ for all $s>0$. Moreover, by carefully
examining our proof one observes that 2.5 can be replaced with a weaker
condition $\int_{0}^{T}\|g_{s}\|_{L^{1}(\mathbb{R}^{d})}ds<\infty$ for all
$T>0$ finite. Finally, we note that 2.5 can be omitted in the case $b\equiv
0$.

######  Remark 2.

The proof of Theorem 2.2 actually gives more, however for purposes of
exposition we chose to present the simplified form 2.6. For a discussion about
further cases (including for instance the wave equation) see, Section 4.

######  Remark 3.

Note that if 2.6 is satisfied for some $\beta_{0}>0$, then it is automatically
valid for all $\beta>\beta_{0}$ as well.

###  2.1. Examples

In this section we outline some examples of quite general linear equations for
which Theorem 2.2 can be applied.

####  2.1.1. Linear operators with Gaussian upper bounds

Many parabolic and hypoelliptic evolution operators satisfy a Gaussian upper
bound for their fundamental solution, namely an estimate of the following type

$\displaystyle G_{s}(x)\leq g(t)e^{-\frac{C|x|_{\mathbb{G}}^{2}}{t}}$ (2.7)

where $\mathbb{G}$ is a homogeneous Lie group. The reason for writing it as
above is that if we consider $|x|_{\mathbb{G}}=|x|$, being the standard
Euclidean norm, we can cover a big class of parabolic operators. But for
hypoelliptic equations, we get a nontrivial norm induced by the corresponding
Lie-group.

Let us now consider some interesting examples satisfying 2.7 with the
Euclidean norm.

Let $L$ in $\mathbb{R}^{n}$ be given as

$\displaystyle
L=\sum_{i,j=1}^{n}a_{ij}(x,t)\partial_{ij}+\sum_{j=1}^{n}b_{i}(x,t)\partial_{i}+c(x,t)$

then if the matrix $a_{ij}(x,t)$ is uniformly elliptic, the functions $a,b,c$
are bounded and uniformly Hölder continuous with exponent $\alpha$, then there
exists positive constants $c_{1},c_{2}$ such that

$\displaystyle|G_{t}(x)|\leq c_{1}t^{-n/2}e^{-c_{2}\frac{|x|^{2}}{t}}$ (2.8)

where $G_{t}(x)$ is the fundamental solution to $\partial_{t}-L$, see [8, 9].
As such Theorem 2.2 is applicable, and since the Fourier transform of a
Gaussian is a Gaussian we get that if

$\displaystyle\int_{\mathbb{R}^{d}}\frac{1}{\beta_{0}+2C(c_{1})|\xi|^{2}}\hat{\gamma}(d\xi)<\infty$
(2.9)

then there exists a mild solution to 2.1. That is, the condition reduces to
the standard one for the heat equation, [5]. For instance, if $\gamma$ is
given by the Riesz kernel, i.e. $\hat{\gamma}=|\xi|^{\lambda-d}d\xi$, then the
above is verified as long as $\lambda<2$. In the case when $\gamma$ is in
$W^{k,1}(\mathbb{R}^{n})$ then

$\displaystyle|\hat{\gamma}|\leq\frac{C}{(1+|\xi|)^{k}},$

from this we see that 2.9 is verified if $k>n-2$. Finally we remark that if
$\gamma=\delta_{0}$ (white noise), then 2.9 only holds in dimension 1.

In fact the above can be extended as follows, if we instead consider the
divergence form operator

$\displaystyle L=\sum_{i,j=1}^{n}\partial_{i}(a_{ij}(x,t)\partial_{j}\cdot),$

with $a_{ij}$ still an elliptic matrix but now only bounded and measurable,
then $\partial_{t}-L$, Aronson [1] proved that 2.8 holds also in this case.

The case for non-divergence form equations with rough coefficients was treated
in [7].

Let us now consider the hypoelliptic setting and consider some relevant
examples. Let us begin with some notation. We consider the hypoelliptic
evolution operator

$\displaystyle M=\sum_{i=1}^{m}X_{i}^{2}+X_{0}-\partial_{t},$

where $X_{i}$ are smooth vector fields in $\mathbb{R}^{n}$ for $i=0,\ldots,m$,
and usually $m<n$. Such an equation induces an interesting geometry. Namely,
let us denote

$\displaystyle Y=X_{0}-\partial_{t},\quad\text{and}\quad\lambda\cdot
X=\lambda_{1}X_{1}+\ldots\lambda_{m}X_{m}$

then a curve in $\gamma:[0,R]\to\mathbb{R}^{n}\times[0,T]$ is $M$-admissible
if it is absolutely continuous and

$\displaystyle\gamma^{\prime}(s)=\lambda(s)\cdot
X(\gamma(s))+Y(\gamma(s)),\quad\text{a.e. in [0,R]}.$

That is, the curve has a direction at each point in the tangent-bundle
described by the vector-fields $X,Y$. If we assume that the vector fields are
such that

  * •

every two points $(x,t),(\xi,\tau)\in\mathbb{R}^{n}\times[0,T]$ can be
connected with an $M$ admissible curve,

  * •

there exists a homogeneous Lie group
$\mathbb{G}=(\mathbb{R}^{n}\times[0,T],\circ,\delta_{\lambda})$ for which
$X,Y$ are left translation invariant. Furthermore, that $X$ is
$\delta_{\lambda}$-homogeneous of degree 1 and $Y$ is
$\delta_{\lambda}$-homogeneous of degree 2.

Then the operator $M$ is hypoelliptic (distributional solutions are smooth),
see [3, 11]. This allows us to write $\mathbb{R}^{n}=V_{1}\oplus\ldots\oplus
V_{l}$, i.e. the vector field $X$ stratifies $\mathbb{R}^{n}$ into a direct
sum of subspaces, such that if we write $x=x^{(1)}+\ldots+x^{(l)}$, where
$x^{(k)}\in V_{k}$, then the dilations becomes simple multiplication in the
sense that

$\displaystyle\delta_{\lambda}(x,t)=(\lambda
x^{(1)}+\ldots+\lambda^{l}x^{(l)},\lambda^{2}t).$

Furthermore if we introduce the $\delta_{\lambda}$ homogeneous norm as

$\displaystyle|x|_{\mathbb{G}}=\max\left\\{|x_{i}^{(k)}|^{\frac{1}{k}},k=1,\ldots,l,i=1,\ldots,m_{k}\right\\}$

Then, for such operators $M$ we have that there exists (see [11]) a positive
constant $C$ such that

$\displaystyle|G_{t}(x)|\leq\frac{C}{t^{\frac{Q-2}{2}}}e^{-\frac{|x|_{\mathbb{G}}^{2}}{Ct}}.$

From this we see that we can apply Theorem 2.2 again in this context, however
now $\gamma$ needs to be adapted to the geometry of $\mathbb{G}$. Furthermore,
for certain hypoelliptic operators we even have positivity of the fundamental
solution, [2].

A prototypical example of a hypoelliptic operator is the following (well known
Kolmogorov operator)

$\displaystyle
K=\sum_{j=1}^{n}\partial_{x_{j}}^{2}+\sum_{j=1}^{n}x_{j}\partial_{y_{j}}-\partial_{t}$

where we are working in $\mathbb{R}^{2n}\times[0,T]$ and we denote the first
$n$ coordinates as $x_{i}$ and the other $n$ as $y_{i}$. The operator $K$ is
hypoelliptic in the above sense [3]. Furthermore, the fundamental solution for
$n=1$ is given by

$\displaystyle G_{t}(x,y)=\frac{\sqrt{3}}{2\pi
t^{2}}\exp\left(-\frac{x^{2}}{t}-\frac{3xy}{t^{2}}-\frac{3y^{2}}{t^{3}}\right)$

Considering the above "heat type" kernel in the finite interval $[0,T]$, it is
clear that we can bound it from above by
$g_{t}(x)=\frac{C}{t^{2}}\exp\left(-\frac{1}{C}\frac{x^{2}+y^{2}}{t}\right)$
for a positive constant $C$. As such the same conclusion as in Section 2.1.1
holds.

Operators of this type often appear in the context of SDE’s where the noise is
only in some directions, like the kinetic equations, for which the density
satisfies the kinetic Fokker-Planck equation. The kinetic Fokker-Planck is
hypoelliptic, see [13] and the references therein.

####  2.1.2. Fractional heat equation

Consider the fractional evolution equation defined as

$\displaystyle\partial_{t}u=-(-\Delta)^{s}u+Cu,\quad 0<s\leq 1,$

where $C=0$ corresponds to the fractional heat equation. The fourier transform
of said equation is

$\displaystyle\partial_{t}\hat{u}=-|\xi|^{2s}\hat{u}+C\hat{u}.$

As such, the fundamental solution is given by
$\hat{G}_{t}(\xi)=e^{-(|\xi|^{2s}-C)t}$, and we can thus apply Theorem 2.2
with $b\equiv 0$ and $\hat{f}\geq-C$, together with Remark 1.

####  2.1.3. Mixture operators

We formally consider an equation of the following type

$\displaystyle\partial_{t}u=L_{1}u+L_{2}u.$

If $G_{1},G_{2}$ are the fundamental solutions of

$\displaystyle\partial_{t}u=L_{1}u,\quad\partial_{t}u=L_{2}u,$

respectively, then $G=G_{1}\ast G_{2}$ is a fundamental solution of the above,
(under some integrability assumptions). Indeed, we have

$\displaystyle\partial_{t}(G_{1}\ast G_{2})$
$\displaystyle=(\partial_{t}G_{1}\ast
G_{2})+(G_{1}\ast\partial_{t}G_{2})=(L_{1}G_{1})\ast
G_{2}+G_{1}\ast(L_{2}G_{2})$ $\displaystyle=(L_{1}+L_{2})(G_{1}\ast G_{2}).$

Thus for instance we can combine Sections 2.1.1 and 2.1.2 and apply Theorem
2.2 to get the existence of a mild solution to 2.1.

##  3\. Proof of Theorem 2.2

###  3.1. Preliminaries on stochastic calculus

In this section we introduce stochastic analysis with respect to the noise
$\dot{W}$.

Denote by $C_{c}^{\infty}\left([0,\infty)\times\mathbb{R}^{d}\right)$ the
class of $C^{\infty}$ functions on $[0,\infty)\times\mathbb{R}^{d}$ with
compact support. We consider a Gaussian family of centered random variables

$\left(W(\varphi),\varphi\in
C^{\infty}_{c}\left([0,\infty)\times\mathbb{R}^{d}\right)\right)$

on some complete probability space $\left(\Omega,\mathcal{F},P\right)$ such
that

$\displaystyle\mathbb{E}[W(\varphi)W(\psi)]$ (3.1) $\displaystyle=$
$\displaystyle\int_{0}^{\infty}\int_{\mathbb{R}^{d}}\int_{\mathbb{R}^{d}}\varphi(s,y)\psi(s,y^{\prime})\gamma(y-y^{\prime})dydy^{\prime}ds:=\langle\varphi,\psi\rangle_{\mathfrak{H}}.$

By taking the Fourier transform, this can equivalently be written as 2.2. Note
that here $\gamma$ is not a function, and hence 3.1 should be understood as

$\mathbb{E}[W(\varphi)W(\psi)]=\int_{0}^{\infty}\int_{\mathbb{R}^{d}}\varphi(s,y)\left[\psi(s,\cdot)\ast\eta\right](y)dyds,$
(3.2)

where $\ast$ denotes the convolution. For the simplicity of our presentation,
we use notation $\gamma(y-y^{\prime})dydy^{\prime}$ in Sections 3.2 and 3.3
from which change of variable transformations are easier to follow.

We denote by $\mathfrak{H}$ the Hilbert space defined as the closure of
$C_{c}^{\infty}\left([0,\infty)\times\mathbb{R}^{d}\right)$ with respect to
the inner product (3.1). As a result, we obtain an isonormal process
$(W(\varphi),\varphi\in\mathfrak{H})$, which consists of a Gaussian family of
centered random variable such that, for every $\varphi,\psi\in\mathfrak{H}$,

$\mathbb{E}[W(\varphi)W(\Psi)]=\langle\varphi,\psi\rangle_{\mathfrak{H}}.$

The filtration $\mathbb{F}=(\mathcal{F}_{t})_{t\geq 0}$ associated to the
random noise $W$ is generated by random variables $W(\varphi)$ for which
$\varphi\in\mathfrak{H}$ has support contained in $[0,t]\times\mathbb{R}^{d}$.

Let us now define the stochastic integral with respect to $\dot{W}$. For every
random field $\\{X(s,y),s\geq 0,y\in\mathbb{R}^{d}\\}$ such that

$\mathbb{E}\|X\|_{\mathfrak{H}}^{2}=\mathbb{E}\int_{0}^{\infty}\int_{\mathbb{R}^{d}}X(s,y)X(s,y^{\prime})\gamma(y-y^{\prime})dydy^{\prime}<\infty,$

we can define the stochastic integral

$\int_{0}^{\infty}\int_{\mathbb{R}^{d}}X(s,y)W(ds,dy)$

in the sense of Dalang-Walsh (see, e.g. [5, 14]). It follows that we have the
Isometry

$\mathbb{E}\left[\int_{0}^{\infty}\int_{\mathbb{R}^{d}}X_{n}(s,y)W(ds,dy)\right]^{2}=\mathbb{E}\|X_{n}\|_{\mathfrak{H}}.$
(3.3)

Moreover, we have the following version of the Burkholder-Davis-Gundy
inequality: for any $t\geq 0$ and $p\geq 2$,

$\displaystyle\left|\left|\int_{0}^{\infty}\int_{\mathbb{R}^{d}}X(s,y)W(ds,dy)\right|\right|_{p}^{2}$
$\displaystyle\leq
c_{p}\int_{0}^{t}\int_{\mathbb{R}^{d}}\int_{\mathbb{R}^{d}}\|X(s,y)X(s,y^{\prime})\|_{\frac{p}{2}}\gamma(y-y^{\prime})dydy^{\prime}ds.$
(3.4)

###  3.2. Auxiliary results

Set

$I(s)=\|g_{s}\|_{L^{1}(\mathbb{R}^{d})}+\int_{\mathbb{R}^{d}}\int_{\mathbb{R}^{d}}g_{s}(y)g_{s}(y^{\prime})\gamma(y-y^{\prime})dydy^{\prime}.$
(3.5)

By taking the Fourier transform, we have

$I(s)=\|g_{s}\|_{L^{1}(\mathbb{R}^{d})}+\int_{\mathbb{R}^{d}}|\hat{g}_{s}(\xi)|^{2}\hat{\gamma}(d\xi)=\|g_{s}\|_{L^{1}(\mathbb{R}^{d})}+\int_{\mathbb{R}^{d}}e^{-2s\hat{f}(\xi)}\widehat{\gamma}(d\xi).$

######  Lemma 3.1.

Suppose that 2.5-2.6 hold for some $\beta_{0}$ and set

$\Upsilon(\beta)=\int_{0}^{\infty}e^{-\beta s}I(s)ds.$

Then $\Upsilon:(\beta_{0},\infty)\rightarrow(0,\infty)$ is well-defined and
non-increasing in $\beta$. Moreover, $\lim_{\beta\to\infty}\Upsilon(\beta)=0$.

###### Proof.

By Tonelli’s theorem, we have

$\int_{0}^{\infty}e^{-\beta
s}\int_{\mathbb{R}^{d}}e^{-2s\hat{f}(\xi)}\widehat{\gamma}(d\xi)ds=\int_{\mathbb{R}^{d}}\int_{0}^{\infty}e^{-\beta
s-2s\hat{f}(\xi)}ds\hat{\gamma}(d\xi)=\int_{\mathbb{R}^{d}}\frac{1}{\beta+2\hat{f}(\xi)}\hat{\gamma}(d\xi)$

for every $\beta$ such that $\beta+2\hat{f}(\xi)>0$. Since $I(s)$ is non-
negative, it is clear that $\Upsilon(\beta)$ is non-increasing which concludes
the proof. ∎

######  Remark 4.

We remark that for our purposes, it would actually suffice to consider

$\Upsilon_{T}(\beta)=\int_{0}^{T}e^{-\beta s}I(s)ds$

for each fixed $T<\infty$. Hence we could replace the condition 2.5 with
$\int_{0}^{T}\|g_{s}\|_{L^{1}(\mathbb{R}^{d})}ds<\infty$, cf. Remark 1.

The following Proposition is the main technical ingredient. The result follows
directly from [5, Lemma 15] adapted to our context.

######  Proposition 3.2.

Let $I(s)$ be given by 3.5 and suppose that 2.5-2.6 hold. Let $\iota>0$,
$\beta>\beta_{0}$, and $T\in(0,\infty)$ be fixed, and let $h_{n}$ be a
sequence of non-negative functions such that $\sup_{t\in[0,T]}h_{0}(t)<\infty$
and, for $n\geq 1$, we have

$h_{n}(t)\leq\iota\int_{0}^{t}h_{n-1}(s)e^{-\beta(t-s)}I(t-s)ds.$ (3.6)

Then the series

$H(\iota,p,t):=\sum_{n\geq 0}h_{n}(t)$ (3.7)

converges uniformly in $t\in[0,T]$.

###  3.3. Proof of Theorem 2.2

###### Proof of Theorem 2.2.

Let $T>0$ be fixed and finite. We consider the standard Picard iterations by
setting $u_{0}(t,x)=\int_{\mathbb{R}^{d}}G_{t}(x-y)u_{0}(y)dy$ and, for $n\geq
1$ and a given $\beta>\beta_{0}$,

$\begin{split}e^{-\beta t}u_{n+1}(t,x)&=e^{-\beta
t}u_{0}(t,x)+\int_{0}^{t}e^{-\beta
t}\int_{\mathbb{R}^{d}}G_{t-s}(x-y)b(s,u_{n}(s,y))dyds\\\
&+\int_{0}^{t}e^{-\beta
t}\int_{\mathbb{R}^{d}}G_{t-s}(x-y)\sigma(s,u_{n}(s,y))W(ds,dy),\hskip
14.22636ptt\geq 0,x\in\mathbb{R}^{d}.\end{split}$

We first prove that, for each $n\geq 1$, $u_{n}(t,x)$ is well-defined and
satisfies

$\sup_{t\in(0,T]}\sup_{x\in\mathbb{R}^{d}}\mathbb{E}\left[e^{-p\beta
t}\left|u_{n}(t,x)\right|^{p}\right]<\infty.$ (3.8)

We first note that if 3.8 holds for some $n\geq 0$, it follows from the
Lipschitz continuity of $\sigma$ that

$\begin{split}\mathbb{E}\left[e^{-p\beta
s}|\sigma(s,u_{n}(s,y))|^{p}\right]&\leq C\left[\mathbb{E}\left[e^{-p\beta
s}|\sigma(s,u_{0}(s,y))|^{p}\right]\right.\\\
&+\left.\mathbb{E}\left[e^{-p\beta
s}|u_{n}(s,y))-u_{0}(s,y)|^{p}\right]\right].\end{split}$

By the boundedness of $\sigma$ on compacts and 2.3, we get

$\sup_{s\in[0,T],y\in\mathbb{R}^{d}}\mathbb{E}\left[e^{-p\beta
s}|\sigma(s,u_{0}(s,y))|^{p}\right]=\sup_{s\in[0,T],y\in\mathbb{R}^{d}}e^{-p\beta
s}|\sigma(s,u_{0}(s,y))|^{p}<\infty$

and consequently,

$\mathbb{E}\left[e^{-p\beta s}|\sigma(s,u_{n}(s,y))|^{p}\right]\leq
C\left(1+\sup_{s\in[0,T],y\in\mathbb{R}^{d}}\mathbb{E}\left[e^{-p\beta
s}|u_{n}(s,y)|^{p}\right]\right)<\infty.$

This in turn gives us, by Hölder’s inequality, that

$\sup_{s\in[0,T]}\sup_{y,y^{\prime}\in\mathbb{R}^{d}}e^{-\beta
s}\left|\left|\sigma(s,u_{n}(s,y))\sigma(s,u_{n}(s,y^{\prime}))\right|\right|_{\frac{p}{2}}<\infty.$

Exactly the same way, we obtain

$\sup_{s\in[0,T]}\sup_{y\in\mathbb{R}^{d}}\mathbb{E}\left[e^{-p\beta
s}|b(s,u_{n}(s,y))|^{p}\right]<\infty.$

In view of 3.3 and 3.4 together with the boundedness of $u_{0}(t,x)$ and
$|G_{t}(y)|\leq g_{t}(y)$, applying the Minkowski integral inequality leads to

$\displaystyle\mathbb{E}\left[e^{-p\beta
t}\left|u_{n+1}(t,x)\right|^{p}\right]$ $\displaystyle\leq$ $\displaystyle
C\left(\left[e^{-\beta t}u_{0}(t,x)\right]^{p}+\Big{\|}\int_{0}^{t}e^{-\beta
t}\int_{\mathbb{R}^{d}}G_{t-s}(x-y)b(s,u_{n}(s,y))dyds\Big{\|}_{p}^{p}\right.$
$\displaystyle+$ $\displaystyle\left.\Big{\|}\int_{0}^{t}e^{-\beta
t}\int_{\mathbb{R}^{d}}\int_{\mathbb{R}^{d}}G_{t-s}(x-y)G_{t-s}(x-y^{\prime})\right.$
$\displaystyle\phantom{kukkuu}\left.\times\sigma(s,u_{n}(s,y))\sigma(s,u_{n}(s,y^{\prime}))\gamma(y-y^{\prime})dy^{\prime}dyds\Big{\|}_{\frac{p}{2}}^{\frac{p}{2}}\right)$
$\displaystyle\leq$ $\displaystyle
C\left[1+\left(\int_{0}^{t}e^{-\beta(t-s)}\int_{\mathbb{R}^{d}}g_{t-s}(x-y)e^{-\beta
s}\|b(s,u_{n}(s,y))\|_{p}dyds\right)^{p}\right.$ $\displaystyle+$
$\displaystyle\left.\left(\int_{0}^{t}e^{-\beta(t-s)}\int_{\mathbb{R}^{d}}\int_{\mathbb{R}^{d}}g_{t-s}(x-y)g_{t-s}(x-y^{\prime})\right.\right.$
$\displaystyle\phantom{kukkuu}\left.\left.\times e^{-\beta
s}\left|\left|\sigma(s,u_{n}(s,y))\sigma(s,u_{n}(s,y^{\prime}))\right|\right|_{\frac{p}{2}}\gamma(y-y^{\prime})dy^{\prime}dyds\right)^{\frac{p}{2}}\right]$
$\displaystyle\leq$ $\displaystyle
C\left[1+\left(\int_{0}^{t}e^{-\beta(t-s)}I(t-s)ds\right)^{p}+\left(\int_{0}^{t}e^{-\beta(t-s)}I(t-s)ds\right)^{\frac{p}{2}}\right]$
$\displaystyle\leq$ $\displaystyle
C\left(1+\Upsilon^{p}(\beta)+\Upsilon^{\frac{p}{2}}(\beta)\right).$

which is finite for $\beta>\beta_{0}$. Since $u_{0}(t,x)$ is uniformly bounded
on $[0,T]\times\mathbb{R}^{d}$, 3.8 for all $n\geq 0$ thus follows from
induction and hence $u_{n+1}$ is well-defined. The same arguments applied to

$H_{n}(t):=\sup_{x\in\mathbb{R}^{d}}\left[e^{-\beta
pt}\mathbb{E}\left[\left|u_{n+1}(t,x)-u_{n}(t,x)\right|^{p}\right]\right]^{\frac{1}{p}}$

gives us

$\displaystyle H_{n}(t)$ $\displaystyle\leq
C\left[\int_{0}^{t}e^{-\beta(t-s)}\int_{\mathbb{R}^{d}}g_{t-s}(x-y)H_{n-1}(s)dyds\right.$
$\displaystyle+\left(\int_{0}^{t}\int_{\mathbb{R}^{d}}\int_{\mathbb{R}^{d}}e^{-2\beta
t}g_{t-s}(x-y)g_{t-s}(x-y^{\prime})\left|\left|\sigma(s,u_{n}(s,y))-\sigma(s,u_{n-1}(s,y))\right|\right|_{p}\right.$
$\displaystyle\phantom{kukkuu}\times\left.\left.\left|\left|\sigma(s,u_{n}(s,y^{\prime}))-\sigma(s,u_{n-1}(s,y^{\prime}))\right|\right|_{p}\gamma(y-y^{\prime})dy^{\prime}dyds\right)^{\frac{1}{2}}\right]$
$\displaystyle\leq
C\left[\int_{0}^{t}e^{-\beta(t-s)}I(t-s)H_{n-1}(s)ds+\left(\int_{0}^{t}e^{-2\beta(t-s)}I(t-s)H^{2}_{n-1}(s)ds\right)^{\frac{1}{2}}\right].$

Here Hölder’s inequality gives

$\left(\int_{0}^{t}e^{-2\beta(t-s)}I(t-s)H^{2}_{n-1}(s)ds\right)^{\frac{1}{2}}\leq
C\int_{0}^{t}e^{-\beta(t-s)}I(t-s)H_{n-1}(s)ds$

leading to

$H_{n}(t)\leq C\int_{0}^{t}e^{-\beta(t-s)}I(t-s)H_{n-1}(s)ds.$

By noting that $b(s,u_{0}(s,y))$ and $\sigma(s,u_{0}(s,y))$ are uniformly
bounded over $[0,T]\times\mathbb{R}^{d}$, we obtain from the above
computations that $\sup_{s\in[0,T]}H_{0}(s)<\infty$. Consequently, it follows
from Proposition 3.2 that $\sum_{n\geq 1}H_{n}(t)$ converges uniformly on
$[0,T]$, and hence the sequence $u_{n}$ converges in $L^{p}(\Omega)$ to some
process $u\in L^{p}(\Omega)$, uniformly on $[0,T]\times\mathbb{R}^{d}$. From
the uniform convergence one can deduce further that $u$ satisfies 2.4, and
thus we have obtained the existence of the mild solution. Proving the
uniqueness in a similar fashion and noting that $T>0$ was arbitrary concludes
the proof. ∎

##  4\. Concluding remarks

In this article we have considered general SPDEs 2.1 and showed that one can
obtain existence and uniqueness of the mild solution by classical arguments,
initiated by Dalang [5], for a very large class of operators $L_{x}$. Our
condition 2.6 is similar to the Dalang condition 1.2 for the stochastic heat
equation and actually, our examples in Section 2.1 reveals that 1.2 is indeed
sufficient for a very large class of operators $L_{x}$. By carefully examining
the above proofs, we observe that even more is true. Indeed, condition 2.6
could be replaced with a condition $|G_{s}(x)|\leq g_{s}(x)$ and, for all
$T>0$,

$\int_{0}^{T}\int_{\mathbb{R}^{d}}e^{-\beta
s}[\hat{g}_{s}(\xi)]^{2}\hat{\gamma}(d\xi)ds<\infty.$

This formulation allows to consider even more general operators $L_{t}$ in the
time variable instead of considering merely $L_{t}=\partial_{t}$. In
particular, with this formulation one can easily recover the existence and
uniqueness result of [5] for the stochastic wave equation (in one or two
dimensions) with $L_{t}=\partial_{tt}$. As another example, we can obtain the
existence and uniqueness result in the case of the fractional power of the
full heat operator $(\partial_{t}-\Delta)^{s}$, studied for instance in [12].

## References

  * [1] T. Aronson. Bounds for the fundamental solution of a parabolic equation.Bulletin of the American Mathematical society 73.6 (1967): 890-896.
  * [2] B. Avelin, L. Capogna, G. Citti and K. Nyström. Harnack estimates for degenerate parabolic equations modeled on the subelliptic p-Laplacian. Adv. Math. 257 (2014), 25–65.
  * [3] A. Bonfiglioli, E. Lanconelli, and F. Uguzzoni. Stratified Lie groups and potential theory for their sub-Laplacians. Springer Science & Business Media, 2007.
  * [4] L. Chen and R. Dalang. Moments, intermittency and growth indices for nonlinear stochastic fractional heat equation. Stoch. Partial Differ. Equ. Anal. Comput., 3(3): 360–397, 2015.
  * [5] R. Dalang. Extending the Martingale Measure Stochastic Integral With Applications to Spatially Homogeneous S.P.D.E.’s. Electron. J. Probab. Volume 4, paper no. 6, 29 pp., 1999.
  * [6] L. Debbi and M. Dozzi. On the solutions of nonlinear stochastic fractional partial differential equations in one spatial dimension. Stoch. Proc. Appl., 115: 1761–1781, 2005.
  * [7] L. Escauriaza. Bounds for the fundamental solutions of elliptic and parabolic equations: In memory of eugene fabes. Communications in Partial Differential Equations, 25.5-6 (2000): 821–845.
  * [8] A. Friedman. Partial differential equations of parabolic type. Courier Dover Publications, 2008.
  * [9] A. Friedman. Stochastic differential equations and applications. Vol. 1. Probability and Mathematical Statistics, Vol. 28. Academic Press [Harcourt Brace Jovanovich, Publishers], New York-London, 1975.
  * [10] M. Foondun and D. Khoshnevisan. On the stochastic heat equation with spatially-colored random forcing. Transactions of the American Mathematical Society, 365: 409–458, 2013.
  * [11] A.E. Kogoj and E. Lanconelli. An invariant Harnack inequality for a class of hypoelliptic ultraparabolic equations. Mediterranean Journal of Mathematics 1.1 (2004): 51–80.
  * [12] K. Nyström and O. Sande. Extension properties and boundary estimates for a fractional heat operator. Nonlinear Analysis 140 (2016): 29–37.
  * [13] C. Villani, Hypocoercivity. No. 949-951. American Mathematical Soc., 2009.
  * [14] J. B. Walsh. An Introduction to Stochastic Partial Differential Equations. In: École d’été de probabilités de Saint-Flour, XIV—1984, 265-439. Lecture Notes in Math. 1180, Springer, Berlin, 1986.


# Resolution Complete In-Place Object Retrieval
given Known Object Models

Daniel Nakhimovich, Yinglong Miao, and Kostas E. Bekris  The authors are with
the Dept. of Computer Science, Rutgers, New Brunswick, NJ. Email: {d.nak,
<EMAIL_ADDRESS>The work is partially supported by NSF awards
1934924 and 2021628. The opinions expressed here are those of the authors and
do not reflect the positions of the sponsor.

###### Abstract

This work proposes a robot task planning framework for retrieving a target
object in a confined workspace among multiple stacked objects that obstruct
the target. The robot can use prehensile picking and in-workspace placing
actions. The method assumes access to 3D models for the visible objects in the
scene. The key contribution is in achieving desirable properties, i.e., to
provide (a) safety, by avoiding collisions with sensed obstacles, objects, and
occluded regions, and (b) resolution completeness ($\tt RC$) - or
probabilistic completeness ($\tt PC$) depending on implementation - which
indicates a solution will be eventually found (if it exists) as the resolution
of algorithmic parameters increases. A heuristic variant of the basic $\tt RC$
algorithm is also proposed to solve the task more efficiently while retaining
the desirable properties. Simulation results compare using random picking and
placing operations against the basic $\tt RC$ algorithm that reasons about
object dependency as well as its heuristic variant. The success rate is higher
for the $\tt RC$ approaches given the same amount of time. The heuristic
variant is able to solve the problem even more efficiently than the basic
approach. The integration of the $\tt RC$ algorithm with perception, where an
RGB-D sensor detects the objects as they are being moved, enables real robot
demonstrations of safely retrieving target objects from a cluttered shelf.

##  I INTRODUCTION

Robotic manipulation has the potential of being integrated into daily lives of
people, such as in household service areas [1, 2]. A useful skill for such
household settings involves the retrieval of a target object from a confined
and cluttered workspace, such as a fridge or a shelf, which may also require
the rearrangement of other objects in the process. In this context, it is
important to consider how to safely retrieve objects while minimizing the time
spent or the amount of pick and place operations, so as to assist humans
efficiently.

Figure 1: (Top) Setup for the real demonstration using an RGB-D sensor,
robotiq gripper, and Yaskawa Motoman robot to retrieve the target bottle.
(Bottom Left) The camera view in which objects are occluded. (Bottom Right)
The corresponding voxel map.

One of the challenging aspects of these problems that requires explicit
reasoning relates to heavy occlusions in the scene, as the sensor is often
mounted on the robot and has limited visibility. These visibility constraints
complicate the task planning process, as rearranging one object can limit
placements for others and can introduce new occlusions. Moreover, real-world
scenes in household setups are often unstructured and involve objects with
complex spatial relationships, such as objects stacked on each other.

Many previous efforts on object retrieval have focused on cases where blocking
objects are extracted from the workspace [3, 4], which simplifies the
challenge as it does not require identifying temporary placement locations for
the objects within the confined space. In-place rearrangement has been
considered in some prior efforts [5]. While this prior method is efficient, it
is not complete as it limits the reasoning on the largest object in the scene
to analyze object traversability [4]. Alternatives use machine learning to
guide the decision making [6, 7, 8], which is an exciting direction but does
not easily allow for performance guarantees, such as resolution completeness.
Setups where object stacking arise have received less attention and most
solutions that do consider stacking are dependent on machine learning for
reasoning [9, 10, 11]. Some works have proposed testbeds [12] that can help
evaluate solutions in this domain.

This work focuses on object retrieval in clutter where occlusions arise and
objects may be initially stacked under the assumption of known object models
(e.g. Figure 2). It aims at a theoretical understanding to show the algorithm
has safety and $\tt RC$ guarantees. A heuristic variant improves the practical
efficiency. Key features of the proposed $\tt RC$ framework are the following:

  * $\bullet$

it employs an adaptive dependency graph data structure inspired by solutions
in object rearrangement with performance guarantees [13] that express a larger
variety of object relationships than previously considered (namely occlusion
dependencies);

  * $\bullet$

it computes the occlusion volume of detected objects as a heuristic to inform
the planning process;

  * $\bullet$

it reasons about the collision-free placement of objects in the confined
workspace efficiently by utilizing a voxelized representation of the space;

  * $\bullet$

it achieves $\tt RC$ (or $\tt PC$) depending on the implementation of the
underlying sampling subroutines;

  * $\bullet$

it provides an early termination criterion when a solution cannot be found for
the given resolution.

Simulation results, using a model of a Yaskawa Motoman manipulator for
rearranging objects on a tabletop as shown in Fig. 1, evaluate the proposed
$\tt RC$ framework against a baseline using random picking and placing
operations. Both variants of the $\tt RC$ reason about object dependencies.
One doesn’t use heuristics and one is heuristically guided while retaining
$\tt RC$. Both of the $\tt RC$ approaches outperform the baseline. The
heuristically guided solution is able to solve the problem more efficiently
than the basic $\tt RC$ solution. The integration of the proposed approach
with perception, where an RGB-D sensor detects the objects as they are being
moved, provides real robot demonstrations of safe object retrieval from a
cluttered shelf.

##  II RELATED WORK

Some works on object retrieval rely on geometric analysis of object occlusion
[3, 4]. They provide theoretical insights but frequently do not limit actions
to in-place rearrangement of blocking objects. Specifically, one method
constructs a dependency graph taking into account objects that jointly occlude
a region and objects that block others [3]. The occlusion volume is used to
estimate belief regarding the target object position and helps to construct an
optimal A* algorithm. An alternative constructs a Traversability graph
(T-graph) [4], where the edges encode if the largest object in the scene can
be moved between two poses. It then constructs an algorithm to extract the
target object, but is limited as the traversability edges are too
constraining. The POMDP formulation is popular for the task [14, 15], which
allows the application of general POMDP solvers. The POMDP formulation was
also adopted by the work that formalizes object retrieval in unstructured
scenes as ”mechanical search” [16].

Alternatives rely on learning-based methods to solve such challenges, such as
reinforcement learning [6] or target belief prediction [7, 9, 17]. They report
good performance but do not provide theoretical guarantees given the black-box
nature of the solutions. In particular, a reinforcement learning solution [6]
uses the rendered top-down projections of the scene to predict the target
poses. A recent follow-up effort [7] on previous work [17] estimates the 1D
position belief of the target object on the shelf via machine learning. It
then constructs a policy based on the distribution change after applying
pushing and suction actions. It incorporates stacking and unstacking actions,
where object stacking is represented by a tree structure. Other works such as
[18] utilize learning for planning grasps to greedily empty bins of complex
and novel objects.

A related work [19] proposes a complete framework to safely reconstruct all
objects in the scene amidst object occlusions. Nevertheless, object retrieval
may not require reconstructing all objects and requires a search procedure
that is more task-driven for efficiency. There are also previous works [20,
19] that construct a voxelization of the environment to model object
occlusions, similar to the current work. This representation is used to
compute an object’s occlusion volume, which provides heuristic guidance.
Object spatial relationships are often represented by scene graphs [10, 11],
or implicitly in machine learning solutions [16, 21, 22].

What stands out in this work is that it proposes a general template for a $\tt
RC$ or $\tt PC$ approach to task retrieval in occluded environments that only
relies on basic motion and perception primitives. This modular nature allows
for quick sim-to-real transfer and passive performance improvement as the
primitives are improved over time. Furthermore, an efficient implementation is
demonstrated utilizing a voxelized representation of the environment for quick
collision filtering of object placements as well as providing an effective
heuristic to rank object manipulations. Thus, the framework enables effective
in-workspace manipulation.

##  III PROBLEM STATEMENT

Consider an environment with a set
$\mathbb{S}=\\{o_{1},...,o_{n}\\}\subset\mathbb{O}$ of $n$ objects for which
there are available 3D models. The objects are stably resting on a support
surface. Objects are allowed to be initially stacked and occlude each other
from the camera view.

The robot has one fixed RGB-D sensor at its disposal. Discovered objects are
those recognized given the observation history. An objects is assumed to be
recognized once an image segmentation process recognizes it as an individual
object in the observed image. Similarly, a perception method for detecting the
target object once observed is assumed. The region of the workspace occluded
by object $o_{i}$ at pose $s_{i}$ is denoted as $O_{i}(s_{i})$. Similarly the
space uniquely occluded by object $o_{i}$ at pose $s_{i}$, called the direct
occlusion space, is denoted as $\tilde{O}_{i}(s_{i})$. The proposed algorithms
gradually removes occlusions and recognizes objects. A motion planner is used
to plan pick-and-place actions.

While objects can start out stacked, they are not re-stacked and are only
placed on the ground surface during actions. While objects can start out
stacked, no reasoning about stability and ability to re-stack objects is
considered. Thus, once an action to pick up a stacked object is taken, that
object will only be placed on the ground surface. For further assumptions and
required properties of the motion planner see section V.

The objective for the object retrieval task is to determine a sequence of pick
and place actions in order to discover and subsequently retrieve the target
object; the target need not be directly visible or pickable from the robot’s
sensors. The corresponding solution should provide desirable guarantees: (a)
safety, by avoiding collisions with sensed obstacles and objects as well as
occluded regions, and (b) resolution completeness (RC) - or alternatively
probabilistic completeness, depending on the implementation of the underlying
motion planner, grasping process and object placement sampling. The
optimization objective is to minimize the number of performed actions until
the target object is retrieved.

##  IV METHOD

The proposed pipeline is detailed in Algorithm 1. First a voxelized
representation of the scene and a dependency graph are computed (lines 3,4,5),
which are detailed in subsection IV-A. The dependency graph contains a belief
state of the current scene based on visibility and reachability constraints.
All the sinks of this directed graph represent likely pickable objects. The
ranks are described later (see subsection IV-C). If the target object is
pickable then it is retrieved and the pipeline is terminated line 6.
Otherwise, a placement is planned (see subsection IV-B) one at a time for each
object in a shuffled order biased by the ranks (TryMoveOne line 7). The first
successful plan is executed.

If no placement is found for any of the pickable objects, then a fallback
procedure is called to try and pick one object and move it temporarily out of
view of the camera; the first successful plan is executed, the scene is re-
sensed, and a new placement is sampled for the object or the object is simply
put back at the same spot (MoveOrPlaceback line 10). At this stage, the
pipeline could restart if a new object was discovered (lines 12-15). Otherwise
the same set of pickable objects are tested for new placements in the scene
(line 16). If these two operations of “moving an object to look behind
it”(MoveOrPlaceback) followed by retrying to “move one of the pickable objects
to a new spot”(TryMoveOne) fail for all objects, then the pipeline can return
and report failure for the current resolution (line 22). If the sequence of
operations succeeds, then the pipeline can restart (line $19\rightarrow 2$).

Algorithm 1 $\tt RC\\_Pipeline$($\tt target$)

1:$\tt failure\leftarrow false$

2:while $\tt failure=false$ do

3: $\tt space\leftarrow UpdateVoxelsFromImage()$

4: $\tt dg\leftarrow{\tt DepGraph(space)}$

5: $\tt sinks,ranks\leftarrow RankSinks(target,dg)$

6: if $\tt target\in sinks$ then $\tt break$

7: if $\tt{\tt TryMoveOne(sinks,ranks)}=false$ then

8: $\tt failure\leftarrow true$

9: for $\tt sink\in sinks$ do

10: if $\tt{\tt MoveOrPlaceback(sink)}=false$ then

11: $\tt continue$

12: $\tt space\leftarrow UpdateVoxelsFromImage()$

13: if $\tt DidDiscoverObject(space)$ then

14: $\tt failure\leftarrow false$

15: $\tt break$

16: if $\tt{\tt TryMoveOne(sinks,\emptyset)}=false$ then

17: $\tt continue$

18: $\tt failure\leftarrow false$

19: $\tt break$

20:if $\tt not~{}failure$ then

21: $\tt Retrieve(target)$

22:$\tt return~{}failure$

\begin{overpic}[width=145.69568pt,trim=0.0pt 140.525pt 40.15pt
22.08249pt,clip]{figures/sim_exp/step1_ws.png}
\put(3.0,10.0){\huge$\displaystyle 1)$} \end{overpic}

\begin{overpic}[width=145.69568pt,trim=0.0pt 140.525pt 40.15pt
22.08249pt,clip]{figures/sim_exp/step2_ws.png}
\put(3.0,10.0){\huge$\displaystyle 2)$} \end{overpic}

\begin{overpic}[width=145.69568pt,trim=0.0pt 140.525pt 40.15pt
22.08249pt,clip]{figures/sim_exp/step3_ws.png}
\put(3.0,10.0){\huge$\displaystyle 3)$} \end{overpic}

\begin{overpic}[width=145.69568pt,trim=0.0pt 140.525pt 40.15pt
22.08249pt,clip]{figures/sim_exp/step4_ws.png}
\put(3.0,10.0){\huge$\displaystyle 4)$} \end{overpic}

\begin{overpic}[width=145.69568pt,trim=0.0pt 140.525pt 40.15pt
22.08249pt,clip]{figures/sim_exp/step5_ws.png}
\put(3.0,10.0){\huge$\displaystyle 5)$} \end{overpic}

\begin{overpic}[width=145.69568pt,trim=0.0pt 140.525pt 40.15pt
22.08249pt,clip]{figures/sim_exp/step6_ws.png}
\put(3.0,10.0){\huge$\displaystyle 6)$} \end{overpic}

\begin{overpic}[width=106.23477pt]{figures/sim_exp/step1_dg_c.png}
\put(-1.0,5.0){\huge a)} \end{overpic}

\begin{overpic}[width=106.23477pt]{figures/sim_exp/step2_dg.png}
\put(-1.0,5.0){\huge b)} \end{overpic}

\begin{overpic}[width=106.23477pt]{figures/sim_exp/step3_dg.png}
\put(-1.0,5.0){\huge c)} \end{overpic}

\begin{overpic}[width=106.23477pt]{figures/sim_exp/step4_dg.png}
\put(-1.0,5.0){\huge d)} \end{overpic}

Figure 2:  Images (1)-(6) show a simulated experiment from initial
configuration to final one action at a time with corresponding camera views in
the top left. The corresponding generated dependency graphs transitioning
between the images (1)-(2), (2)-(3), (3)-(4), and (4)-(5) are shown in images
(a)-(d). The colors of nodes corresponds with the objects in the scene and the
red object (labeled ‘T’) is the target object for the trial. The last graph
between images 5-6 is not shown since it is trivial having no dependencies
between any objects.

###  IV-A Voxel Map and Dependency Graph

A 3D occlusion voxel grid within the workspace is constructed from RGB-D
images. First the point cloud (in world frame) is generated using the RGB-D
image and the inverse of the camera projection. These points are down-sampled
into a voxelization of the scene. Given segmentation image information, each
object is associated with a portion of the voxel grid. Object geometry is used
to label voxels as occupied and remaining associated voxels as occluded. The
occluded regions of objects may intersect when jointly occluded.

The dependency graph is a directed graph where each node represents a visible
(or target) object and a labeled edge $(o_{i},o_{j},r)$ represents a relation
$r$ from object $o_{i}$ to $o_{j}$ that necessitates $o_{j}$ to be picked and
placed before $o_{i}$ could be picked. Valid relations in this work include
“below”, “grasp blocked by”, and – for the prediction of the target object –
“hidden by”. See Figure 2 for a sequence of such dependency graphs generated
during an example experiment.

Object x is defined to be “below” object y ($x\xrightarrow{below}y$) if object
x touches object y and the z-coordinate of the center of mass of x is less
than that of y. Note that this isn’t guaranteed to capture all intuitive cases
of one object being below another for non-convex objects. This relation is
computed using object models and poses given by the perception system.

Object x has its “grasp blocked by” y ($x\xrightarrow{blocked}y$) if there are
no collision free grasp poses for object x and the arm is in collision with y
for one or more valid grasp poses. (Grasp poses are sampled and tested by
inverse kinematics (IK) for discovered objects) Grasp poses are sampled using
inverse kinematics (IK) to discovered objects. (or, if the grasping pose
collides with objects, an edge to each collided object is added) If there
exists a collision free grasp for an object, no blocking edges are added;
otherwise, an edge to each object that has a collision with the arm is added.
Note that although this relation guarantees that the source object blocks the
target, it doesn’t capture all such reachability dependencies. This however is
not an issue for completeness as Algorithm 1 will eventually try to grasp all
objects in the case of motion planning failure.

The target object t is possibly “hidden by” x ($t\xrightarrow{hidden}x$) if
the target isn’t sensed in the scene and object x is touching the table. This
relation is used to keep track of the belief state of where the target is.
Each edge is assigned a probability based on the volume of the occluded space
behind object x (see subsection IV-C).

###  IV-B Placement Sampling

A valid placement is one that doesn’t collide with another object or any
undiscovered area of the workspace. Instead of randomly sampling x,y
coordinates for an object and checking for collision at that point, we create
a grid matching the horizontal extents of the workspace and add a collision
mask which is the shadow of the object occupancy and occlusion voxel grid
looking from a birds eye view. This mask is then convolved with the shadow of
the object that is to be placed. The occupied pixel indices indicate
collision-free placements and are converted to world coordinates. Object
orientations can be enumerated by rotating the object shadow.

###  IV-C Target Object Prediction

Intelligent object location prediction is achieved by applying a heuristic
which ranks the pickable objects determined by the dependency graph (line 4 in
Algorithm 1). This is done by augmenting the dependency graph edges with a
weight $p\in(0,1]$ estimating the probability for the relation
$o_{i}\overset{r}{\longrightarrow}o_{j}$ to be true. To rank a pickable
object, the sum of products of edge weights of all simple paths between the
target and the object is computed. When sampling from the list of pickable
objects, this rank is used as a probability weight.

For the “below” relation $p=1$ since object segmentation is assumed to be
reliable. For the “grasp blocked by” relation p is equal to the fraction of
total sampled grasps for which that object is in collision with the arm. Note,
the weights of grasp blocking edges coming out of any object need not add to
one (hence don’t truly represent probabilities) since the arm could be
colliding with multiple objects for any single grasp. For the “hidden by”
relation, the goal is to encourage knowledge gain of the environment. This is
done by normalizing the volume of the direct occlusion region of each stack of
objects and assigning the inverse as the probability estimate that the target
is hidden behind each stack. This heuristic biases the pipeline towards
discovering large volumes of occluded workspace.

From an algorithmic point of view, there is technically no reason to normalize
the output of the heuristic, however, representing the heuristics as
probabilities is insightful since - without prior knowledge - the probability
the hidden object to be in a larger volume is larger than the probability of
it being in a smaller volume. Furthermore, modeling the dependency graph edges
with probabilities as opposed to non-normalized weights is conducive for
exploring future work which might seek to combine the probabilities based on
the proposed volumetric-heuristics with priors based on the semantics of the
objects involved or from an additional human instruction (see section VII
work).

##  V RESOLUTION COMPLETENESS

Given a (formally) complete motion planner (finds solution in finite time if
exists), a continuous space grasp sampler, and continuous placement sampler,
the proposed pipeline would be $\tt PC$ Given a $\tt RC$ motion planner, a
discrete space grasp sampler, and discrete placement sampler, the proposed
pipeline would be $\tt RC$.

To show the $\tt PC$ or $\tt RC$ of the algorithm proposed in this work, a
simpler version of the algorithm is analysed first. Without loss of
generality, the actual proposed algorithm will be likewise proven complete.

Consider a much simpler algorithm that, at every iteration, tries to pick an
object at random and, if is not the target object, subsequently place it
randomly in the explored region of the workspace; call this RAND-ACT.

###### Lemma V.1

RAND-ACT is $\tt PC$ or $\tt RC$ (depending on the planning and sampling
subroutines).

###### Proof:

At every iteration, RAND-ACT attempts to perform a random action.
Consequently, RAND-ACT executes a random walk on the space of all actions.
Since pick and place actions are reversible, if a sequence of such valid
actions exists, this algorithm will eventually perform it (or an augmented
version of the sequence - i.e. placing an object back to where it was picked
from) in the limit (or finite time for RC subroutines). ∎

###### Corollary V.2

Algorithm 1 is $\tt PC$ or $\tt RC$ (depending on the planning and sampling
subroutines)).

###### Proof:

Indeed Algorithm 1 is really a fancy implementation of RAND-ACT. At every
iteration, the dependency graph is used to identify (and heuristically rank)
the currently pickable objects. One of these is chosen randomly for a pick and
place action via the $\tt TryMoveOne$ subroutine. The $\tt MoveOrPlaceback$
subroutine acts as a fallback in case there are no new discovered placements
found; it delays placement sampling till after the environment is re-sensed
with the picked object moved out of the way. Notice further, that even if the
pickable objects are sampled weighted according to their ranking rather than
uniformly, the action space is still explored entirely because each action has
a positive probability of being sampled. Thus, w.l.o.g. Algorithm 1 is $\tt
PC$ or $\tt RC$ as well. ∎

Failure Detection: The implementation in this work uses the $\tt RC$ approach.
In addition to $\tt RC$, Algorithm 1 takes a step closer towards achieving
general completeness by actually detecting certain unsolvable cases within the
completeness constraints of the motion and sampling subroutines. The
detectable unsolvable instances for which the algorithm will return failure in
finite time are as follows.

  * $\bullet$

No object can be grasped. This could happen if two objects each block the
grasp of the other.

  * $\bullet$

No objects can be placed anywhere (except its current spot). This could happen
in a highly cluttered scene where the only valid placement for each object is
to just put it back where it was.

Thus, Algorithm 1 has stronger guarantees than $\tt RC$ but is not formally
complete since it may run forever by juggling two objects between two
placements.

A Caveat: A fundamental assumption in the argument presented is that actions
are reversible. This is not always true in practice depending on the
implementation of the sampling subroutines. And indeed, the implementation of
the placement sampling process proposed in subsection IV-B implies
irreversible actions for scenes with stacked objects because it does not
consider the possibility of re-stacking objects. Thus, the proposed pipeline
(as implemented) is only resolution complete on the any sub-task where the
objects are no longer stacked. Implementation of stacking actions is planned
for future work.

##  VI EXPERIMENTS

Figure 3: Execution on the real robot : (a) Initial scene where the red bottle
is hidden. (b) The robot moves the yellow bottle, which occludes the most
space. (c) The robot moves the second yellow bottle, revealing the red bottle.
(d and e) The robot moves the green and the blue bottles to reach the red
bottle. (f) Target is now reachable.

Simulated experiments and the real demonstrations are performed with the
Yaskawa Motoman sda10f, with a robotiq 85 gripper attached on the right arm.

The simulated trials are randomly generated by picking random objects,
dimensions, and collision-free placements within the specified workspace. 20
scenes with each of 6, 8, 10, 12, and 14 objects were used, all of which
contain objects occluded from the camera. Each of the 100 trials is a unique
scene. The target object in each scene was selected to be the hidden object
with the most objects above it, if any.

All tested algorithms were given 20 minutes to run before being terminated. A
trial run is considered successful only if the target object was retrieved
within the time limit. Discovering but failing to pick up the object was still
considered a failure.

Comparisons of the success rate, number of actions for solved trials, total
run-time for solved trials, and number of timed out trials are shown in Figure
4 for 3 algorithms. The algorithms compared are the baseline random action
approach (blue in figure), the proposed resolution complete pipeline without
the object ranking heuristic (orange) and with it (green). For the resolution
complete approaches timing-out is not the only failure mode since they could
detect certain infeasible problems (section V: Failure Detection)

The success rate of the resolution complete approaches is higher than that of
the random baseline. Although its not directly apparent from the plotted
results in Figure 4, the resolution complete approaches (as expected) always
found a solution when the random baseline found a solution; however in one
such trial the resolution complete approach without heuristic exceeded the 20
min time-limit. It is also clear that the heuristic approach has better
success than the non heuristic approach even though they are both complete.
Looking at the data for timed-out experiments, it becomes clear that the
increased success of the heuristic approach is due to timing out less
frequently. This also coincides with the data showing that the heuristic
approach overwhelmingly finds solutions faster and with fewer object
manipulations. In fact, while the non-heuristic RC approach started timing out
linearly with increase in the number of objects, the heuristic approach had
virtually no issue until the scenes got very cluttered with 14 objects.

It is clear that for all methods, success rate starts dropping off
significantly at around 14 objects. This marks the difficulty level for the
given industrial Motoman robot and the workspace. A more compact robot with a
streamlined end-effector (such as the “bluction” tool [7]) could scale to more
cluttered scenes.

Figure 4: On top are shown the graphs of the: success rate (left) and number
of timed out trials (right). On the bottom are the number of actions and the
total runtime for the subset of trials in which all algorithms were
successful.

###  VI-A Integration with Perception & Real Robot Demonstration

The pipeline is directly transferable to scenarios on the real robot to
retrieve a target red bottle from a cluttered shelf. Due to time constraints,
a simple implementation of a perception system is used which only segments and
detects colored cylinders without stacking. Despite the simplifications, a
scene with significant object occlusion is still demonstrated with a
successful retrieval. The proposed pipeline (with heuristic) is run online and
communicates with the robot controller and the RGBD camera for execution and
sensing. The camera extrinsic matrix is estimated by a classical robot-camera
calibration procedure using ArUco markers [23]. For object recognition, the
perception component is implemented via plane fitting, DBSCAN segmentation
[24], and cylinder fitting using Open3D [25]. The plane fitting algorithm
extracts the boundaries of the workspace, which is used to construct the
collision geometries in MoveIt [26]. The inliers of each segmented cylinder
are used to produce the segmentation mask for the RGBD image, which is used to
label occlusion correspondence for each object. To ensure safety, additional
cubic collision geometries are added to the planning scene to avoid collisions
between the robot and the camera. 111Videos can be found at
https://sites.google.com/scarletmail.rutgers.edu/occluded-obj-retrieval.
Extensive experiments of the proposed pipeline were not performed on the real
robot but the demonstration presented was performed a few times and the
pipeline was observed to have qualitatively similar performance as in
simulated experiments; however, calibration and perception issues were
observed to lead to pipeline failure.

##  VII DISCUSSION

It’s worth mentioning that the physical execution accounts for over 60% of
time used for the trials. This shows that there could be room for performance
improvement by performing scene perception asynchronously, since a lot can
still be sensed while the robot is moving. Further performance improvement can
be found by parallelizing the planning of picks and placements for multiple
objects as well.

While this work applies heuristics for selecting objects based on the
occlusion volume, additional information regarding effective placements can
also improve practical performance. In order to solve a larger variety of
problems it would be useful to adapt the placement primitive to allow placing
objects on top of others when there is limited space on the workspace surface.

Another direction is to integrate the task planner with human instructions.
For instance, it would be helpful to use human language to identify the target
as well as influence the search at some regions over others. Additional
heuristics can also be obtained from semantic reasoning of the scene when
objects of the same category tend to be placed closer [27]. Since current
experiments only include simple geometries, such as cylinders and rectangular
prisms, future work can investigate more complex objects where state-of-the-
art perception algorithms are necessary. This would also be necessary for
realistic human-robot integration.

## References

  * [1] L. L. Wong, L. P. Kaelbling, and T. Lozano-Pérez, “Manipulation-based active search for occluded objects,” in 2013 IEEE International Conference on Robotics and Automation, pp. 2814–2819, IEEE, 2013.
  * [2] C. R. Garrett, C. Paxton, T. Lozano-Pérez, L. P. Kaelbling, and D. Fox, “Online replanning in belief space for partially observable task and motion problems,” in 2020 IEEE International Conference on Robotics and Automation (ICRA), pp. 5678–5684, IEEE, 2020.
  * [3] M. R. Dogar, M. C. Koval, A. Tallavajhula, and S. S. Srinivasa, “Object search by manipulation,” Autonomous Robots, vol. 36, no. 1, pp. 153–167, 2014\.
  * [4] C. Nam, S. H. Cheong, J. Lee, D. H. Kim, and C. Kim, “Fast and resilient manipulation planning for object retrieval in cluttered and confined environments,” IEEE Transactions on Robotics, vol. 37, no. 5, pp. 1539–1552, 2021.
  * [5] J. Ahn, J. Lee, S. H. Cheong, C. Kim, and C. Nam, “An integrated approach for determining objects to be relocated and their goal positions inside clutter for object retrieval,” in 2021 IEEE International Conference on Robotics and Automation (ICRA), pp. 6408–6414, IEEE, 2021.
  * [6] W. Bejjani, W. C. Agboh, M. R. Dogar, and M. Leonetti, “Occlusion-aware search for object retrieval in clutter,” in 2021 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), pp. 4678–4685, IEEE, 2021\.
  * [7] H. Huang, M. Danielczuk, C. M. Kim, L. Fu, Z. Tam, J. Ichnowski, A. Angelova, B. Ichter, and K. Goldberg, “Mechanical search on shelves using a novel ”bluction” tool,” IEEE International Conference on Robotics and Automation (ICRA), 2022.
  * [8] A. Price, L. Jin, and D. Berenson, “Inferring occluded geometry improves performance when retrieving an object from dense clutter,” in The International Symposium of Robotics Research, pp. 376–392, Springer, 2019.
  * [9] H. Huang, L. Fu, M. Danielczuk, C. M. Kim, Z. Tam, J. Ichnowski, A. Angelova, B. Ichter, and K. Goldberg, “Mechanical search on shelves with efficient stacking and destacking of objects,” arXiv preprint arXiv:2207.02347, 2022\.
  * [10] Y. Zhu, J. Tremblay, S. Birchfield, and Y. Zhu, “Hierarchical planning for long-horizon manipulation with geometric and symbolic scene graphs,” in 2021 IEEE International Conference on Robotics and Automation (ICRA), pp. 6541–6548, IEEE, 2021.
  * [11] K. N. Kumar, I. Essa, and S. Ha, “Graph-based cluttered scene generation and interactive exploration using deep reinforcement learning,” in 2022 International Conference on Robotics and Automation (ICRA), pp. 7521–7527, IEEE, 2022.
  * [12] Z. Liu, W. Liu, Y. Qin, F. Xiang, M. Gou, S. Xin, M. A. Roa, B. Calli, H. Su, Y. Sun, et al., “Ocrtoc: A cloud-based competition and benchmark for robotic grasping and manipulation,” IEEE Robotics and Automation Letters, vol. 7, no. 1, pp. 486–493, 2021.
  * [13] R. Wang, K. Gao, D. Nakhimovich, J. Yu, and K. E. Bekris, “Uniform Object Rearrangement: From Complete Monotone Primitives to Efficient Non-Monotone Informed Search,” in 2021 IEEE International Conference on Robotics and Automation (ICRA), pp. 6621–6627, May 2021\.  ISSN: 2577-087X.
  * [14] W. Zhao and W. Chen, “Hierarchical pomdp planning for object manipulation in clutter,” Robotics and Autonomous Systems, vol. 139, p. 103736, 2021.
  * [15] Y. Xiao, S. Katt, A. ten Pas, S. Chen, and C. Amato, “Online planning for target object search in clutter under partial observability,” in 2019 International Conference on Robotics and Automation (ICRA), pp. 8241–8247, IEEE, 2019.
  * [16] M. Danielczuk, A. Kurenkov, A. Balakrishna, M. Matl, D. Wang, R. Martín-Martín, A. Garg, S. Savarese, and K. Goldberg, “Mechanical search: Multi-step retrieval of a target object occluded by clutter,” in 2019 International Conference on Robotics and Automation (ICRA), pp. 1614–1621, IEEE, 2019.
  * [17] H. Huang, M. Dominguez-Kuhne, V. Satish, M. Danielczuk, K. Sanders, J. Ichnowski, A. Lee, A. Angelova, V. Vanhoucke, and K. Goldberg, “Mechanical search on shelves using lateral access x-ray,” in 2021 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), pp. 2045–2052, IEEE, 2021.
  * [18] A. Zeng, S. Song, K.-T. Yu, E. Donlon, F. R. Hogan, M. Bauza, D. Ma, O. Taylor, M. Liu, E. Romo, et al., “Robotic pick-and-place of novel objects in clutter with multi-affordance grasping and cross-domain image matching,” The International Journal of Robotics Research, vol. 41, no. 7, pp. 690–705, 2022.
  * [19] Y. Miao, R. Wang, and K. Bekris, “Safe, occlusion-aware manipulation for online object reconstruction in confined spaces,” in International Symposium on Robotics Research (ISRR), 2022.
  * [20] M. Gupta, T. Rühr, M. Beetz, and G. S. Sukhatme, “Interactive environment exploration in clutter,” in 2013 IEEE/RSJ International Conference on Intelligent Robots and Systems, pp. 5265–5272, IEEE, 2013.
  * [21] J. Poon, Y. Cui, J. Ooga, A. Ogawa, and T. Matsubara, “Probabilistic active filtering for object search in clutter,” in 2019 International Conference on Robotics and Automation (ICRA), pp. 7256–7261, IEEE, 2019.
  * [22] T. Novkovic, R. Pautrat, F. Furrer, M. Breyer, R. Siegwart, and J. Nieto, “Object finding in cluttered scenes using interactive perception,” in 2020 IEEE International Conference on Robotics and Automation (ICRA), pp. 8338–8344, IEEE, 2020.
  * [23] R. Y. Tsai, R. K. Lenz, et al., “A new technique for fully autonomous and efficient 3 d robotics hand/eye calibration,” IEEE Transactions on robotics and automation, vol. 5, no. 3, pp. 345–358, 1989.
  * [24] M. Ester, H.-P. Kriegel, J. Sander, X. Xu, et al., “A density-based algorithm for discovering clusters in large spatial databases with noise.,” in kdd, vol. 96, pp. 226–231, 1996.
  * [25] Q.-Y. Zhou, J. Park, and V. Koltun, “Open3d: A modern library for 3d data processing,” arXiv preprint arXiv:1801.09847, 2018.
  * [26] S. Chitta, “Moveit!: an introduction,” in Robot Operating System (ROS), pp. 3–27, Springer, 2016.
  * [27] J. K. Li, D. Hsu, and W. S. Lee, “Act to see and see to act: Pomdp planning for objects search in clutter,” in 2016 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), pp. 5701–5707, IEEE, 2016\.