# Graph-based Detection of Multiuser
Impulse Radio Systems

Yehuda Shen Bahar and Ofer Amrani
School of Electrical Engineering,
Dept. Elec. Eng. - systems,
Tel-Aviv University, Tel-Aviv, Israel
<EMAIL_ADDRESS>+972-3-6407766

Abstract \- Impulse-Radio (IR) is a wideband modulation technique that can
support multiple users by employing random Time-Hopping (TH) combined with
repeated transmissions. The latter is aimed at alleviating the impact of
collisions. This work employs a graphical model for describing the multiuser
system which, in turn, facilitates the inclusion of general coding schemes.
Based on factor graph representation of the system, several iterative
multiuser detectors are presented. These detectors are applicable for any
binary linear coding scheme. The performance of the proposed multiuser
detectors is evaluated via simulations revealing large gains with low
complexity.

Keywords -

##  I Introduction

Impulse-radio systems use short baseband pulses. In Time-Hopping Impulse Radio
(TH-IR) systems the information is encoded in the polarity of the pulses
(Pulse amplitude modulation - PAM), or in the position of the pulses (Pulse
position modulation - PPM). To support multiple access, additional delay is
introduced per pulse. The added delay, which is unique for each user, is
(pseudo) random, and is assumed to be known to the receiver. Additionally, in
order to improve performance, each information bit is typically transmitted
several times. As a result of the above measures, the probability of
catastrophic collision (multiple-user interference) is minimized. In-depth
treatment of Impulse Radio communications is given by Win and Scholtz [1].
While impulse radio is typically associated with Radio-Frequency (RF)
communications, Visible Light Communications (VLC) is yet another data
communication technique that can employ short (light) pulses for signalling,
and can hence benefit from the algorithms and results presented herein.

TH-IR was analyzed in the past for several channel models and interference
types, see e.g. [2],[3], [4]. Scholtz [5] suggested employing this
communication scheme for supporting multiple access. The Multiuser
interference (MUI) was modeled as Gaussian noise and the receiver employed a
matched filter to detect a specific user [5], [1]. Two different receiver
architectures were defined: AIRMA (Analog impulse radio multiple access) and
DIRMA (Digital IRMA) receivers [6]. In another work [7] the choice of a
different time hopping sequence was discussed, a pseudo-chaotic sequence,
aimed at alleviating the MUI problem. Deterministic sequences designed to
mitigate MUI altogether were also proposed [8],[9]. To improve system
performance and better deal with MUI in the framework of single-user
detection, several authors proposed more general coding and modulation
schemes, see e.g. [10], [11],[12],[13].

Rather than treating the signals received from the many users as interference,
one can benefit by performing MultiUser Detection (MUD), thus extracting the
data associated with all users. In general, MUD is a computationally intensive
task. To alleviate this problem, Fishler and Poor [14] suggested using an
iterative multiuser detector for the simple repetition-based code in order to
achieve good performance with low-complexity. Wang et. al. [15] suggested
using the same iterative detector as [14] with different message passing aimed
at reducing computational complexity. Chen et. al. [16] used the conventional
TH-IR scheme concatenated with low-density parity-check (LDPC) coding for
improving system performance. Their receiver employed the standard TH-IR
detector followed by an LDPC decoder. Sathish et. al. [17] proposed using the
standard TH-IR scheme concatenated with convolutional coding. The receiver
employed an iterative soft-input soft-output, 3-stage multiuser detector. The
aforementioned contributions relating to MUD employ concatenation of the
standard TH-IR scheme with coding. It would be interesting to introduce
advanced coding [18] and multiuser detection techniques for further improving
system performance.

LDPC codes were originally introduced by Gallager [19, 20] in 1962, and
"reintroduced" in recent years. It was the introduction of Turbo convolutional
Codes with efficient iterative decoding (which exhibit notably low error
probabilities at low SNR) that triggered the search for such codes - including
LDPC codes on this family of codes (see e.g. [21] [22] [23]).

In this paper we present several iterative multiuser receivers for the
original TH-IR system. Then, we study a modified system where all users employ
arbitrary linear coding, and generalize the above multiuser receiver for this
setting. The proposed receivers are based on factor graph representation of
the complete system. Comparison with the classical repetition-based systems
reveal the significant improvements possible with the proposed practical
approach. The rest of this manuscript is organized as follows. Section II
reviews the system model. In Section III we present iterative multiuser
detectors. Finally, Section IV provides simulation results for the detectors
introduced in the previous section using several codes including LDPC codes.
Some of the results reported herein appeared in [24].

##  II System Model

A typical UWB TH-IR signal can be described by the following general model:

$s_{tr}^{(k)}(t)=\sum_{j=-\infty}^{\infty}b^{k}_{\left\lfloor
j/N_{f}\right\rfloor}w_{tr}(t-jT_{f}-c_{j}^{k}T_{c}),$ (1)

where $S_{tr}^{(k)}$ is the transmitted signal of the k’th user; $T_{f}$ is
the nominal pulse repetition time; $w_{tr}$ is the transmitted pulse shape;
$b_{i}^{k}$ is the $i$th symbol transmitted by user $k$; $c_{j}^{k}$ is the
time hopping sequence used by user $k$; $N_{f}$ is the number of frames in
which a symbol is transmitted; $T_{c}$ is the chip size.

A user repeats every information symbol in $N_{f}$ different frames, where
each frame consists of $N_{c}$ chips, also referred to as slots. The time
hopping sequence employed by each user is a set of values chosen randomly from
$\\{0,1,\ldots,N_{c}-1\\}$. Usually $N_{c}<T_{f}/T_{c}$ for avoiding inter-
symbol-interference (ISI). In this work the information symbols are assumed to
be binary digits (bits) and the signaling is binary-phase shift keying (BPSK).

The system consists of $K$ transmitting users and one receiver, where the
different users are centrally coordinated and synchronized by the receiver
(base station) [25], [26], [27].

The channel over which the signal is transmitted can be frequency selective;
it is assumed that the combined response of the channel and pulse shape is
such that the inter-symbol interference is negligible, and the channel
characteristics are slowly varying in time. The received signal is perturbed
by additive white Gaussian noise (AWGN).

The received signal is given by

$r(t)=\sum_{k=1}^{K}A_{k}\sum_{j=-\infty}^{\infty}b^{k}_{\left\lfloor
j/N_{f}\right\rfloor}w_{rx}(t-jT_{f}-c_{j}^{k}T_{c})+n\left(t\right),$ (2)

where $K$ denotes the number of users in the system; $w_{rx}$ is the received
waveform associated with one transmitted pulse; $A_{k}$ is the amplitude
associated with user $k$; $n\left(t\right)$ denotes an additive white Gaussian
noise process.

In particular, this channel model can be encountered in scenarios where the
dominant propagation path is the line-of-sight, and $T_{c}$ is chosen to
satisfy $T_{c}>Sup\\{w_{rx}\\}$, where $Sup\\{.\\}$ is the support of the
received pulse $w_{rx}$. See also [14] and the references therein.

###  II-A Discrete Time model

The receiver tracks the gains associated with the different users. The
received signal (2) goes through a matched filter whose output is sampled
every $T_{c}$ seconds. Denote by $\mathbf{r}[i]$ the vector of samples at the
output of the matched filter corresponding to the _i_ th information symbol.
The size of this vector is $N_{f}N_{c}$ and it can be described by the
following equation

$\mathbf{r}[i]=\mathbf{S}[i]\mathbf{Ab}[i]+\mathbf{n}[i],$ (3)

where $\mathbf{S}[i]$ is a $[N_{c}N_{f}\times{K}]$ matrix describing the slots
used by the different users for transmitting their _i_ ’th information symbol.
The matrix is defined by

$\mathbf{S}[i]_{lk}=\left\\{\begin{array}[]{ll}1&\textrm{if
$c^{k}_{(i-1)N_{f}+\lfloor\frac{l}{N_{c}}\rfloor
N_{c}}=l-\lfloor\frac{l}{N_{c}}\rfloor N_{c}$}\\\
0&\textrm{otherwise}.\end{array}\right.$ (4)

$\mathbf{A}$ is the gain between the _k_ ’th transmitter and the receiver,
$\mathbf{A}=diag(A_{1},A_{2},\ldots,A_{k})$; $\mathbf{b}[i]$ is the
information vector transmitted by the $K$ users at the _i_ ’th interval,
$\mathbf{b}[i]=[b_{i}^{1},b_{i}^{2},\ldots,b_{i}^{K}]^{T}$; $\mathbf{n}$ is a
zero-mean Gaussian random vector with correlation matrix $\sigma_{n}^{2}I$,
where $\sigma_{n}^{2}=\frac{N_{0}}{2}$.

Without loss of generality we can examine the first information symbol and
therefore omit the time index $i$, consequently (3) can be rewritten as

$\mathbf{r}=\mathbf{SAb}+\mathbf{n}.$ (5)

Figure 1 depicts a toy example of transmitting one information symbol in a
2-user system with the following parameters: $N_{f}=3$ frames, $N_{c}=5$ chips
per frame. The rectangles represent the slots used by the users.

Figure 1: System illustration

##  III Iterative multiuser detection for TH-IR systems

A major drawback of optimum Multi User Detection (MUD) is implementation
complexity. In this section we try to alleviate this problem by presenting
several iterative detectors of reduced complexity.

We first introduce a simple iterative detector for the original TH-IR system
settings presented in Section II. Then, we propose a graph-based description
of the same system, and develop the corresponding MAP-based iterative
detector. Finally, we introduce coding into the original system and develop
the appropriate graph-based iterative multi-user detector.

###  III-A ID detector

The first detector we present is a (very) low complexity iterative detector
based on intuition, rather than mathematical arguments. It is therefore
referred to as ID (Intuition-driven) detector. It is tailored for detecting
repetition-based transmissions, and may be characterized as a Gallager-type
decoder.

Returning to the toy example of Figure 1, the corresponding $\mathbf{S}$
matrix (4) is given by

$\mathbf{S}^{T}=\left[\begin{array}[]{ll}0\ 1\ 0\ 0\ 0\ \ 1\ 0\ 0\ 0\ 0\ \ 0\
0\ 0\ 1\ 0\\\ 0\ 0\ 1\ 0\ 0\ \ 1\ 0\ 0\ 0\ 0\ \ 1\ 0\ 0\ 0\ 0\\\
\end{array}\right],$ (6)

and the graph we associate with the ID detector is given in Figure 2.

Figure 2: ID decoder graph: $K=2$ users, $N_{f}=3$ frames, $N_{c}=5$ chips per
frame.

The detector consists of two stages performed in an iterative manner. In the
first stage, the estimated value of a specific pulse is calculated given the
received signal and the information transmitted by other colliding users.
Associated with the first stage of the decoding is the left-hand side of the
graph consisting of at most $N_{f}N_{c}$ nodes, one for each of the outputs
$r$ of the matched filter. We shall hence refer to these nodes as _input
nodes_.

In the second stage, the estimated value of a bit associated with a specific
user is calculated given the information from the first stage. The second
stage uses the fact that all the pulses of the _k_ th user correspond to the
same bit. The right hand side of the graph is associated with the second stage
of the decoding. It consists of $K$ nodes, referred to as _check nodes_ , one
for each of the $K$ users.

The connection between the two sets of nodes is determined by the matrix
$\mathbf{S}$ \- _input node_ $i$ is connected to _check node_ $j$, if
$S_{ij}=1$. The messages passed between the nodes will are "hard" (binary
valued), rather than "soft" (e.g. LLR).

(a) ID decoder - input node

(b) ID decoder - check node

Figure 3: ID decoder nodes.

We shall now provide a quantitative description of the ID detector beginning
with the input nodes, one of which is depicted in Figure 3(a). Without loss of
generality consider an output message passed via the first edge. Let $m_{0}$
denote the output message, while $m_{j},j\neq 0$, denote all the input
messages, and $r_{i}$ is the output of the matched filter connected to input
node $i$. Note that an input node is connected to $k$ check nodes, $0\leq
k\leq K$ being the number of users choosing to transmit in chip position _i_

$k=\sum_{l=1}^{K}S_{i,l}.$ (7)

The output message passed to the check node is defined as

$m_{0}=r_{i}-\sum_{l=1}^{k-1}A_{i,l}*m_{l},$ (8)

where $A_{i,l}$ is the amplitude of user _l_ as seen by input node $i$.

With the aid of Figure 3(b), we proceed to describe the check nodes. Any check
node is connected to exactly $N_{f}$ input nodes (due to the repetition nature
of the transmission). The message computed by a check node (to be passed back
to an input node) is defined as

$m_{0}=sgn(\sum_{l=1}^{N_{f}-1}m_{l}).$ (9)

Note that for the first iteration, the output from all check nodes is
initialized to zero. Finally, after the last iteration, the output of the
detector is taken from the check nodes, yet it includes all the inputs as
follows

$\hat{b_{k}}=sgn(\sum_{l=0}^{N_{f}-1}m_{l}).$ (10)

###  III-B 3-stage Factor-Graph (FG3) detector

In this section we present an iterative (soft) detector based on a 3-stage
factor-graph description of the system. While the proposed detector is aimed
at decoding the repetition-based transmission, it is of great interest as it
lays the ground for introducing arbitrary linear coding into the system.

Henceforth, we let $\mathbf{y}$ represent the received vector $\mathbf{r}$,
i.e. $\mathbf{y}=\mathbf{r}$. The output of the MAP decoder (for the $k$th
bit) is given by

$\begin{split}\hat{b}^{k}&=\operatornamewithlimits{argmax}_{b^{k}=\pm
1}\\{p(b^{k}|\mathbf{y})\\}\\\ &=\operatornamewithlimits{argmax}_{b^{k}=\pm
1}\left\\{\sum_{-b^{k}}p(\mathbf{y}|\mathbf{b})*\frac{p(\mathbf{b})}{p(\mathbf{y})}\right\\},\end{split}$
(11)

where $\sum_{-b^{k}}$ denotes summing over all values of the vector
$\mathbf{b}$ excluding those containing $-b^{k}$. Since all input vectors,
$\\{\mathbf{b}\\}$ are equiprobable,

$\begin{split}\hat{b}^{k}&=\operatornamewithlimits{argmax}_{b^{k}=\pm
1}\left\\{\sum_{-b^{k}}p(\mathbf{y}|\mathbf{b})\right\\}\\\
&=\operatornamewithlimits{argmax}_{b^{k}=\pm
1}\left\\{\sum_{-b^{k}}p(\mathbf{y}|b^{1},\ldots,b^{K})\right\\}.\\\
\end{split}$ (12)

Using the following definition

$\mathbb{E}(x_{1},\ldots,x_{i})=\left\\{\begin{array}[]{ll}1&ifx_{1}=\ldots=x_{i}\\\
0&otherwise\end{array},\right.$ (13)

one can write

$\begin{split}\hat{b}^{k}_{j}=\operatornamewithlimits{argmax}_{b^{k}_{j}=\pm
1}\\{\sum_{-b^{k}}&p(\mathbf{y}|b_{0}^{1},\ldots,b_{N_{f}-1}^{1},\ldots,b_{K}^{1},\ldots,b_{N_{f}-1}^{K})\\\
&\cdot\prod_{k=1}^{K}\mathbb{E}(b_{0}^{k},\ldots,b_{N_{f}-1}^{k})\\}.\end{split}$
(14)

Note that the subscript $j$, denoting the frame index, has been added though
one expects all $N_{f}$ estimations $\hat{b}_{j}^{k},j\in\\{0,1,...N_{f}-1\\}$
to produce the same value.

Since the channel is memoryless

$p(\mathbf{y}|b_{0}^{1},\ldots,b_{N_{f}-1}^{1},\ldots,b_{K}^{1},\ldots,b_{N_{f}-1}^{K})=\prod_{i=1}^{N_{f}N_{c}-1}p(y_{i}|b_{yi}),$
(15)

where $b_{yi}$ represents the transmitted bits from all users employing chip
slot $i$. Finally, we have

$\hat{b}^{k}_{j}=\operatornamewithlimits{argmax}_{b^{k}_{j}=\pm
1}\left\\{\sum_{-b^{k}}\prod_{i=1}^{N_{f}N_{c}-1}p(y_{i}|b_{yi})\prod_{k=1}^{K}\mathbb{E}(b_{0}^{k},\ldots,b_{N_{f}-1}^{k})\right\\}.$
(16)

This, so-called, "sum-of-products" function can be calculated using a massage
passing algorithm operating iteratively on a bipartite graph, termed factor
graph [28],[29]. A factor graph typically consists of two types of nodes -
function nodes and variable nodes. In Figure 4(a) we depict a graph
corresponding to Equation (16).

(a) FG3 based detector

(b) CFG3 based detector

Figure 4: Factor graph representation used for multiuser detection -
Repetition code example.

The graph consists of three types of nodes. The squares represent function
nodes, which are divided into two types. The first type of function nodes are
associated with the function $p(y_{i}|b_{yi})$; there are exactly $N_{f}N_{c}$
such nodes. We shall denote these as _P function nodes_. The second type of
function nodes are associated with
$\mathbb{E}(b_{0}^{k},\ldots,b_{N_{f}-1}^{k})$; there are exactly $K$ such
nodes. Denote this type of nodes as _E function nodes_. Nodes indicated by
circles represent the variable nodes for $b_{i}^{k}$; there are exactly
$KN_{f}$ such nodes. Denote these nodes as _b variable nodes_. Note that the
edges connecting the _P function nodes_ with the _b variable nodes_ are
completely defined by the matrix $S$. Finally, as will be shown later on, by
omitting the _b variable nodes_ , the graph reduces to the one shown for the
ID detector.

The messages passed by the algorithm shall be denoted by $\mu(x)$, $x=\pm 1$.
$\mu_{p\rightarrow v}$ represents a message passed from a _P function node_ to
a _b variable node_ , $\mu_{v\rightarrow E}$ represents a message from _b
variable node_ to _E function node_ , while $\mu_{E\rightarrow v}$ is a
message from _E function node_ to _b variable node_. Next, we provide explicit
description of the messages associated with the different types of nodes.

####  III-B1 _P function nodes_

The $\mu_{p\rightarrow v}$ messages, calculated at the _P function nodes_ ,
are

$\begin{split}\mu_{p\rightarrow
v}(+1)&=\sum_{x_{1},\ldots,x_{J}}p(y|+1,x_{1},\ldots,x_{J})\prod_{j=1}^{J}\mu_{j}(x_{j});\\\
\mu_{p\rightarrow
v}(-1)&=\sum_{x_{1},\ldots,x_{J}}p(y|-1,x_{1},\ldots,x_{J})\prod_{j=1}^{J}\mu_{j}(x_{j}),\end{split}$
(17)

where $x_{i}$ is a binary value conveyed by an edge connected to an adjacent
variable node. W.l.o.g. $\mu(\pm 1)$ denotes the output message passed to the
variable node of index $0$, while $\mu_{j}(x_{j})$ denotes an input message
originating at an adjacent _b variable node_. (There are $J$ input edges, with
$j$ being the edge index.) We shall henceforth omit the cumbersome directive
pointers, $p\rightarrow v$ and $v\rightarrow p$, as it will always be easy to
realize the correct flow.

Rather than passing two messages, one may use a single message in the form of
the log likelihood ratio (LLR). Let $r$ be defined as follows

$\begin{split}r&\equiv\frac{\mu(+1)}{\mu(-1)}\\\
&=\frac{\sum_{x_{1},\ldots,x_{J}}p(y|1,x_{1},\ldots,x_{J})\prod_{j=1}^{J}\mu_{j}(x_{j})}{\sum_{x_{1},\ldots,x_{J}}p(y|-1,x_{1},\ldots,x_{J})\prod_{j=1}^{J}\mu_{j}(x_{j})}\\\
&=\frac{\sum_{x_{1},\ldots,x_{J}}p(y|1,x_{1},\ldots,x_{J})\prod_{j=1}^{J}\mu_{j}(x_{j})\frac{1}{\mu_{j}(-1)}}{\sum_{x_{1},\ldots,x_{J}}p(y|-1,x_{1},\ldots,x_{J})\prod_{j=1}^{J}\mu_{j}(x_{j})\frac{1}{\mu_{j}(-1)}}.\end{split}$
(18)

Note that

$\frac{\mu_{j}(x_{j})}{\mu_{j}(-1)}=r_{j}^{\frac{x_{j}+1}{2}}.$ (19)

Finally, using (19) the LLR, $l=log(r)$, is given by

$l=\log{\frac{\sum_{x_{1},\ldots,x_{J}}p(y|1,x_{1},\ldots,x_{J})\prod_{j=1}^{J}{e^{l_{j}}}^{\frac{x_{j}+1}{2}}}{\sum_{x_{1},\ldots,x_{J}}p(y|-1,x_{1},\ldots,x_{J})\prod_{j=1}^{J}{e^{l_{j}}}^{\frac{x_{j}+1}{2}}}}.$
(20)

####  III-B2 _b variable nodes_

It is easy to see that each of the _b variable nodes_ is connected to exactly
one _P function node_ on the left and one _E function node_ on the right. In
this case the output messages are the same as the input messages.

####  III-B3 _E function nodes_

For this type of nodes we have

$\begin{split}r&=\frac{\sum_{x_{1},\ldots,x_{N_{f}-1}}\mathbb{E}(1,x_{1},\ldots,x_{N_{f}-1})\prod_{j=1}^{N_{f}-1}\mu_{j}(x_{j})}{\sum_{x_{1},\ldots,x_{N_{f}-1}}\mathbb{E}(-1,x_{1},\ldots,x_{N_{f}-1})\prod_{j=1}^{N_{f}-1}\mu_{j}(x_{j})}\\\
&=\frac{\sum_{x_{1},\ldots,x_{N_{f}-1}}\mathbb{E}(1,x_{1},\ldots,x_{N_{f}-1})\prod_{j=i}^{N_{f}-1}r_{j}^{\frac{x_{j}+1}{2}}}{\sum_{x_{1},\ldots,x_{N_{f}-1}}\mathbb{E}(-1,x_{1},\ldots,x_{N_{f}-1})\prod_{j=i}^{N_{f}-1}r_{j}^{\frac{x_{j}+1}{2}}.}\end{split}$
(21)

Recalling the definition of $\mathbb{E}$, (13), it is easily verified that

$\begin{split}r&=\frac{\prod_{j=1}^{N_{f}-1}r_{j}^{\frac{1+1}{2}}}{\prod_{j=1}^{N_{f}-1}r_{j}^{\frac{-1+1}{2}}}\\\
&=\prod_{j=1}^{N_{f}-1}r_{j},\end{split}$ (22)

and the LLR is simply

$l=log(r)=\sum_{j=1}^{N_{f}-1}l_{j}.$ (23)

After the last iteration is performed, the output of the detector is taken
from the _E function node_ :

$\hat{b}^{k}=sgn(\sum_{j=1}^{N_{f}}l_{j}).$ (24)

We described a multiuser detector based on factor graph representation of the
system. Recall that Fishler and Poor (FP) [14] presented an iterative
multiuser detector for the same system that follows the turbo principle.
Interestingly, the FG3 detector turns out to be the same as the FP detector
although a different model is employed for describing the system. This
assertion may not be obvious by simply comparing at the massage passing
equations. However, by using simple mathematical manipulations one can move
from the set of equations presented herein to those used by FP. It can be
shown that the FP and FG3 detectors are not MAP-achieving since the associated
graphs are not cyclic-free.

###  III-C CFG3 detector

The FG3 proposed model and detection technique shall now serve as the basis
for the introduction of a coded system. The receiver to be developed in this
section is aimed at providing a solution for a system employing arbitrary
linear coding. It shall hence be denoted as coded-FG3 (CFG3).

A linear code is typically defined by its parameters $[n,k,d]$, where $n$ is
the code length, $k$ is the dimension of the code (not to be confused with the
number of users), and $d$ denotes the minimum Hamming distance of the code.
Hence, in a coded TH-IR system, the number of frames $N_{f}$ will satisfy
$N_{f}=n$, and the system rate is $\frac{k}{n}$.

Referring to Equation (16), the function $\mathbb{E}(\cdot)$ represents the
fact that the bits associated with a specific user must all be the same. When
a code $[n,k,d]$ is used, we shall replace $\mathbb{E}(\cdot)$ with a new
function, $\mathbb{PC}$, representing a Parity Check equation.

The parity check matrix $\mathbf{H}$ is an $(n-k)\times n$ binary matrix

$\mathbb{PC}=\prod_{i=1}^{n-k}PC_{i},$ (25)

where $PC_{i}$ is defined as follows

$\mathbb{PC}_{i}=\left\\{\begin{array}[]{ll}1&\textrm{if
($\sum_{j,h_{i,j}=1}b_{j})=1$ where the sum is over GF2}\\\
0&\textrm{otherwise}\end{array}.\right.$ (26)

As an example, Figure 4(b) depicts a factor graph of a linear code when the
latter is a repetition code whose parity check matrix $\mathbf{H}$ is given by

$\mathbf{H}=\left(\begin{array}[]{ccccc}1&0&0&\ldots&1\\\ 0&1&0&\ldots&1\\\
\vdots&\vdots&\ddots\\\ 0&0&\ldots&1&1\\\ \end{array}\right).$

Although this graph is different from Figure 4(a), which means that two
different detectors are to be used, both graphs target the same system and
code.

The _b variable nodes_ have one input (form the left), but typically several
outputs (connected to PC function nodes). The messages to be calculated and
passed to the output are

$\begin{split}\mu(1)&=\prod_{i}\mu_{i}(1)\\\
\mu(-1)&=\prod_{i}\mu_{i}(-1),\end{split}$ (27)

and hence

$\begin{split}r=\frac{\mu(1)}{\mu(-1)}=\prod_{i}r_{i}\\\
l=log(r)=\sum_{i}l_{i}.\end{split}$ (28)

The messages to be calculated at the _PC function nodes_ are associated with
single parity check equations. These are the same messages used in LDPC
decoding:

$l=2\tanh^{-1}(\prod\tanh\frac{l_{j}}{2}).$ (29)

After iterating all messages through the graph, the final marginalization is
performed at the _b variable nodes_

$\hat{b}^{k}=sgn(\sum_{j}l_{j}).$ (30)

It can be shown that by appropriately choosing the code, the associated graph
can be made cyclic-free, and hence CFG3 is MAP-achieving.

##  IV Simulation Results

###  IV-A Conventional repetition-based system

We first study the original system employing a repeated transmission scheme.
We shall consider two different detectors: the FG3 detector and the CFG3
detector. In all presented simulations the number of slots per frame is
$N_{c}=20$, while the number of users varies, $K=\\{3,10,30\\}$.

Figure 5: Repetition-based system, FG3 and CFG3 detectors, $Nc=20$

Figures 5 hold the simulation results. The performance of the FG3 and the CFG3
detectors are practically the same. Therefore, whenever using repetition-based
coding, we shall only consider the FG3 detectors.

###  IV-B Repetition vs. LDPC-based systems

We next compare the performance of a repeated transmission system (using ID
and FG3 detectors) with an LDPC-based coded system (using CFG3 detector). The
parameters of the LDPC code chosen are $(n=120,k=56,R=0.4667)$ \- Mackay code
120.64.3.109 [30]. It has been chosen because of its relatively short
length111Recall that the number of frames, $N_{f}$, in our system model is
actually the code length. Since the total number of time slots in the system
is $N_{f}N_{c}$, reasonable simulation times are when using shorter codes.; no
attempt has been made to optimize it for our system. Repetition code of rate
$\frac{1}{3}\rightarrow N_{f}=3$ have been chosen for comparison.

Figures 6 and 7 pertain to the repetition code with ID and FG3 detectors,
respectively, while Figure 8 pertains to the LDPC code-based system employing
the CFG3 detector. The ID and FG3 detectors exhibit an error floor which is
irreducible due to the fact that system performance is limited by Multi User
Interference (MUI). To support this assertion note how the error floor
increases with the number of users. The FG3 detector exhibits better error
floor performance than the ID detector, as expected.

Clearly, the LDPC-based system with the CFG3 detector behaves differently.
First, the BER curve is much sharper as might be expected of a coded system.
Second, as demonstrated in Figure 8, when the number of users increases, the
$E_{b}/N_{0}$ threshold-point also increases. Still, once the threshold is
passed, the slopes associated with all cases are similar. Increasing the
number of users amounts to adding more noise, which leads us to the next
observation: another threshold that can be clearly identified corresponds to
the number of users - beyond a certain number of users, the system collapses.

In general, we argue that the LDPC-based system handles MUI much better than
the alternative approaches mainly because it employs a large number of frames
($N_{f}$). Consequently, "catastrophic" collisions among different users are
much less common.

Complexity-wise we argue as follows. ID performs only simple addition
operations. Calculation of the message at the input node (8) requires $O(K`)$
additions, where $K`$ is the number of colliding users. Check node calculation
(9) requires $O(N_{f})$ additions. Recall that the graphs associated with the
ID and FG3 systems are the same. The two detectors differ only in the
calculations carried out at the input nodes. Calculating (20), instead of (8),
increases complexity by $O(2^{K`})$ [14]. Moving from FG3 to CFG3, as the
input nodes remain the same, the added complexity is that of decoding. In our
example, LDPC coding has been employed, and therefore the additional
complexity is given by $C_{LDPC}\cdot K$, where $C_{LDPC}$ is the complexity
required for decoding the LDPC code used.

Compared to ID, FG3 better treats MUI for the price of increased complexity at
the input side (left-hand side of the graph). CFG3 add coding gain for the
price of increased complexity at the coding side (right-hand side of the
graph). In conclusion, among the three described schemes (for the same
$E_{b}/N_{0}$), the LDPC-based system provides the best performance both in
terms of BER and overall system throughput.

Figure 6: ID detector, Repetition, $N_{f}=3$, $N_{c}=20$ Figure 7: FG3
detector, Repetition, $N_{f}=3$, $N_{c}=20$ Figure 8: CFG3 detector, LDPC,
($n=120$,$k=56$,$R=0.4667$), $N_{c}=20$

###  IV-C Performance dependency on the number of users

In a multi-user environment users often join and leave the system. Analyzing
the performance of the system for different numbers of users is therefore of
interest.

Figure 9: Multi-User Performance, $E_{b}/N_{0}=20dB$ Figure 10: Multi-User
Performance, $E_{b}/N_{0}=5dB$

A system with $N_{c}=20$ slots per frame is considered: ID and FG3 detectors
are employed for a repetition-based transmission scheme with $N_{f}=3$
$(rate=1/3)$; CFG3 detector with LDPC code $(n=120$,$k=56$,$R=0.4667)$ is also
considered.

Figures 9 and 10 present the simulation results for $E_{b}/N_{0}=20dB$ and
$E_{b}/N_{0}=5dB$, respectively. For $E_{b}/N_{0}=20dB$, MUI is the dominant
impairment. The performance of the ID and FG3 detectors gradually deteriorate
as the number of users is increased from 3 to 40. FG3 exhibits better behavior
throughout. The LDPC coded system with CFG3 detector performs quite
differently. No errors are observed until the number of users reaches a
certain threshold, 33 in this simulation. Once the threshold is passed,
performance degradation is steep. Exact position of the threshold and its
behavior depend on system parameters.

In the second simulation, shown in Figure 10, the additive noise is not
negligible. The CFG3 detector behaves much like before: as long as the number
of users is below a certain threshold, performance degradation is quite
moderate. Also notice that the threshold moved to the left due to the
existence of noticeable noise. The ID and FG3 detectors perform much worse
than the CFG3 detector with FG3 being only slightly superior to the ID
detector.

###  IV-D Performance dependency on iterations

The influence of the number of iterations on performance has so far been
overlooked; a fixed number of $8$ iterations was used throughout. Notably,
increasing the number of iterations benefits the LDPC-based system more than
the alternatives.

As an illustrative example Figure 11 presents the BER performance of the three
detectors, as a function of $E_{b}/N_{0}$, for $K=20$ users and different
number of iterations.

Figure 11: Performance dependency on number of iterations, $K=20$ users

The performance of the ID and FG3 detectors converge in only 4 iterations.
More iterations are required for the LDPC-based system to converge. In the
single-user case (simulation results not shown), the number of iterations
required for the LDPC-based system to converge is smaller than in the multi-
user case. One may argue that this follows from the fact that the graph
associated with the multi-user case is much more involved.

###  IV-E Some after simulation comments

In a single-user system, the user can clearly employ all the time slots for
its own usage. In the case of high SNR, with two-level modulation, the channel
capacity is upper bounded by $1[bit/chip]$. Since each frame provides $N_{c}$
transmit opportunities, the capacity of the single-user system is upper
bounded by $C=N_{c}[bits/frame]$. Let us now consider the CFG3-based system
with $K$ users. The code rate used was $0.46$ and $N_{c}=20$. The CFG3
detector hardly produced any bit errors as long as the number of users
satisfied $K<33$. The overall system throughput in this case is
$C=14.72[bits/frame]$ as compared to the above upper bound of
$C=N_{c}=20[bits/frame]$. With an optimized LDPC code of greater length, the
obtained throughput is expected to grow closer to this bound.

##  V Conclusions

This work is concerned with iterative receivers for a multiuser TH-IR system.
In particular, we study the gain achieved by introducing coding into the
originally uncoded system. Several iterative multi-user detectors have been
presented based on factor graph representation of the complete system. These
detectors are general and can support any binary linear coding scheme. Yet
another strength of the proposed approach is that the graph used may be
extended to account for multipath components (owing to the nature of the UWB
channel). This is an interesting topic for future work.

Simulation results for the above mentioned detectors are presented for several
codes including LDPC codes. It is demonstrated that for $20$ users, BER of
$10^{-3}$ can be achieved with higher system rate and more than $4dB$ gain in
$E_{b}/N_{0}$ when using an LDPC-coded system. Furthermore, it is shown that
the achievable performance of the original system quickly saturates at
relatively poor BER.

## References

  * [1] M.Z. Win and R.A. Scholtz.  Impulse radio: how it works.  Communications Letters, IEEE, 2(2):36–38, Feb 1998.
  * [2] I. Bergel, E. Fishler, and H. Messer.  Narrowband interference suppression in time-hopping impulse-radio systems.  In Ultra Wideband Systems and Technologies. Digest of Papers. 2002 IEEE Conference on, pages 303–307, 2002.
  * [3] Li Zhao and A.M. Haimovich.  Performance of ultra-wideband communications in the presence of interference.  Selected Areas in Communications, IEEE Journal on, 20(9):1684–1691, Dec 2002.
  * [4] F. Ramirez-Mireles.  On the performance of ultra-wide-band signals in gaussian noise and dense multipath.  Vehicular Technology, IEEE Transactions on, 50(1):244–249, Jan 2001\.
  * [5] R. Scholtz.  Multiple access with time-hopping impulse modulation.  In Military Communications Conference, 1993. MILCOM ’93. Conference record. ’Communications on the Move’., IEEE, volume 2, pages 447–450 vol.2, Oct 1993.
  * [6] M.Z. Win and R.A. Scholtz.  Ultra-wide bandwidth time-hopping spread-spectrum impulse radio for wireless multiple-access communications.  Communications, IEEE Transactions on, 48(4):679–689, Apr 2000.
  * [7] D.C. Laney, G.M. Maggio, F. Lehmann, and L. Larson.  Multiple access for uwb impulse radio with pseudochaotic time hopping.  Selected Areas in Communications, IEEE Journal on, 20(9):1692–1700, Dec 2002.
  * [8] C.J. Le Martret and G.B. Giannakis.  All-digital impulse radio with multiuser detection for wireless cellular systems.  Communications, IEEE Transactions on, 50(9):1440–1450, Sep 2002\.
  * [9] Liuqing Yang and G.B. Giannakis.  Block-spreading codes for impulse radio multiple access through isi channels.  In Communications . ICC 2002. IEEE International Conference on, volume 2, pages 807–811 vol.2, 2002.
  * [10] Xiaoli Chu and R.D. Murch.  Multidimensional modulation for ultra-wideband multiple-access impulse radio in wireless multipath channels.  Wireless Communications, IEEE Transactions on, 4(5):2373–2386, Sept. 2005.
  * [11] Yafei Tian and Chenyang Yang.  Coded modulation in uwb communications.  In Personal, Indoor and Mobile Radio Communications, 2005. PIMRC 2005\. IEEE 16th International Symposium on, volume 1, pages 491–495, Sept. 2005\.
  * [12] Youn Seok Kim, Won Mee Jang, and Lim Nguyen.  Self-encoded th-ppm uwb system with iterative detection.  In Advanced Communication Technology, 2006. ICACT 2006. The 8th International Conference, volume 1, pages 5 pp.–714, Feb. 2006.
  * [13] H. Zhang, W. Li, and T.A. Gulliver.  Pulse position amplitude modulation for time-hopping multiple-access uwb communications.  Communications, IEEE Transactions on, 53(8):1269–1273, Aug. 2005\.
  * [14] E. Fishler and H.V. Poor.  Low-complexity multiuser detectors for time-hopping impulse-radio systems.  Signal Processing, IEEE Transactions on, 52(9):2561–2571, Sept. 2004.
  * [15] Bao-Yun Wang and Wei Xing Zheng.  Multiuser detection for time-hopping impulse-radio systems using a low-complextiy iterative algorithm.
  * [16] Fangni Chen, Shiju Li, and Lei Shen.  Performance evaluation of a novel low-density parity-check coded ultra wide band impulse radio (uwb-ir) system.  In Wireless Communications, Networking and Mobile Computing, 2006\. WiCOM 2006.International Conference on, pages 1–4, Sept. 2006.
  * [17] A.D. Sathish and S.K. Jayaweera.  Iterative low-complexity multiuser detection and decoding for coded uwb systems.  In Spread Spectrum Techniques and Applications, 2006 IEEE Ninth International Symposium on, pages 173–177, Aug. 2006.
  * [18] Chia-Chin Chong, F. Watanabe, and H. Inamura.  Potential of uwb technology for the next generation wireless communications.  In Spread Spectrum Techniques and Applications, 2006 IEEE Ninth International Symposium on, pages 422–429, Aug. 2006.
  * [19] R. Gallager.  Low-density parity-check codes.  Information Theory, IRE Transactions on, 8(1):21–28, January 1962\.
  * [20] R. Gallager.  Low Density Parity Check Codes.  M.I.T. Press, Cambridge, Massachusetts, 1963.
  * [21] D.J.C. MacKay.  Good error-correcting codes based on very sparse matrices.  Information Theory, IEEE Transactions on, 45(2):399–431, March 1999\.
  * [22] T.J. Richardson and R.L. Urbanke.  The capacity of low-density parity-check codes under message-passing decoding.  Information Theory, IEEE Transactions on, 47(2):599–618, February 2001.
  * [23] T.J. Richardson, M.A. Shokrollahi, and R.L. Urbanke.  Design of capacity-approaching irregular low-density parity-check codes.  Information Theory, IEEE Transactions on, 47(2):619–637, February 2001.
  * [24] Y.S. Bahar and O. Amrani.  Iterative multiuser detection in coded uwb.  In Electrical and Electronics Engineers in Israel. IEEEI 2008. IEEE 25th Convention of, pages 450–451, Dec. 2008.
  * [25] Chee-Cheon Chui and R.A. Scholtz.  A synchronizing scheme for an impulse network.  In Military Communications Conference, MILCOM 2004. IEEE, volume 2, pages 635–641 Vol. 2, Oct.-3 Nov. 2004.
  * [26] C. Carbonelli, S. Franz, U. Mengali, and U. Mitra.  Semi-blind ml synchronization for uwb transmitted reference systems.  In Signals, Systems and Computers, 2004. Conference Record of the Thirty-Eighth Asilomar Conference on, volume 2, pages 1491–1495 Vol.2, Nov. 2004.
  * [27] E.A. Homier and R.A. Scholtz.  Rapid acquisition of ultra-wideband signals in the dense multipath channel.  In Ultra Wideband Systems and Technologies. Digest of Papers. 2002 IEEE Conference on, pages 105–109, 2002.
  * [28] F.R. Kschischang, B.J. Frey, and H.-A. Loeliger.  Factor graphs and the sum-product algorithm.  Information Theory, IEEE Transactions on, 47(2):498–519, Feb 2001\.
  * [29] H.-A. Loeliger.  An introduction to factor graphs.  Signal Processing Magazine, IEEE, 21(1):28–41, Jan. 2004.
  * [30] D.J.C. MacKay.  Encyclopedia of Sparse Graph Codes.  http://www.inference.phy.cam.ac.uk/mackay/codes/data.html.


# Learning from Demonstrations for Autonomous Soft-tissue Retraction*

Ameya Pore1,2, Eleonora Tagliabue1, Marco Piccinelli1, Diego Dall’Alba1,
Alicia Casals2, Paolo Fiorini1 1 Department of Computer Science, University of
Verona, Italy. Email<EMAIL_ADDRESS>Automatic Control and Computer
Engineering Department, Universitat Politècnica de Catalunya, Barcelona, Spain
*This project has received funding from the European Union’s Horizon 2020
research and innovation programme under the Marie Skłodowska-Curie (grant
agreement No. 813782 ”ATLAS”) and (grant agreement No. 742671 ”ARS”)

###### Abstract

The current research focus in Robot-Assisted Minimally Invasive Surgery
(RAMIS) is directed towards increasing the level of robot autonomy, to place
surgeons in a supervisory position. Although Learning from Demonstrations
(LfD) approaches are among the preferred ways for an autonomous surgical
system to learn expert gestures, they require a high number of demonstrations
and show poor generalization to the variable conditions of the surgical
environment. In this work, we propose an LfD methodology based on Generative
Adversarial Imitation Learning (GAIL) that is built on a Deep Reinforcement
Learning (DRL) setting. GAIL combines generative adversarial networks to learn
the distribution of expert trajectories with a DRL setting to ensure
generalisation of trajectories providing human-like behaviour. We consider
automation of tissue retraction, a common RAMIS task that involves soft
tissues manipulation to expose a region of interest. In our proposed
methodology, a small set of expert trajectories can be acquired through the da
Vinci Research Kit (dVRK) and used to train the proposed LfD method inside a
simulated environment. Results indicate that our methodology can accomplish
the tissue retraction task with human-like behaviour while being more sample-
efficient than the baseline DRL method. Towards the end, we show that the
learnt policies can be successfully transferred to the real robotic platform
and deployed for soft tissue retraction on a synthetic phantom.

##  I INTRODUCTION

Robot-Assisted Minimally Invasive Surgery (RAMIS) is a consolidated paradigm
in the field of medical robotics. RAMIS is a viable alternative to other
surgical approaches that reduce several patient complications such as
excessive intraoperative blood loss, post-operative trauma, and a high amount
of mortality associated with more invasive surgeries [1]. One of the widely
adopted platforms for RAMIS is the da Vinci Surgical System (Intuitive
Surgical Inc.), a teleoperated system that enables dexterous control with
enhanced precision, stability and accuracy.

Figure 1: The proposed methodology of LfD for the tissue retraction surgical
gesture. (a) Expert demonstrations are performed and recorded using the dVRK
console (b) Robotic agent is trained within a simulated environment. (c) The
learnt policy is translated to the real robotic system.

A significant portion of RAMIS procedures is spent in mobilizing and
manipulating tissues to reach the region of interest. The tissue is
repetitively grasped and retracted to expose the underlying anatomical area
[2]. This gesture of tissue retraction occurs in multiple phases of surgery
and requires extensive interaction with soft tissues having heterogeneous
physical and geometric properties, such as stiffness and viscoelasticity, with
high inter and intra-subject variability. During RAMIS, tissue retraction task
may require the surgeon either to switch robotic arms during surgery with a
different set of visuomotor feedback and limited perception or to instruct an
assistant with the desired motion [3]. This involves risks such as instrument
collision or tissues damage that can negatively affect the procedure.
Therefore, soft tissue retraction is an ideal candidate for the automation of
surgical gestures since it would help to reduce the cognitive and physical
workload of the surgeons and place them in a supervisory position.

Standard motion planning techniques such as potential fields have been
proposed for the automation of surgical gestures [4]. Such methods have shown
to be successful in static conditions, but fail to generalize to dynamic
environments like the surgical one, where human-level dexterity and
adaptability are required. Hence, Learning from Demonstrations (LfD) is a
preferred way to learn human gestures. However, a significant drawback of LfD
methods is that they require a huge number of demonstrations to be trained
properly, which is unfeasible in clinical settings considering the time,
resources and ethical constraints. Moreover, LfD methods are affected by
several limitations, such as errors in the acquisition and poor generalisation
performance [5]. Using LfD, the robot can only become as good as the human’s
demonstrations. There is no additional information for improving the learnt
behaviour.

By contrast, Deep Reinforcement Learning (DRL) allows the robot to discover
new control policies through free exploration of the state-action space. DRL
has shown generalisation capabilities in learning adaptable behaviours for
complex and diverse scenarios such as dexterous manipulation and grasping [6].
However, DRL methods often take a long time to converge and require a well-
shaped and carefully designed reward function to learn simple goal-directed
behaviours. Approaches that combine LfD and DRL aim at exploiting the strength
of both to overcome their respective drawbacks. Particularly, demonstrations
can be used to guide the exploration done during learning, reducing the time
required to find an improved control policy, which may depart from the
demonstrated behaviour.

In this work, we address the limitations of LfD approaches by proposing a
training methodology for surgical robots based on Generative Adversarial
Imitation Learning (GAIL). GAIL relies on generative adversarial networks to
generate trajectories that are similar to the acquired expert demonstrations
[7]. To automate soft tissue retraction, we train GAIL in a simulated replica
of the real environment using a small set of demonstrations acquired through
the real robotic platform via teleoperation. This work represents an initial
step for utilizing expert demonstrations to learn a control policy and to
transfer the motion skills to a real robotic manipulator. Hence, our
contribution is the introduction of an LfD training methodology to learn
human-level surgical gestures using a small set of task demonstrations.

In Sec. II, we provide an overview of the related works. We summarise the
mathematical formulation of the DRL methods used in Sec. III. In Sec. IV, we
detail the training methodology and the learning setup. In Sec. V, we describe
experiments conducted to validate our methodology. Further, we present the
results in Sec. VI along with discussion. In the final Sec. VII, we elaborate
the conclusions and highlight our future works.

##  II RELATED WORKS

Previous works to automate soft-tissue retraction use motion planning
algorithms such as probabilistic roadmaps for optimization based objectives
[2]. This method works well for pre-operative planning in a static
environment. Nagy et al. developed an approach for tissue retraction based on
images, where three methods based on proportional control, Hidden Markov
Models and fuzzy logic are validated [8]. Further, Attanasio et al. developed
a trajectory planner based on coordinates extracted from images [3]. Both of
these studies require hand-crafting control strategies and movement sequences.
The execution of complex non-linear trajectories and behaviours may be
challenging using these methods.

DRL for surgical tasks: Prior works have used DRL to learn the tensioning
policy for tissue manipulation [9]. Pedram et al. developed a Q-learning RL
algorithm based on visual features to learn a control policy for soft-tissue
manipulation. Shin et al. developed an adaptive Model Predictive Controller
(MPC) to learn the tissue dynamics [10]. Human demonstrations are used to
pretrain and initialize the tissue dynamics model. The learnt tensioning
policy is used to develop a visual model-based RL algorithm that outputs
actions to manipulate the tissue points to pre-specified locations. In the
former work, the tissue dynamics is learnt implicitly in the DRL setting,
whereas in the latter, the dynamics is learnt via an MPC. These methods would
be well suited for learning dynamics of deformable objects with similar
physical properties. However, variability in soft tissues properties would
require complete re-training of these methods for each subject. Hence, we base
our hypothesis on learning task features of tissue manipulation from expert
demonstrations in contrast to learning tissue dynamics. Task features would
inherently include the surgeon’s knowledge in manipulating a variety of
tissues and would demonstrate better generalisation.

In the context of evaluation platform, Richter et al. proposed a framework for
training DRL policies in simulation for the dVRK system and showed that the
learnt policy could be transferred to a real robotic system [11]. The proposed
environment does not support deformable objects, which is a major limitation
when simulating the surgical scenario. Recently, Tagliabue et al. proposed a
virtual framework called UnityFlexML for simulating deformable tissues that is
well suited to train DRL methods [12]. We use this simulation framework to
develop our LfD approach.

LfD for surgical tasks: Reiley et al. proposed a demonstration based framework
using Gaussian Mixture Models (GMM) for motion generation [13]. Recently, a
similar approach of GMM has been used to learn dynamic motion primitives from
the demonstrations obtained from expert surgeons [14]. Osa et al. introduced
an iterative technique to learn a single reference trajectory for knot tying
[15]. However, a single demonstration does provide enough consistency to model
a manipulation skill. Schulman et al. used a trajectory transfer algorithm to
learn from demonstrations for the task of suturing [16]. Murali et al. devised
a method to segment demonstrations into motion sequences [17]. Standard
approaches for imitation learning such as Behaviour Cloning (BC) and Inverse
Reinforcement Learning (IRL), tend to be successful with large amounts of
data. However, they suffer from compounding error caused by covariant shift
and are extremely expensive to train [6], hence not suited for the surgical
paradigm. An emerging derivative of IRL is represented by GAIL [7]. Contrary
to IRL, which learns a cost function, GAIL methods learn how to act by
directly learning the policy. GAIL has shown an empirical improvement in
reducing the number of demonstrations required to successfully learn the task,
compared to other imitation learning methods. Therefore we selected GAIL as
the LfD method, since it has been successfully applied in other surgical
domains (e.g., endovascular manipulators [18]) but never used in the RAMIS
context.

##  III BACKGROUND

A general Reinforcement Learning (RL) problem can be formulated in terms of a
Markov Decision Process (MDP), where an agent learns by interacting with the
environment. An MDP $\mathcal{M}$ is defined as a tuple
$(\mathcal{S},\mathcal{A},r,\mathcal{P},\gamma,T)$, where $\mathcal{S}$ is a
set of possible states, $\mathcal{A}$ is the set of actions, $\mathcal{P}$ is
the transition probability distribution, $r$ is the reward function,
$\gamma\in[0,1]$ is the discount factor and T is the time horizon per episode.
In RL, at each timestep $t$, the environment produces a state observation
$s_{t}\in\mathcal{S}$. The agent then samples an action $a_{t}\sim\pi(s_{t})$,
$a_{t}\in\mathcal{A}$ and applies the action to the environment. As a
consequence, the agent transitions to a new state $s_{t+1}$ sampled from the
transition function $p(s_{t+1}|s_{t},a_{t})$, $p\in\mathcal{P}$ or terminates
the episode at state $s_{T}$. The agent’s goal is to learn a stochastic
behaviour policy $\pi$ parametrised by $\phi$,
$\pi_{\phi}:\mathcal{S}\rightarrow\mathcal{P}(\mathcal{A})$ to maximise the
expected future discounted reward $E[\sum_{i=0}^{T-1}\gamma^{i}r_{i}]$.

###  III-A Proximal Policy Optimisation (PPO)

PPO is an on-policy RL algorithm capable of dealing with both continuous and
discrete action spaces. PPO alternates between collecting new observations and
improving the policy, while approximating the value function as well [19]. The
update function for the PPO policy is the following

$L(s_{t},a_{t},\theta_{k},\theta)=min\left(\frac{\pi_{\theta}(a_{t}|s_{t})}{\pi_{\theta_{k}}(a_{t}|s_{t})}\hat{A}^{\pi_{\theta_{k}}}(s_{t},a_{t}),g(\epsilon,\hat{A}^{\pi_{\theta_{k}}}(s_{t},a_{t}))\right)$
(1)

where $\theta_{k}$ are the parameters of the old policy, and $g$ is defined
as:

$g(\epsilon,\hat{A})=\begin{cases}(1+\epsilon)\hat{A}_{t},&\hat{A}_{t}\geq
0\\\ (1-\epsilon)\hat{A}_{t},&\hat{A}_{t}<0\end{cases}$ (2)

Where $\hat{A}_{t}$ is the advantage estimator function at timestep $t$ and
$\epsilon$ is a hyperparameter. The idea behind PPO is to limit the impact of
the policy update using the min operator, so that the improvement in the
policy is stable.

###  III-B Generative Adversarial Imitation Learning (GAIL)

GAIL is an imitation learning algorithm [7]. The principle of GAIL is based on
the generative adversarial networks, which consists of a discriminator
$D_{\phi^{\prime}}$ and a policy generator $G_{\phi}$, where $\phi$ denotes
the parameters associated with each network. The policy network generates
exploration trajectories which are used by the discriminator to compute a
surrogate function measuring the similarity between the generated policy and
the expert policy. This similarity metric acts as reward proxy for the RL
step. Unlike IRL techniques, GAIL directly generates policies instead of the
reward function. The discriminator is trained to minimize the loss function:

$L_{GAIL}=E_{\tau_{\phi}}[\log(D_{\phi^{\prime}}(s_{t},a))]+E_{\tau_{E}}[\log(1-D_{\phi^{\prime}}(s_{t},a))]$
(3)

where $\tau_{\phi}$ are the trajectories generated by $G_{\phi}$ and
$\tau_{E}$ are the expert trajectories. The policy generator $G_{\phi}$ is
often adopted from methods based on stochastic policy such as PPO. There are
two reasons why PPO is used for GAIL: first, PPO uses smooth policy update for
stable learning and second, PPO generates diversified trajectories that act as
a wide sampling range for the discriminator in GAIL.

##  IV METHODS

The objective of our task is to efficiently and optimally train the robotic
agent to expose the tumour region by grasping and retracting the fat tissue
surrounding it. In particular, our DRL agent is represented by the End-
Effector (EE) of the da Vinci Patient Side Manipulator (PSM) arm, which learns
to move from an initial position $\mathbf{p}_{0}$ to a position close to the
tumour $\mathbf{q}$, grasp the tissue and retract it to reach the desired
target position $\mathbf{p}_{T}$. For the sake of simplicity, the PSM
orientation is kept constant. The initial state of the anatomical environment,
as well as the tumour centroid position $\mathbf{q}$, are assumed to be known
from pre-operative data. In order to make the learnt motion primitives robust
to different initial configurations, the EE starts from a different position
after each episode (i.e. 2500 timesteps in our case) at training time. The
considered state-space leverages solely on kinematics information defining the
current robot state and the environment:

$S_{t}=[\mathbf{p}_{t},\mathbf{q},\mathbf{p}_{T},||\mathbf{p}_{t}-\mathbf{q}||,||\mathbf{p}_{t}-\mathbf{p}_{T}||,g_{t}]\\\
$ (4)

where $||.||$ represents the Euclidean distance, $p_{t}$ denotes the position
of the EE at time $t$ and $g_{t}\in\\{0,1\\}$ represents the gripper state
(open/closed). At each time, the action the agent can take is defined by an
increment motion of $0.5\beta\,mm$ in each spatial dimension, where $\beta$
varies in $\\{0,-1,+1\\}$, corresponding respectively to the conditions of no
motion, backward motion and forward motion.

The presented algorithms are trained within UnityFlexML111Project-
https://gitlab.com/altairLab/unityflexml framework, a platform to create
simulation environments supporting deformable objects and dVRK kinematics.

###  IV-A Deep Reinforcement Learning (DRL) setup

Standard DRL algorithms are trained in a simulation environment that is an
exact replica of the real one. We consider PPO as the standard baseline DRL
algorithm [19]. For the training phase, we rely on a reward function which
varies with the gripper state, to encourage EE motion towards the tumour if
the tissue has not been grasped yet (i.e. gripper is still open), and tissue
retraction if the tissue has been already grasped:

$r(s_{t})=\begin{dcases*}-||\mathbf{p}_{t}-\mathbf{q}||*k-0.5,&\text{before
grasp}\\\ -||\mathbf{p}_{t}-\mathbf{p}_{T}||*k,&\text{after
grasp}\end{dcases*}$ $None$

The scalar quantity of $-0.5$ is added to restrict the reward in the range
$(-1.0,-0.5)$ before grasping and $(-0.5,0)$ after grasping. The normalization
constant $k$ is introduced to allow re-scaling of the trajectories to
different working spaces and is inversely proportional to the maximum distance
the PSM can move.

###  IV-B Learning from Demonstrations (LfD) setup

Since we aim at finding a better policy than the provided demonstrations, the
training is based on the linear combination of the respective DRL and GAIL
losses:

$L_{Total}=\alpha L_{DRL}+\beta L_{GAIL}\\\ $ (5)

where $\alpha$ and $\beta$ represent weighting factors for the two loss
functions. For a PPO agent, $\alpha=1$ and $\beta=0$. For the GAIL agent, our
initial investigation on hyper-parameters tuning yielded best performance for
$\alpha=0.2$ and $\beta=0.8$. Other values of $\alpha$ and $\beta$ yielded
slower convergence.

Training of a GAIL agent requires the collection of trajectory demonstrations.
In this work, task demonstrations are acquired on the real dVRK and
transferred to the UnityFlexML framework. The acquired trajectories consist of
repetitive fat lifting task, performed by an expert user. Since the expert
user is well aware of the final objective of tumour exposure, the grasp
position is near the tumour for all the demonstrations. Moreover, the expert
user is instructed to diversify the trajectories by starting each
demonstration from a different initial position above the fat surface.

Figure 2: Network architecture of GAIL and PPO. PPO consists of a policy
(actor) and Value network. The policy network acts as Generator for GAIL.
Generated trajectories and expert trajectories are passed to the
Discriminator. Discriminator learns a probability function which classifies
the generator trajectory as expert or non-expert. The network layer details
are depicted inside each box in the format (hidden units, activation)
respectively.

Acquisition of task demonstrations from the real environment leverages the
communication pipeline provided by UnityFlexML (Fig. 1). Registration between
the simulated and real environment is guaranteed following the same
registration process described in [12]. The joint values are sent to
UnityFlexML through UDP sockets and the desired configuration is reached with
direct kinematics. Each recorded demonstration consists of the set of
kinematic observations that define the state space (Sec. IV) and the
corresponding actions at each timestep. An important aspect of this
implementation is the challenge associated when we reset each episode. In the
simulation, as soon as the target position is reached, the grasp is released
and the episode resets. The position of the EE is then immediately teleported
to the next initial point. This reset strategy has been adapted to cope with
the real robotic system, during the recording of expert demonstrations. In
particular, a delay of some timesteps has been added between the moment when
the grasp is released and the beginning of the next episode, to allow
repositioning. We make use of 35 continuous episodes recordings. Although our
simulation framework supports demonstration recordings using a keyboard or a
joystick, the established communication pipeline between dVRK and UnityFlexML
is crucial since it helps to acquire demonstrations directly with the real
robotic system, thus without deviating from the surgical workflow. The network
architecture used for our proposed GAIL and PPO training method is depicted in
Fig. 2.

Figure 3: Sequence of action frames for task completion in simulation (top)
and reality (bottom). From left to right: approach, grasp, retract, tumour
exposure. (top) Perspective of the simulated camera is overlaid on the bottom
left of the simulator frames. (bottom) The real camera is placed in the same
position as in the simulation (which do not correspond to the viewpoint used
to take these pictures) The colour segmentation using a circular mask is
illustrated for all the configurations (a) 0% exposure (b) 0% exposure
(c)$\sim$15% exposure (b) 100% exposure.

##  V EXPERIMENTS

We consider a tissue retraction task in the context of a partial nephrectomy
procedure. In particular, we aim at exposing a kidney tumour that is hidden by
perinephric fat tissue. The real robotic setup consists of a synthetic kidney
phantom covered with silicone fat tissue (Fig. 3). Methods presented in Sec.
IV are trained in a simulation environment which is an exact replica of the
phantom, where the deformation properties of the fat tissue have been
optimized to mimic those of the real synthetic tissue used in the experiments
[12]. All the simulation experiments, including DRL training and dVRK control,
are executed on a workstation equipped with an AMD Ryzen 3700X processor and
NVIDIA TitanX GPU.

The considered algorithms are tested both in simulation and reality based on
two different criteria: sample efficiency and the optimality of the
accomplished task, i.e. the ability of each method to expose the tumour.
Sample efficiency is estimated by the number of time steps required by each
algorithm to reach high reward values. Secondly, we estimate the optimality of
the learnt behaviour by a Tumour Exposure (TE) metric. TE is assessed as the
normalized percentage of tumour surface which can be seen from a camera
positioned in front of the kidney, both for the simulated and the real setup.
This evaluation allows us to understand if the quality of the exposure can be
maximized regardless of the PSM starting position. The behaviour learnt in the
simulated and real scenario is depicted in Fig. 3 top and Fig. 3 bottom,
respectively.

###  V-A Simulation experiments

The high level of realism of the simulated environment created within
UnityFlexML does not only allow us to train the standard DRL methods with a
sim-to-real approach but also provides a platform for testing the presented
methods in realistic settings. In order to assess whether the behaviour learnt
by the agent is robust to different starting EE positions, we perform an
experiment where the trained agent has to perform the task starting from 49
different positions uniformly sampled on a 7x7 regular grid above the portion
of the fat tissue. We evaluate TE each time the EE reaches $\mathbf{p}_{T}$.

###  V-B Real robotic experiments

Our real experimental setup (Fig. 3 bottom) is precisely aligned with the
simulated scene (Fig. 3 top) by mapping the poses of the PSM in the two
environments, as described in [20]. Accurate registration is essential to
prevent inconsistencies between the two environments, especially considering
that all the movements of the dVRK arm in the real system are controlled
through the simulated robot, including the grasping action. To compute the TE
metric, we select a circular region of interest around the tumour (exploiting
the fact that its position is fixed) and we extract the visible portion by
applying a mask with HSV bounds matching tumour colour (Fig. 3 bottom).

##  VI RESULTS AND DISCUSSION

The results obtained will be presented in two separate sections, the first
dedicated to simulation results while the second is dedicated to the results
on the real setup.

###  VI-A Simulation

Figure 4: The obtained learning curve for GAIL and PPO. Cumulative reward is
normalised in the range $[-1,0]$. The shaded area spans the range of values
obtained when training the agent starting from three different initialization
seeds. Figure 5: Simulation experiments: TE from the camera at different
initial positions, of the PSM for (a) PPO, (b) GAIL. The colour of each
subregion is related to the percentage of visible tumour area when
$\textbf{p}_{0}$ belongs to that subregion. The fat boundary from the top view
is depicted in red dashed lines whereas the fat attachment is shown in the
solid red line

Learning curves obtained with the considered learning configurations are
showed in Fig. 4. GAIL is more efficient than PPO and shows a monotonous and
smooth learning pattern. GAIL learning curve begins to increase towards high-
reward values and diverges from PPO around 1 million steps. PPO shows a
modular reward trend: it requires 2.5 million steps to learn the approach
behaviour and 1 million steps to learn the retract behaviour. This experiment
shows that incorporating human demonstrations makes learning sample-efficient
compared to baseline PPO. This result verifies our hypothesis that
incorporating human knowledge can provide initial prior reference and takes
fewer timesteps to learn the behaviour.

Figure 6: Real grasp experiments: TE from the camera when starting from
different initial positions of the EE, using (a) PPO (b) GAIL. The portion of
fat tissue which is not considered for the experiments is coloured in grey.

The plot in Fig. 5 shows the TE from the simulated camera depending on the
starting position of the PSM arm above and outside the boundary of the fat
tissue and points. In the case of PPO, when the starting EE position is close
to the fat attachment, the agent tends to grasp near the proximity of the fat
attachment, thus causing little or no tumour exposure (Fig. 5a). Note that,
for PPO, the reward function is handcrafted such that before the grasp, the
agent is encouraged to approach the known position of the tumour, because
grasping near the tumour is the best strategy to maximise tumour exposure. Low
TE indicates that manually tuning the reward function to encode complex task
objective such as tumour exposure can be challenging. We performed a
preliminary evaluation of a scenario where we design a reward function that
incorporates a TE-dependent term. However, it did not show significant
improvements in the results. This might be due to the fact that before
grasping the tumour exposure is always zero and represents a sparse reward
scenario. We plan to further investigate this in future works.

On the contrary, when human demonstrations are incorporated, GAIL is able to
grasp closer to the tumour and expose the tumour regardless of the starting
position. Note that the strategy adopted by the user while acquiring
demonstrations is to move and grasp towards points close to the tumour to
maximize the exposure. The difference between the behaviour learnt by PPO and
GAIL is in the robust selection of the grasping point when the starting
position varies. When the starting position is above the attached area, PPO
tends to grasp near the fat attachment, thus leading to low TE, whereas GAIL
learns to grasp closer to the tumour, obtaining a higher TE.

###  VI-B Real robotic setup

We have been able to successfully replicate the learned behaviour from the
simulated to the real environment without any appreciable inconsistency. The
da Vinci EE successfully gets in contact with the fat tissue for all the
different initial positions, and it is always able to reach the target point.
The tumour exposure percentage starting from various points is illustrated in
Fig. 6. For PPO, we did not initialize the EE positions near the attachment
(represented as the unattempted grey region in Fig. 6a) because grasping near
the attachment led to fat tissue tear. When comparing the results obtained for
GAIL and PPO, it emerges that GAIL is not only able to reach higher overall
exposure, but it is also more robust to changes in the initial PSM position.
In particular, tumour exposure is achieved also when starting from points that
were unattempted for PPO (Fig. 6b), thus suggesting an overall improvement in
performances, due to grasping closer to the tumour position.

This observation indicates that the initial PSM position has a great impact on
the performances in the case of PPO, whereas GAIL is able to reach optimal
performance regardless of the starting position, confirming results obtained
in the simulated experiments. In terms of the Average TE (ATE) computed
considering all trials from different starting points, PPO obtains an ATE of
0.33 in a simulated environment and 0.38 in the real, while GAIL obtains an
ATE of 0.84 in simulation and 0.90 in the real domain. We can infer that
performance using demonstrations in both simulation and real-world is robust
and outperforms PPO.

##  VII CONCLUSIONS & FUTURE WORK

In this work, we present an LfD methodology for automation of tissue
retraction based on GAIL. The proposed methodology can be trained in a
simulated replica of the real scene using a small set of real robotic
demonstrations and deployed in the real environment. The method builds on a
consolidated DRL architecture and can learn generalised human-like
trajectories in a sample-efficient way. Experiments in simulation and real
environment show that, while both baseline DRL methods and GAIL can accomplish
the task, the latter boosts the learning process by reducing the number of
steps required and learning near-human trajectories. The learnt policy has
been deployed on the dVRK and the tissue retraction task has been successfully
completed.

This work has some limitations. The underlying hypothesis is based on the
assumption of knowing the target positional coordinates (e.g., tumour
position) pre-operatively. However, the tissue retraction surgical gesture can
be carried out as an exploratory subtask without a known target. Hence, our
future work will be directed in utilizing visual information to estimate the
kinematic coordinates of the various image features, as described in [3, 21].
Further experimentation will be carried out to assess the impact of the
quality and the number of demonstrations required to learn optimal behaviour
with different surgical gestures, involving experts from various levels of
expertise. Another limitation to consider is the safety issues due to free
exploration of the state space that might lead to dangerous movements. Hence,
in subsequent work, we will incorporate safety constraints through a Safe-
Reinforcement learning technique [22].

In conclusion, surgical subtask automation is an anticipated research
direction that will lead to a higher level of autonomy and put the surgeons in
a supervisory role. This approach will enable a reduction in surgical workload
and improved patient outcomes.

## ACKNOWLEDGEMENT

The authors would like to thank Enrico Magnabosco, Sanat Ramesh for their
support in the initial development of the simulator and discussion
respectively. We also acknowledge the support of NVIDIA Corporation for the
donation of the Titan Xp GPU used in this research.

## References

  * [1] N. Simaan, R. M. Yasin, and L. Wang, “Medical technologies and challenges of robot-assisted minimally invasive intervention and diagnostics,” _Annual Review of Control, Robotics, and Autonomous Systems_ , vol. 1, pp. 465–490, 2018.
  * [2] S. Patil and R. Alterovitz, “Toward automated tissue retraction in robot-assisted surgery,” in _2010 IEEE International Conference on Robotics and Automation_. IEEE, 2010, pp. 2088–2094.
  * [3] A. Attanasio, B. Scaglioni, M. Leonetti, A. F. Frangi, W. Cross, C. S. Biyani, and P. Valdastri, “Autonomous tissue retraction in robotic assisted minimally invasive surgery–a feasibility study,” _IEEE Robotics and Automation Letters_ , vol. 5, no. 4, pp. 6528–6535, 2020.
  * [4] M. Ginesi, D. Meli, A. Calanca, D. Dall’Alba, N. Sansonetto, and P. Fiorini, “Dynamic movement primitives: Volumetric obstacle avoidance,” in _2019 19th International Conference on Advanced Robotics (ICAR)_. IEEE, 2019, pp. 234–239.
  * [5] A. Hussein, M. M. Gaber, E. Elyan, and C. Jayne, “Imitation learning: A survey of learning methods,” _ACM Computing Surveys (CSUR)_ , vol. 50, no. 2, pp. 1–35, 2017.
  * [6] J. Ibarz, J. Tan, C. Finn, M. Kalakrishnan, P. Pastor, and S. Levine, “How to train your robot with deep reinforcement learning: lessons we have learned,” _The International Journal of Robotics Research_ , p. 0278364920987859, 2021\.
  * [7] J. Ho and S. Ermon, “Generative adversarial imitation learning,” in _Advances in neural information processing systems_ , 2016, pp. 4565–4573.
  * [8] T. D. Nagy, M. Takács, I. J. Rudas, and T. Haidegger, “Surgical subtask automation—soft tissue retraction,” in _2018 IEEE 16th World Symposium on Applied Machine Intelligence and Informatics (SAMI)_. IEEE, 2018, pp. 000 055–000 060.
  * [9] N. D. Nguyen, T. Nguyen, S. Nahavandi, A. Bhatti, and G. Guest, “Manipulating soft tissues by deep reinforcement learning for autonomous robotic surgery,” in _2019 IEEE International Systems Conference (SysCon)_. IEEE, 2019, pp. 1–7.
  * [10] C. Shin, P. W. Ferguson, S. A. Pedram, J. Ma, E. P. Dutson, and J. Rosen, “Autonomous tissue manipulation via surgical robot using learning based model predictive control,” in _2019 International Conference on Robotics and Automation (ICRA)_. IEEE, 2019, pp. 3875–3881.
  * [11] F. Richter, R. K. Orosco, and M. C. Yip, “Open-sourced reinforcement learning environments for surgical robotics,” _arXiv preprint arXiv:1903.02090_ , 2019\.
  * [12] E. Tagliabue, A. Pore, D. Dall’Alba, E. Magnabosco, M. Piccinelli, and P. Fiorini, “Soft tissue simulation environment to learn manipulation tasks in autonomous robotic surgery,” in _2020 IEEE International Conference on Intelligent Robots and Systems (IROS)_. IEEE, 2020.
  * [13] C. E. Reiley, E. Plaku, and G. D. Hager, “Motion generation of robotic surgical tasks: Learning from expert demonstrations,” in _2010 Annual international conference of the IEEE engineering in medicine and biology_. IEEE, 2010, pp. 967–970.
  * [14] H. Su, A. Mariani, S. E. Ovur, A. Menciassi, G. Ferrigno, and E. De Momi, “Toward teaching by demonstration for robot-assisted minimally invasive surgery,” _IEEE Transactions on Automation Science and Engineering_ , 2021\.
  * [15] T. Osa, N. Sugita, and M. Mitsuishi, “Online trajectory planning and force control for automation of surgical tasks,” _IEEE Transactions on Automation Science and Engineering_ , vol. 15, no. 2, pp. 675–691, 2017.
  * [16] J. Schulman, A. Gupta, S. Venkatesan, M. Tayson-Frederick, and P. Abbeel, “A case study of trajectory transfer through non-rigid registration for a simplified suturing scenario,” in _2013 IEEE/RSJ International Conference on Intelligent Robots and Systems_. IEEE, 2013, pp. 4111–4117.
  * [17] A. Murali, S. Sen, B. Kehoe, A. Garg, S. McFarland, S. Patil, W. D. Boyd, S. Lim, P. Abbeel, and K. Goldberg, “Learning by observation for surgical subtasks: Multilateral cutting of 3d viscoelastic and 2d orthotropic tissue phantoms,” in _2015 IEEE International Conference on Robotics and Automation (ICRA)_. IEEE, 2015, pp. 1202–1209.
  * [18] W. Chi, G. Dagnino, T. M. Kwok, A. Nguyen, D. Kundrat, M. E. Abdelaziz, C. Riga, C. Bicknell, and G.-Z. Yang, “Collaborative robot-assisted endovascular catheterization with generative adversarial imitation learning,” in _2020 IEEE International Conference on Robotics and Automation (ICRA)_. IEEE, 2020, pp. 2414–2420.
  * [19] J. Schulman, F. Wolski, P. Dhariwal, A. Radford, and O. Klimov, “Proximal policy optimization algorithms,” _arXiv preprint arXiv:1707.06347_ , 2017\.
  * [20] N. Piccinelli, A. Roberti, E. Tagliabue, F. Setti, G. Kronreif, R. Muradore, and P. Fiorini, “Rigid 3d registration of pre-operative information for semi-autonomous surgery,” in _2020 International Symposium on Medical Robotics (ISMR)_. IEEE, 2020.
  * [21] Y. Li, F. Richter, J. Lu, E. K. Funk, R. K. Orosco, J. Zhu, and M. C. Yip, “Super: A surgical perception framework for endoscopic tissue manipulation with surgical robotics,” _IEEE Robotics and Automation Letters_ , vol. 5, no. 2, pp. 2294–2301, 2020.
  * [22] R. Cheng, G. Orosz, R. M. Murray, and J. W. Burdick, “End-to-end safe reinforcement learning through barrier functions for safety-critical continuous control tasks,” in _Proceedings of the AAAI Conference on Artificial Intelligence_ , vol. 33, no. 01, 2019, pp. 3387–3395.