LastName1, LastName2, and LastNameLastAuthor

# Optimizing Data Intensive Flows for Networks on Chips

Junwei Zhang
Uber Technoogies Inc
1191 2nd Ave 1200  Seattle  WA 98101
<EMAIL_ADDRESS>
Yang Liu
Uber Technoogies Inc
555 Market Street  San Francisco  CA  USA
<EMAIL_ADDRESS>
Li Shi
Snap Inc
Ocean Front Walk  Venice  CA
<EMAIL_ADDRESS>

Thomas G. Robertazzi
IEEE fellow
Department of Electrical and Computer Engineering
Stony Brook University
100 Nicolls Rd  Stony Brook  NY 11794
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
## Abstract

Data flow analysis and optimization is considered for homogeneous rectangular
mesh networks. We propose a flow matrix equation which allows a closed form
characterization of the nature of the minimal time solution, speedup and a
simple method to determine when and how much load to distribute to processors.
We also propose a rigorous mathematical proof about the flow matrix optimal
solution existence and that the solution is unique. The methodology introduced
here is applicable to many interconnection networks and switching protocols
(as an example we examine toroidal networks and hypercube networks in this
paper). An important application is improving chip area and chip scalability
for networks on chips processing divisible style loads.

Keywords: Divisible Load Theory, Network on Chip (NOC) Interconnection
Network, Mesh, Toroidal Networks, Hypercube Networks, Performance analysis,
Data Intensive Load.

##  1 Introduction

###  1.1 Background

Networks on chips (NOC) represent the smallest networks that have been
implemented to date [Robertazzi (2017)]. A popular choice for the
interconnection network on such networks on chips is the rectangular mesh. It
is straightforward to implement and is a natural choice for a planar chip
layout. Data to be processed can be inserted into the chip at one or more so-
called “injection points”, that is node(s) in the mesh that forward the data
to other nodes. Beyond NOCs, injecting data into a parallel processor’s
interconnection network has been done for some time, for instance in IBM’s
Bluegene machines [Krevat et al. (2002)].

In this paper it is sought to determine, for a single injection point on a
homogeneous rectangular mesh, how to optimally assign load to different
processors/links in a known timed pattern so as to process a load of data in a
minimal amount of time (i.e. minimize makespan). In this paper we succeed in
presenting an optimal technique for single point injection in homogeneous
meshes that involves no more complexity than linear equation solution. The
methodology presented here can be applied to a variety of interconnection
networks and switching/scheduling protocols besides those directly covered in
this paper. As examples, toroidal and hypercube networks are also considered
in this paper. A companion paper examines this problem with multiple sources
of load [Zhang (2018)].

Crucial to our success is the use of divisible load scheduling theory
[Bharadwaj et al. (2003)][Bharadwaj et al. (1996)]. Developed over the past
few decades, it assumes load is a continuous variable that can be arbitrarily
partitioned among processors and links in a network. Use is made of the
divisible load scheduling’s optimally principle [Bharadwaj et al. (1996)][Sohn
and Robertazzi (1996)], which says makespan is minimized when one forces all
processors to stop at the same time [Bharadwaj et al. (1996)] [Sohn and
Robertazzi (1996)] (intuitively otherwise one could transfer load from busy to
idle processors to achieve a better solution). This leads to a series of
chained linear flow and processing equations that can be solved by linear
equation techniques, often yielding recursive and even closed form solutions
for quantities such as makespan and speedup.

In this paper, the use of virtual cut-through switching [Kermani and Kleinrock
(1979)] and a modified version of store and forward switching is investigated.
These are one of many switching protocols that the methodology described here
applies to. In the virtual cut-through environment, a node can begin relaying
the first part of a message (packet) along a transmission path as soon as it
starts to arrive at the node, that is, it doesn’t have to wait to receive the
entire message before it can begin forwarding the message. In pure store and
forward switching, messages must be completely be received before being
forwarded.

More specifically, first, an equivalent processor (and makespan, speedup and
processor load fractions) is found for a $2\times 2$ homogeneous mesh network,
which can be generalized to a homogeneous $2\times n$ mesh network. After
that, we analyze the more general case of a homogeneous $m\times n$ mesh
network and obtain a general closed-form matrix representation yielding a
processor with equivalent processor speed, makespan, speedup and
processor/link load allocation. Different single data injection point
positions, such as the corner, boundary and inner grid are also discussed. In
addition, a rigorous mathematical proof about the flow matrix solution’s
existence and uniqueness is presented.

In summary, in this work, a flow matrix quantitative model, which tells one
how to deploy the data fractions to each processor in a homogeneous mesh in a
makespan optimal manner is proposed. The complexity of the technique is no
more than that of linear equation solution complexity. This work has relevance
to mesh interconnection networks used in parallel processing in general and to
meshes used in Networks on Chips in particular. An important application is
improving chip area and chip scalability for networks on chips processing
divisible style loads.

###  1.2 Related Work

In path breaking work in the 1990’s, Drozdowksi and others and created models
and largely recursive solutions for single source distribution in 2D
[Błażewicz and Drozdowski (1996)] and 3D meshes [Drozdowski and Głazek (1999)]
[Głazek (2003)], toroidal meshes [Błażewicz et al. (1999)] and hypercubes
[Błazewicz and Drozdowski (1995)]. For 2D meshes [Błażewicz and Drozdowski
(1996)] recursive solutions and closed form asymptotic results were found.
This was extended to 3D meshes with recursive solutions for load fractions
[Głazek (2003)] [Drozdowski and Głazek (1999)]. Recursive solutions for
toroidal networks and hypercubes were also found [Błażewicz et al. (1999)]
[Błazewicz and Drozdowski (1995)]. The hypercube work included a closed form
expression for speedup in terms of a fundamental load fraction assignment.

###  1.3  Our Contribution

This work is distinct form earlier work in providing matrix based solutions
(created through induction) for 2D meshes, toroidal networks and hypercubes.
Also different injection point locations in finite 2D meshes (corner, boundary
and center) are considered. Extensive simulation results based on this
modeling are presented in [Zhang (2018)].

##  2 Flow Matrix Model

###  2.1 Definitions and Assumptions

######  Definition 1.

Equivalence Computation
Equivalence computation is a technique, which consists of combining a cluster
of processors as one whole processor with equivalent processing capabilities.

The following assumptions are used throughout the paper:

  * •

Virtual cut-through [Kermani and Kleinrock (1979)] switching and store and
forward switching is used to transmit the assigned workload between
processors.

    * –

Under virtual cut-through switching, a node can relay the beginning bits of a
message (packet) before the entire message is received.

    * –

Under store and forward switching, a message must be completely received by a
node before it can be relayed to the next node along its transmission path.

  * •

For simplicity, return communication is not considered.

  * •

The communication delays are taken into consideration.

  * •

The time taken by computation and communication are assumed to be linear
function of the data size.

  * •

The network environment is homogeneous, that is, all the processors have the
same computation capacity. The link speeds between any two unit cores are
identical.

  * •

The number of outgoing ports in each processor is limited.

  * •

Single Path Communication : data transfer between two nodes follows a single
path.

The optimization objective functions is as follows :

  * •

Equivalence computation : the problem’s objective function is how to partition
and schedule the workloads among the processors to obtain the minimum makespan
(finish time).

The minimum time solution is obtained by forcing the processors over a network
to stop processing simultaneously. Intuitively, this is because the solution
could be improved by transfer load from some busy processors to idle ones
[Bharadwaj et al. (1996)] [Sohn and Robertazzi (1996)].

Processor equivalence is discussed in [Robertazzi (1993)] [Liu et al. (2007)]
and figure 1 are examples.

Figure 1: A m*n mesh network(m = 5, n = 5)

###  2.2 Notations

The following notations and definitions are utilized:

  * •

$P_{i}$: The $i$th processor. $0\leq i\leq m*n-1$.

  * •

$L$: The work load.

  * •

$D_{i}$: The minimum number of hops from the processor $P_{i}$ to the data
load injection site $L$.

  * •

$\alpha_{0}$: The load fraction assigned to the root processor.

  * •

$\alpha_{i}$: The load fraction assigned to the $i$th processor.

  * •

$\hat{\alpha_{i}}$: The load fraction assigned to each processor on the $i$th
layer $i\in 0\cdots(k-1)$.

  * •

$\omega_{i}$: The inverse computing speed on the $i$th processor.

  * •

$\omega_{eq}$: The inverse computing speed on an equivalent node collapsed
from a cluster of processors.

  * •

$r$: The rank of the flow matrix.

  * •

$z_{i}$: The inverse link speed on the $i$th link.

  * •

$T_{cp}$: Computing intensity constant. The entire load is processed in time
$\omega_{i}T_{cp}$ seconds on the $i$th processor.

  * •

$T_{cm}$: Communication intensity constant. The entire load is transmitted in
time $z_{i}T_{cm}$ seconds over the $i$th link.

  * •

$\hat{T_{f}}$: The finish time of the whole processor network. Here
$\hat{T_{f}}$ is equal to $\omega_{eq}T_{cp}$.

  * •

$T_{f}$: The finish time for the entire divisible load solved on the root
processor. Here $T_{f}$ is equal to $1\times\omega_{0}T_{cp}$, that is
$\omega_{0}T_{cp}$.

  * •

$T_{f,i}$: The finish time for the $i$th processor, $i\in 0\cdots(m*n-1)$.

  * •

$\sigma=\frac{zT_{cm}}{\omega T_{cp}}$: The ratio between the communication
speed to the computation speed, $0<\sigma<1$ [Bharadwaj et al. (1996)] [Hung
and Robertazzi (2004)].

  * •

$\sum_{i=0}^{m*n-1}\alpha_{i}=1$

  * •

$Speedup=\frac{T_{f}}{\hat{T_{f}}}=\frac{\omega T_{cp}}{\alpha_{0}\omega
T_{cp}}=\frac{1}{\alpha_{0}}$

##  3 Virtual Cut-through Switching Scenario

In the virtual cut-through environment, a node can begin relaying the first
part of a message (packets) along a transmission path as soon as it starts to
arrive at the node , that is, it doesn’t have to wait to receive the entire
message before it can begin forwarding the message.

First we consider the $2*2$ mesh network, which can be generalized to a $2*n$
mesh network. We then analyze a $m*n$ mesh network and obtain a general
closed-form matrix presentation. Finally, we give a key methodology to address
this type of question. In addition, different single data injection positions,
such as the corner, boundary and inner grid are also discussed.

###  3.1 Data Injection on The Corner Processor

####  3.1.1 2*2 Mesh Network

The load $L$ is assigned on the corner processor $P_{0}$ (figure 2). The whole
load is processed by four processors $P_{0}$, $P_{1}$, $P_{2}$, $P_{3}$
together.

Figure 2: The 2*2 mesh network and the root processor is $P_{0}$

The processor $P_{0}$, $P_{1}$ and $P_{2}$ start to process its respective
load fraction at the same time. This includes $P_{1}$ and $P_{2}$ as they are
relayed load in virtual cut-through mode at $t=0$. Because we assume a
homogeneous network (in processing speed and communication speed),
$\alpha_{1}=\alpha_{2}$ and $P_{1}$ and $P_{2}$ stop processing at the same
time. The processor $P_{3}$ starts to work when the $\alpha_{1}$ and
$\alpha_{2}$ complete transmission. That is, the link $0-1$ and $0-2$ are
occupied transmitting load to processor $1$ and $2$, respectively and only
transmission to $3$ when that is finished.

According to the divisible load theory [Bharadwaj et al. (2003)], we obtain
the timing diagram figure 3.

Figure 3: The timing diagram for 2*2 mesh network with virtual cut-through and
the root processor is $P_{0}$

Here in the Gantt-like timing diagram communication appears above each axis
and computations appears below the each axis. Let’s assume that all processors
stop computing at the same time in order to minimize the makespan [Sohn and
Robertazzi (1996)].

Based on the timing diagram, we obtain a group of linear equations to find the
fraction workload assigned to each processor $\alpha_{i}$ :

$\displaystyle\alpha_{0}\omega T_{cp}=T_{f,m}$ (1)
$\displaystyle\alpha_{1}\omega T_{cp}=T_{f,m}$ (2)
$\displaystyle\alpha_{2}\omega T_{cp}=T_{f,m}$ (3)
$\displaystyle\alpha_{1}zT_{cm}+\alpha_{3}\omega T_{cp}=T_{f,m}$ (4)
$\displaystyle\alpha_{0}+\alpha_{1}+\alpha_{2}+\alpha_{3}=1$ (5)
$\displaystyle\sigma=\frac{zT_{cm}}{\omega T_{cp}}$ (6) $\displaystyle
0<\sigma<1$ (7) $\displaystyle 0<\alpha_{0}\leq 1$ (8) $\displaystyle
0\leq\alpha_{1},\alpha_{2},\alpha_{3}<1$ (9)

The group of equations are represented by the matrix form:

${\left[\begin{array}[]{ccc}1&2&1\\\ 1&-1&0\\\
0&\sigma-1&1\end{array}\right]}\times\left[\begin{array}[]{c}\alpha_{0}\\\
\alpha_{1}\\\ \alpha_{3}\end{array}\right]=\left[\begin{array}[]{c}1\\\ 0\\\
0\end{array}\right]$ (10)

The matrix is represented as $A\times\alpha=b$. $A$ is named as the flow
matrix. Here because of symmetry $\alpha_{1}=\alpha_{2}$, so $\alpha_{2}$ is
not listed in the matrix equations.

Finally, the explicit solution is:

$\displaystyle\sigma=\frac{zT_{cm}}{\omega T_{cp}}$ (11)
$\displaystyle\alpha_{0}=\frac{1}{4-\sigma}$ (12)
$\displaystyle\alpha_{1}=\frac{1}{4-\sigma}$ (13)
$\displaystyle\alpha_{3}=\frac{1-\sigma}{4-\sigma}$ (14)

The simulation result is illustrated:

Figure 4: 2*2 mesh network. $\alpha_{0}$, $\alpha_{1}$, $\alpha_{2}$,
$\alpha_{3}$ value curve

In figure 4, the three processors $P_{0}$, $P_{1}$, $P_{2}$ have the same data
fraction workload, so the curve of $\alpha_{0}$ and $\alpha_{1}$ coincide. The
figure says that as $\sigma$ grows, the value $\alpha_{3}$ drops. In other
words, as the communication speed decreases, there is less data workload
assigned to $P_{3}$. Further, it means it will be economical to keep the load
local on $P_{0}$ $P_{1}$ $P_{2}$ and not distribute it, to other processors.
Thus for slow communication $\alpha_{0}=\alpha_{1}=\alpha_{2}=\frac{1}{3}$.

The equivalence inverse speed of a a single processor is $w_{eq}$, that can
replace the original network as

$\hat{T_{f}}=1*w_{eq}*T_{cp}$ $w_{eq}=\alpha_{0}*w$
$Speedup=\frac{T_{f}}{\hat{T_{f}}}=\frac{\omega T_{cp}}{\alpha_{0}\omega
T_{cp}}=\frac{1}{\alpha_{0}}=4-\sigma$

For a fast communication ($\sigma\approx 0$), the speedup is $4$.

####  3.1.2 2*n Mesh Network

The $2*n$ figure 5 homogeneous mesh network processes load $L$ and $L$
originates $P_{0}$.

Figure 5: 2*n (n = 10) mesh network and the workload happens on $P_{0}$

Load a distribution from $P_{0}$ to $P_{1}$ and $P_{2}$ via virtual cut-
through. After $P_{1}$ and $P_{2}$ finish receiving load from link $0-1$ and
$0-2$, they will be used to forward load to $P_{3}$ and $P_{4}$ and so on.

Similarly to the analysis of figure 3, the timing diagram for figure 5 is
shown in figure 6

Figure 6: The timing diagram for 2*10 mesh network and the data injection
happens on $P_{0}$ for virtual cut-through

The equations are presented as:

$\displaystyle\alpha_{0}\omega T_{cp}=T_{f,m}$ (15)
$\displaystyle\alpha_{1}\omega T_{cp}=T_{f,m}$ (16)
$\displaystyle\alpha_{2}\omega T_{cp}=T_{f,m}$ (17)
$\displaystyle\alpha_{1}zT_{cm}+\alpha_{3}\omega T_{cp}=T_{f,m}$ (18)
$\displaystyle\alpha_{2}zT_{cm}+\alpha_{4}\omega T_{cp}=T_{f,m}$ (19)
$\displaystyle(\alpha_{1}+\alpha_{3})zT_{cm}+\alpha_{5}\omega T_{cp}=T_{f,m}$
(20) $\displaystyle\vdots$ (21)
$\displaystyle(\alpha_{1}+\cdots+\alpha_{2\times n-1})zT_{cm}+\alpha_{2\times
n-1}\omega T_{cp}=T_{f,m}$ (22)
$\displaystyle\alpha_{0}+\cdots+\alpha_{2\times n-1}=1$ (23)
$\displaystyle\sigma=\frac{zT_{cm}}{\omega T_{cp}}$ (24) $\displaystyle
0<\sigma<1$ (25) $\displaystyle 0<\alpha_{0}\leq 1$ (26) $\displaystyle
0\leq\quad\alpha_{1}\quad\alpha_{2}\quad\cdots\quad\alpha_{2\times n-1}<1$
(27)

The flow matrix is shown:

${\left[\begin{array}[]{ccccccc}1&2&2&\cdots&2&2&1\\\ 1&-1&0&\cdots&0&0&0\\\
0&\sigma-1&1&\cdots&0&0&0\\\ 0&\sigma-1&\sigma&1&0&\cdots&0\\\
0&\sigma-1&\sigma&\sigma&1&0&0\\\ \vdots&\vdots&\vdots&\vdots&\ddots&\ddots\\\
0&\sigma-1&\sigma&\cdots&\sigma&\sigma&1\end{array}\right]}\times\left[\begin{array}[]{c}\alpha_{0}\\\
\alpha_{1}\\\ \alpha_{3}\\\ \alpha_{5}\\\ \vdots\\\ \alpha_{2\times n-3}\\\
\alpha_{2\times n-1}\end{array}\right]=\left[\begin{array}[]{c}1\\\ 0\\\ 0\\\
0\\\ \vdots\\\ 0\\\ 0\end{array}\right]$ (28)

According to the Cramer’s rule,the explicit solution for the group of
equations is:

$\displaystyle\alpha_{i}=\left|\frac{\det A^{\star}_{i}}{\det A}\right|$ (29)

where $A^{\star}_{i}$ is the matrix formed by replacing the $i$-th column of A
by the column vector b.
Specifically,

${A^{\star}_{0}=\left[\begin{array}[]{ccccccc}1&2&2&\cdots&2&2&1\\\
0&-1&0&\cdots&0&0&0\\\ 0&\sigma-1&1&\cdots&0&0&0\\\
0&\sigma-1&\sigma&1&0&\cdots&0\\\ 0&\sigma-1&\sigma&\sigma&1&0&0\\\
\vdots&\vdots&\vdots&\vdots&\ddots&\ddots\\\
0&\sigma-1&\sigma&\cdots&\sigma&\sigma&1\end{array}\right]}$ (30)

$\alpha_{0}=\left|\frac{\det A^{\star}_{0}}{\det A}\right|$ $\det
A^{\star}_{0}=-1$

The equivalence inverse processing speed :

$\hat{T_{f}}=1*w_{eq}*T_{cp}$ $w_{eq}=\alpha_{0}*w$

Finally, the speedup is:

$Speedup=\frac{T_{f}}{\hat{T_{f}}}=\frac{\omega T_{cp}}{\alpha_{0}\omega
T_{cp}}=\frac{1}{\alpha_{0}}=\left|-\det A\right|$

.

Further, we prove the matrix $\det A\neq 0$.

${C=\left[\begin{array}[]{cccccc}-1&0&\cdots&0&0&0\\\
\sigma-1&1&\cdots&0&0&0\\\ \sigma-1&\sigma&1&0&\cdots&0\\\
\sigma-1&\sigma&\sigma&1&0&0\\\ \vdots&\vdots&\vdots&\vdots&\ddots&\ddots\\\
\sigma-1&\sigma&\cdots&\sigma&\sigma&1\end{array}\right]}$ (31)

$C$ is a lower triangular matrix and the diagonal elements are not $0$. So $C$
is non-degenerate, that is, the matrix is column linear independence.

After a series of column reduction and row reduction actions, we get

${A=\left[\begin{array}[]{ccccccc}1&2&2&\cdots&2&2&1\\\ 1&-1&0&\cdots&0&0&0\\\
0&\sigma-1&1&\cdots&0&0&0\\\ 0&\sigma-1&\sigma&1&0&\cdots&0\\\
0&\sigma-1&\sigma&\sigma&1&0&0\\\ \vdots&\vdots&\vdots&\vdots&\ddots&\ddots\\\
0&\sigma-1&\sigma&\cdots&\sigma&\sigma&1\end{array}\right]}$

${\xrightarrow[\text{Reduction}]{\text{Column}}\\\
\left[\begin{array}[]{ccccccc}1&0&0&\cdots&0&0&0\\\ 1&-3&-2&\cdots&-2&-2&-1\\\
0&\sigma-1&1&\cdots&0&0&0\\\ 0&\sigma-1&\sigma&1&0&\cdots&0\\\
0&\sigma-1&\sigma&\sigma&1&0&0\\\ \vdots&\vdots&\vdots&\vdots&\ddots&\ddots\\\
0&\sigma-1&\sigma&\cdots&\sigma&\sigma&1\end{array}\right]}$

${\xrightarrow[\text{Reduction}]{\text{Row}}\\\
\left[\begin{array}[]{ccccccc}1&0&0&\cdots&0&0&0\\\ 0&-3&-2&\cdots&-2&-2&-1\\\
0&\sigma-1&1&\cdots&0&0&0\\\ 0&\sigma-1&\sigma&1&0&\cdots&0\\\
0&\sigma-1&\sigma&\sigma&1&0&0\\\ \vdots&\vdots&\vdots&\vdots&\ddots&\ddots\\\
0&\sigma-1&\sigma&\cdots&\sigma&\sigma&1\end{array}\right]}$

Considering the matrix $\hat{C}$

${\hat{C}=\left[\begin{array}[]{cccccc}-3&-2&\cdots&-2&-2&-1\\\
\sigma-1&1&\cdots&0&0&0\\\ \sigma-1&\sigma&1&0&\cdots&0\\\
\sigma-1&\sigma&\sigma&1&0&0\\\ \vdots&\vdots&\vdots&\vdots&\ddots&\ddots\\\
\sigma-1&\sigma&\cdots&\sigma&\sigma&1\end{array}\right]}$ (32)

, which is still column linear independence. Considering $0<\sigma<1$, the
flow matrix is full rank. So $\det A\neq 0$. This proof can be generalized to
$m\times n$ case.

####  3.1.3 m*n Mesh Network

Considering a general $m*n$ mesh network, such as figure 7 and figure 1.

Figure 7: 3*8 mesh network. The data injection position is $P_{0}$

Utilizing the previous methodology, we obtain the flow matrix equations for
figure 7:

$\sigma-1={\sigma}^{\star}$

${\left[\begin{array}[]{cccccccccc}1&2&3&3&3&3&3&3&2&1\\\
1&-1&0&0&0&0&0&0&0&0\\\ 0&{\sigma}^{\star}&1&0&0&0&0&0&0&0\\\
0&{\sigma}^{\star}&\sigma&1&0&0&0&0&0&0\\\
0&{\sigma}^{\star}&\sigma&\sigma&1&0&0&0&0&0\\\
0&{\sigma}^{\star}&\sigma&\sigma&\sigma&1&0&0&0&0\\\
0&{\sigma}^{\star}&\sigma&\sigma&\sigma&\sigma&1&0&0&0\\\
0&{\sigma}^{\star}&\sigma&\sigma&\sigma&\sigma&\sigma&1&0&0\\\
0&{\sigma}^{\star}&\sigma&\sigma&\sigma&\sigma&\sigma&\sigma&1&0\\\
0&{\sigma}^{\star}&\sigma&\sigma&\sigma&\sigma&\sigma&\sigma&\sigma&1\\\
\end{array}\right]}\times\left[\begin{array}[]{c}\alpha_{0}\\\ \alpha_{1}\\\
\alpha_{3}\\\ \alpha_{6}\\\ \alpha_{9}\\\ \alpha_{12}\\\ \alpha_{15}\\\
\alpha_{18}\\\ \alpha_{21}\\\
\alpha_{23}\end{array}\right]=\left[\begin{array}[]{c}1\\\ 0\\\ 0\\\ 0\\\ 0\\\
0\\\ 0\\\ 0\\\ \vdots\\\ 0\end{array}\right]$ (33)

Also, the flow matrix equations for figure 1:

${\left[\begin{array}[]{ccccccccc}1&2&3&4&5&4&3&2&1\\\ 1&-1&0&0&0&0&0&0&0\\\
0&{\sigma}^{\star}&1&0&0&0&0&0&0\\\ 0&{\sigma}^{\star}&\sigma&1&0&0&0&0&0\\\
0&{\sigma}^{\star}&\sigma&\sigma&1&0&0&0&0\\\
0&{\sigma}^{\star}&\sigma&\sigma&\sigma&1&0&0&0\\\
0&{\sigma}^{\star}&\sigma&\sigma&\sigma&\sigma&1&0&0\\\
0&{\sigma}^{\star}&\sigma&\sigma&\sigma&\sigma&\sigma&1&0\\\
0&{\sigma}^{\star}&\sigma&\sigma&\sigma&\sigma&\sigma&\sigma&1\\\
\end{array}\right]}\times\left[\begin{array}[]{c}\alpha_{0}\\\ \alpha_{1}\\\
\alpha_{3}\\\ \alpha_{6}\\\ \alpha_{10}\\\ \alpha_{15}\\\ \alpha_{19}\\\
\alpha_{22}\\\ \alpha_{24}\end{array}\right]=\left[\begin{array}[]{c}1\\\ 0\\\
0\\\ 0\\\ 0\\\ 0\\\ 0\\\ \vdots\\\ 0\end{array}\right]$ (34)

We use the similar method to prove $\det A\neq 0$. The equivalence inverse
processing speed :

$\hat{T_{f}}=1*w_{eq}*T_{cp}$ $w_{eq}=\alpha_{0}*w$

so the speedup is:

$Speedup=\frac{T_{f}}{\hat{T_{f}}}=\frac{\omega T_{cp}}{\alpha_{0}\omega
T_{cp}}=\frac{1}{\alpha_{0}}=\left|-\det A\right|$

.

The first row in the flow matrix describes the number of cores on each
$D_{i}$.

For example, there is 1 core with $0$ hop distance ($D_{0}$) from load site
$L$. There are 2 cores with $1$ hop distance ($D_{1}$) from load site $L$.
There are 3 cores with $2$ hops distance ($D_{2}$) from load site $L$, and so
on.

The number of rows means the number of different type processor data
fractions.

After these cases’ investigation, we find a crucial fact:

$\forall D_{i}=D_{j},\quad then\quad\alpha_{i}=\alpha_{j},\quad 0\leq i,j\leq
m*n-1$

##  4 Other Networks

###  4.1 Mesh Networks

##  5 Conclusion

In this work a significant problem is addressed: optimal single source load
distribution in mesh, toroidal and hypercube networks. This is done by way of
example for virtual cut-through switching and a modified version of store and
forward switching. However the approach outlined here is applicable to a wide
variety of switching/load distribution strategies and architectural
parameters. We propose a flow matrix equation method to characterize the
nature of the minimal time solution and a simple method to determine when and
how much load to distribute to processors. This work demonstrates that
mathematical modeling and tractable solutions can be an aid to designing and
evaluating scheduling strategies in parallel systems. Parallel systems will be
with us for some time so this and related work is likely to be of enduring
value.

## Acknowledgments

“This work has been submitted to the IEEE for possible publication. Copyright
may be transferred without notice, after which this version may no longer be
accessible.”

## References

  * Bharadwaj et al. (1996) Bharadwaj, V., D. Ghose, V. Mani, and T. G. Robertazzi. 1996.  Scheduling divisible loads in parallel and distributed systems, Volume 8.  John Wiley & Sons.
  * Bharadwaj et al. (2003) Bharadwaj, V., D. Ghose, and T. G. Robertazzi. 2003.  “Divisible load theory: A new paradigm for load scheduling in distributed systems”.  Cluster Computing vol. 6 (1), pp. 7–17.
  * Błazewicz and Drozdowski (1995) Błazewicz, J., and M. Drozdowski. 1995.  “Scheduling divisible jobs on hypercubes”.  Parallel Computing vol. 21 (12), pp. 1945–1956.
  * Błażewicz and Drozdowski (1996) Błażewicz, J., and M. Drozdowski. 1996.  “The performance limits of a two dimensional network of load-sharing processors”.  Foundations of Computing and Decision Sciences vol. 21 (1), pp. 3–15.
  * Błażewicz et al. (1999) Błażewicz, J., M. Drozdowski, F. Guinand, and D. Trystram. 1999.  “Scheduling a divisible task in a two-dimensional toroidal mesh”.  Discrete Applied Mathematics vol. 94 (1-3), pp. 35–50.
  * Drozdowski and Głazek (1999) Drozdowski, M., and W. Głazek. 1999.  “Scheduling divisible loads in a three-dimensional mesh of processors”.  Parallel Computing vol. 25 (4), pp. 381–404.
  * Głazek (2003) Głazek, W. 2003.  “A multistage load distribution strategy for three-dimensional meshes”.  Cluster Computing vol. 6 (1), pp. 31–39.
  * Hung and Robertazzi (2004) Hung, J. T., and T. G. Robertazzi. 2004.  “Switching in sequential tree networks”.  IEEE Transactions on Aerospace and Electronic Systems vol. 40 (3), pp. 968–982.
  * Kermani and Kleinrock (1979) Kermani, P., and L. Kleinrock. 1979.  “Virtual cut-through: A new computer communication switching technique”.  Computer Networks (1976) vol. 3 (4), pp. 267–286.
  * Krevat et al. (2002) Krevat, E., J. G. Castaños, and J. E. Moreira. 2002.  “Job scheduling for the BlueGene/L system”.  In Workshop on Job Scheduling Strategies for Parallel Processing, pp. 38–54.  Springer.
  * Liu et al. (2007) Liu, X., H. Zhao, and X. Li. June 25-28, 2007.  “Scheduling Divisible Workloads from Multiple Sources in Linear Daisy Chain Networks”.  Proceedings of the International Conference on Parallel and Distributed Processing Techniques and Applications vol. 2.
  * Robertazzi (1993) Robertazzi, T. G. 1993.  “Processor equivalence for daisy chain load sharing processors”.  IEEE Transactions on Aerospace and Electronic Systems vol. 29 (4), pp. 1216–1221.
  * Robertazzi (2017) Robertazzi, T. G. 2017.  Introduction to Computer Networking.  Springer Science.
  * Sohn and Robertazzi (1996) Sohn, J., and T. G. Robertazzi. 1996.  “Optimal divisible job load sharing for bus networks”.  IEEE Transactions on Aerospace and Electronic Systems vol. 32 (1), pp. 34–40.
  * Zhang (2018) Zhang, J. 2018.  “Data Distribution Equivalence for Data Intensive Interconnection Networks”.  Ph.D thesis Stony Brook University.

## Author Biographies

JUNWEI ZHANG received the PhD degree from the Applied Mathematics and
Statistics Department of Stony Brook University in 2018. His research
interests include parallel computing optimization, computational geometry and
applied machine learning. His email address is .

YANG LIU received his PhD degree from Department of Electrical and Computer
Engineering at Stony Brook University, Stony Brook, NY, in 2017. Previously,
he received his B.E. degree from Department of Electrical and Computer
Engineering at University of Electronic Science and Technology of China,
Chengdu, China, in 2011. His research interests are in the area of
distributed/parallel computing, networking, and load balancing algorithms. He
is currently working on divisible load theory and heterogeneous system
applications. His email address is .

LI SHI received his Ph.D. degree from Department of Electrical and Computer
Engineering at Stony Brook University, Stony Brook, NY, in 2016. Previous, he
received his B.E. degree in electrical and computer engineering from Shanghai
Jiao Tong University, Shanghai, China, in 2010. He is working at Snap Inc,
Venice, CA. His email address is .

THOMAS ROBERTAZZI is a Professor of Electrical and Computer Engineering at
Stony Brook University. He is an IEEE Fellow. He received the PhD from
Princeton University and the B.E,E, from the Cooper Union. He has published
extensively in areas such as scheduling, performance evaluation and
networking. His email address is .


# Render-in-the-loop aerial robotics simulator: Case Study on Yield Estimation
in Indoor Agriculture

Antun Ivanovic, Marsela Polic, Jelena Tabak, Matko Orsag Authors are with
Faculty of Electrical Engineering and Computing, University of Zagreb
[antun.invanovic, marsela.polic, jelena.tabak , matko.orsag] @fer.hrThis work
has been supported by Croatian Science Foundation under the project Specularia
UIP-2017-05-4042

###### Abstract

Inspired by recent promising results in sim-to-real transfer in deep learning
we built a realistic simulation environment combining a Robot Operating System
(ROS)-compatible physics simulator (Gazebo) with Cycles, the realistic
production rendering engine from Blender. The proposed simulator pipeline
allows us to simulate near-realistic RGB-D images. To showcase the
capabilities of the simulator pipeline we propose a case study that focuses on
indoor robotic farming. We developed a solution for sweet pepper yield
estimation task. Our approach to yield estimation starts with aerial robotics
control and trajectory planning, combined with deep learning-based pepper
detection, and a clustering approach for counting fruit. The results of this
case study show that we can combine real time dynamic simulation with near
realistic rendering capabilities to simulate complex robotic systems.

##  I Introduction

Today, as we witness the new dawn of Artificial Intelligence, it is hard to
imagine Unmanned Aerial Vehicle (UAV) applications without some form of
machine or deep learning. When UAVs are deployed on surveillance, search and
rescue, crop monitoring, or infrastructure inspection missions, Convolutional
Neural Networks (CNNs) are an ubiquitous part of the solution. In robotics,
simulation is the first step towards successfully deployment in real world
applications. Real time simulators like Gazebo [1], accurately and efficiently
simulate the physics and dynamics of robots, however, the photo-realistic
characteristics of images rendered this way leave a great deal to be desired.
Arguably, other game based simulators, like Xplane [2], or AirSim [3], provide
more photo-realistic virtual environments, nevertheless, their photo realism
is no match for a dedicated production renderer like Cycles [4]. Combining
physics simulation and realistic visual appearance enables us to test the
complete aerial robotics solution, starting from the low level control all the
way towards the high level mission planning and Artificial Intelligence (AI).
The difference between the two render approaches is shown in Fig 1.

We chose Cycles, from Blender, since it is an open source cross platform
solution with an integrated Python API, which makes it straightforward to
interface with the ROS environment. Our approach goes hand in hand with recent
sim-to-real attempts in deep learning, where photorealistic renderers are used
to build large synthetic datasets used to train CNNs. The proposed simulator
can be used in different indoor and outdoor scenarios, ranging from
infrastructure inspection to crop monitoring. In this work, we focus on the
yield estimation for the indoor farming. Indoor farming is becoming an
important tool to mitigate the problems farmers face with unpredictable and
often extreme weather patterns induced by global warming. Traditional indoor
agriculture enables farmers to provide their crops with optimal weather
conditions, but outdoor impact, predominantly sunlight, still plays an
important role in the plants’ growth. Modern indoor farming completely
eliminates the outside weather conditions, emulating optimal weather
conditions controlling both the climate and sunlight.

Figure 1: A segment of the greenhouse rendered with the Gazebo simulator, and
its respective counterpart rendered in Blender using the Cycles engine. The
obvious difference in image quality allows us to apply the same deep learning
based pepper detection algorithm we deploy in real world applications. When
put together, the two programs enable us to simulate complex mission
scenarios, including both low level control and high level vision based
algorithms.

While greenhouse climate control can be considered as a solved problem,
replacing manual labor is an attractive research topic for robotics. From the
robotics point of view, indoor farming provides a level of structure in the
environment which makes it easier to deploy a fully autonomous operation. The
farming industry is extremely labor intensive, and the job requirements often
fit the category of dangerous, dull, and dirty, making them ideal for
automation. Labor is even more important when considering the organic farming.
To reduce the use of pesticides, the organic agriculture requires a lot more
manual care, with a comparably smaller agricultural output. The obvious
economical consequence of such a production system is a higher cost of organic
food. Turning to robotics and automation can help alleviate the costs.

The heterogeneous team of robots in our project consists of an aerial
inspection robot, mobile robots that carry plants around the farm, and a
robotic manipulator which treats the plants [5]. However, in this paper we
focus on the aerial robot trajectory planning in a structured greenhouse
environment and fruit counting through the aforementioned Gazebo and Blender
pipeline.

#### Contributions

In this paper we focus on the yield estimation problem and present a pipeline
of technologies used to execute such a task and deploy it on the proposed
simulation framework. We present an aerial manipulator capable of capturing
images close to the plants without harming them. We disseminate the control
and trajectory planning framework necessary to deploy such an aerial
manipulator. We modified the existing deep learning based sweet pepper
detection algorithm to exploit both 2D detection results and 3D scene
information acquired with an on-board depth module. We provide the details of
the simulation environment used to simulate missions in various conditions.
Finally, we provide the results of the simulation analysis using a realistic
scenario of the yield estimation for the sweet pepper harvesting.

#### Paper organization

In the following chapter we describe the kinematics, dynamics and trajectory
planning of the aerial manipulator used in the study. Next we present the deep
learning based approach in yield estimation. Finally, we show how the proposed
simulation pipeline combines the low level and high level control. We conclude
the paper with simulation results demonstrating the effectiveness of the
proposed render-in-the-loop simulation setup.

##  II Aerial manipulator

To start, we first present an aerial manipulator equipped with an RGB-D camera
at the end-effector shown in Fig 2. Such a vehicle can fly through a
structured greenhouse environment while taking images and depth information
about the crops, without harming the plants. The RGB-D data gathered during
flight is later used for yield estimation. The structured indoor farming
environment allows us to assume the position of the UAV is accurately
measured. For example, this can be achieved using ultra-wideband technology
fused together with simultaneous localisation and mapping [6].

###  II-A Kinematics

To develop a full kinematic chain of the aerial manipulator, we first define
relevant coordinate systems, as shown in Fig. 2. The $L_{\textsc{w}}$ is the
inertial frame, $L_{\textsc{b}}$ is the body frame placed at the UAV center of
mass, $L_{0}$ is the manipulator base frame, and $L_{ee}$ is the end-effector
frame. Let $\mathbf{T}_{\textsc{a}}^{\textsc{b}}$ denote a homogeneous
coordinate transformation between the coordinate systems $L_{\textsc{a}}$ and
$L_{\textsc{b}}$. Then we can write a full kinematic chain of the aerial
manipulator as:

$\mathbf{T}_{\textsc{w}}^{ee}=\mathbf{T}_{\textsc{w}}^{\textsc{b}}\cdot\mathbf{T}_{\textsc{b}}^{\textsc{0}}\cdot\mathbf{T}_{0}^{ee},$
(1)

where $\mathbf{T}_{\textsc{w}}^{\textsc{b}}$ is the transformation between the
world and the UAV body, $\mathbf{T}_{\textsc{b}}^{\textsc{0}}$ is a fixed
transformation depending on the manipulator mounting point, and
$\mathbf{T}_{0}^{ee}$ can be obtained through the direct kinematics of the
manipulator.

Figure 2: Coordinate systems of the aerial manipulator equipped with an RGB-D
camera at the end-effector. The z-axes of the manipulator joints all point
into the image.

In our case, we consider the UAV to have six Degrees of Freedom (DoF),
$\mathbf{q}_{\textsc{b}}=\begin{bmatrix}\mathbf{p}^{T}&\mathbf{\Theta}^{T}\end{bmatrix}^{T}\in\mathbb{R}^{6}$.
The $\mathbf{p}\in\mathbb{R}^{3}$ denotes the position of the UAV in the world
frame, and the
$\mathbf{\Theta}=\begin{bmatrix}\phi&\theta&\psi\end{bmatrix}^{T}\in\mathbb{R}^{3}$
is the attitude vector containing roll, pitch and yaw angles, respectively.
The UAV used in this paper is a multirotor with all rotors placed in one
plane. From a control viewpoint, this an underactuated system because only
four DoF can be controlled simultaneously. Therefore, we define reduced
generalized coordinates of the UAV as
$\mathbf{q}_{\textsc{b}}^{*}=\begin{bmatrix}\mathbf{p}^{T}&\psi\end{bmatrix}^{T}\in\mathbb{R}^{4}$,
where $\psi$ denotes the yaw angle of the UAV. The reduced generalized
coordinates will be used later in the planning section II-C.

We also consider a serial chain manipulator with three DoF attached to the
body of the UAV. The generalized coordinates of the manipulator can be written
as
$\mathbf{q}_{\textsc{m}}=\begin{bmatrix}q_{1}&q_{2}&q_{3}\end{bmatrix}^{T}$,
where $q_{i}$ denotes the rotation of each joint. Finally, the full state of
the aerial manipulator can be written as
$\mathbf{q}=\begin{bmatrix}\mathbf{q}_{\textsc{b}}^{T}&\mathbf{q}_{\textsc{m}}^{T}\end{bmatrix}^{T}\in\mathbb{R}^{9}$.

###  II-B Dynamics and control

In previous section we defined $\mathbf{q}_{\textsc{b}}$ as generalized
coordinates of the UAV. If we take the time derivative we get
$\mathbf{\dot{q}}_{\textsc{b}}=\begin{bmatrix}\mathbf{\dot{p}}^{T}&\bm{\omega}^{T}\end{bmatrix}$,
where $\mathbf{\dot{p}}$ is the linear velocity and
$\bm{\omega}=\begin{bmatrix}\dot{\phi}&\dot{\theta}&\dot{\psi}\end{bmatrix}^{T}$
is the angular velocity, both expressed in the world frame. Note that the
angular velocity is measured in the body frame of the UAV and can be expressed
as
$\bm{\omega}_{\textsc{b}}=p\mathbf{x}_{\textsc{b}}+q\mathbf{y}_{\textsc{b}}+r\mathbf{z}_{\textsc{b}}$,
where $p$, $q$ and $r$ are angular velocities in the body frame. There exist a
simple mapping between the angular velocities in the body frame to the world
frame, for more details see [7].

In this work we consider using a multirotor vehicle with four rotors. Each
rotor produces the thrust force $F_{i}=k_{\textsc{f}}\cdot\Omega_{i}^{2}$ and
the moment $M_{i}=k_{\textsc{m}}\cdot\Omega_{i}^{2}$ along the
$\mathbf{z}_{\textsc{b}}$ axis, and $\Omega_{i}$ is the angular velocity of
each rotor. As we already mentioned, due to the underactuated nature of the
UAV, only four DoF are controllable. Therefore, we can define the UAV control
vector
$\mathbf{u}_{\textsc{uav}}=\begin{bmatrix}u_{1}&u_{2}&u_{3}&u_{4}\end{bmatrix}^{T}$,
where $u_{1}$ represents the net thrust and $u_{2}$, $u_{3}$ and $u_{4}$
represent moments around body axes. These control inputs can be written as
$\mathbf{u}_{\textsc{uav}}=\mathbf{K}\cdot\bm{\Omega}_{sq}$, where
$\mathbf{K}\in\mathbb{R}^{4\times 4}$ is the input mapping matrix, and
$\bm{\Omega}_{sq}=\begin{bmatrix}\Omega_{1}^{2}&\Omega_{2}^{2}&\Omega_{3}^{2}&\Omega_{4}^{2}\end{bmatrix}^{T}$
is the vector the squared rotors’ angular velocities.

To account for the dynamics of the manipulator, we can write the time
derivative of joint angles
$\mathbf{\dot{q}}_{\textsc{m}}=d\mathbf{q}_{\textsc{m}}/dt$. The control input
for the manipulator is a vector of torques for each joint
$\mathbf{u}_{\textsc{m}}=\begin{bmatrix}\tau_{1}&\tau_{2}&\tau_{3}\end{bmatrix}^{T}$.
Combining that with the UAV velocity vector gives the full system velocity
$\mathbf{\dot{q}}=\begin{bmatrix}\mathbf{\dot{q}}_{\textsc{b}}^{T}&\mathbf{\dot{q}}_{\textsc{m}}^{T}\end{bmatrix}^{T}\in\mathbb{R}^{9}$.
The acceleration vector can be obtained as the time derivative of the velocity
$\mathbf{\ddot{q}}=d\mathbf{\dot{q}}/dt$. The full system dynamics can now be
written as:

$\mathbf{M}({\mathbf{q}}){\mathbf{\ddot{q}}}+{\mathbf{c}}({\mathbf{q}},{\mathbf{\dot{q}}})+{\mathbf{g}}({\mathbf{q}})={\mathbf{u}},$
(2)

where $\mathbf{M}({\mathbf{q}})\in\mathbb{R}^{9\times 9}$ is the inertia
matrix, $\mathbf{c}({\mathbf{q}},{\mathbf{\dot{q}}})\in\mathbb{R}^{9}$ is the
Coriolis vector, $\mathbf{g}({\mathbf{q}})\in\mathbb{R}^{9}$ is the gravity
vector, and
$\mathbf{u}=\begin{bmatrix}\mathbf{u}_{\textsc{uav}}^{T}&\mathbf{u}_{\textsc{m}}^{T}\end{bmatrix}^{T}$
is the aerial manipulator control input. Having kinematics and dynamics of the
aerial manipulator derived gives us the foundation for the motion planning and
simulation discussed in following sections.

###  II-C Aerial manipulator motion planning

To estimate the number of fruit in the greenhouse, it is necessary to scan and
inspect the plants. In this paper, an aerial manipulator with an RGB-D camera
mounted on the end-effector is considered for the task. Since the layout of
the greenhouse is a-priori known, it is possible to scan a single plant
container unit or a whole structured row. The gathered data is then used in
the fruit counting algorithm.

####  II-C1 Scan waypoints

Each plant in a greenhouse row has to be scanned from multiple angles to
obtain informative data for fruit counting. An example of planned waypoints
and trajectory for a single row is depicted in Fig. 3. In short, the aerial
manipulator moves around the row and scans each plant with an elliptical
trajectory.

Precisely scanning each plant requires a carefully planned end-effector
motion, which mostly relies on the dimensions of growth containers. The
elliptical shape of the waypoints is based on the plant position and
dimensions. At each waypoint of the ellipse, the yaw of the UAV directs the
end-effector towards the plant’s centroid, pointing the end-effector upwards
or downwards to get a better overview of the plant.

One row of the greenhouse consists of multiple plants. To perform a scan of a
single row, it is necessary to scan each plant from both sides of the row.
Therefore, an elliptical set of waypoints is planned for each plant, based on
the dimensions and the number of plants in the row. Additional waypoints for
navigating around the row are also included to yield a smooth trajectory.

Planned waypoints are considered to form a path of $n$ waypoints:

$\mathscr{P}=\left\\{\mathbf{q}_{i}^{*}\mid\mathbf{q}_{i}^{*}\in\mathbb{R}^{7},i\in(1,2,\dots,n)\right\\},$
(3)

where
$\mathbf{q}_{i}^{*}=\begin{bmatrix}(\mathbf{q}_{\textsc{b}}^{*})^{T}&\mathbf{q}_{\textsc{m}}^{T}\end{bmatrix}^{T}$
is a single waypoint containing the position and orientation of the UAV, as
well as joint positions of the three DoF dimensional manipulator. As discussed
in previous sections, both waypoints and path are expressed in terms of
reduced generalized coordinates, due to the underactuated nature of the UAV.

Figure 3: An example of the planned elliptical trajectory used to scan a
single row of the structured greenhouse. Red lines and spheres denote the
trajectory and waypoints of the UAV body, while the yellow color denotes the
end-effector trajectory.

####  II-C2 Scan trajectory

Based on the generated path from equation (3), a trajectory that respects
dynamic constraints of the system is planned. Namely, the Time Optimal Path
Parameterization by Reachability Analysis (TOPP-RA) algorithm [8], that
operates on the numerical integration approach, is employed. As an input, the
algorithm requires a set of positions with velocity and acceleration
constraints for each DoF. The constraints are dependent on the aerial robot’s
dynamics and control. The output is a smooth trajectory $\mathscr{T}$:

$\mathscr{T}=\left\\{\mathbf{t}(t)\mid\mathbf{t}(t)\in\mathbb{R}^{21},t\in(0,t_{end})\right\\},$
(4)

where
$\mathbf{t}(t)=\begin{bmatrix}(\mathbf{q}^{*})^{T}&(\mathbf{\dot{q}}^{*})^{T}&(\mathbf{\ddot{q}}^{*})^{T}\end{bmatrix}^{T}$
denotes a trajectory point containing position, velocity and acceleration of
each DoF, and $t_{end}$ denotes the duration of the trajectory.

##  III Enabling technologies

An important drive to deploy robots in agriculture is just-in-time detection.
This is important both for agronomic and economic reasons. Farmers must have
an accurate yield estimate of their crops to plan and organize the harvest. In
this paragraph we outline the key AI technologies developed to successfully
deploy the aerial robot on a sweet pepper yield estimation task.

###  III-A Fruit detection

As most state of the art solutions, in this work, we rely on a commercial
RGB-D camera [9]. The red peppers are detected in 2D RGB images using a deep
learning model. A MobileNet based Single Shot Detector (SSD), pretrained on
the COCO dataset [10], is trained for the object detection task in 2D RGB
images. For the purposes of network training, a synthetic dataset is
procedurally generated in Blender, in order to mitigate the cost of labeling a
large training dataset [11, 12]. Increasingly popular synthetic datasets like
[13], have recently found applications in agriculture for various crops and
cultures [14, 15, 16], including a synthetic dataset for the C. annuum
semantic segmentation tasks [17]. Transfer learning for the network first
conducted on the synthetic dataset is followed by additional fine tuning on a
small dataset of real, manually labeled images. This encourages us to simulate
the entire mission in the proposed simulator setup, before the real world
deployment.

The 2D detection pipeline produces bounding boxes that denote positions of
detected peppers in the image, but the precise spatial information can be
obtained combining both detection and depth information of the RGB-D cameras.
Various methods have been developed over the recent years for 3D pose
estimation, such as surface normal estimation for grasp position optimisation
[18], and peduncle model fitting in harvesting [19]. In geometric model
fitting methods, the detected fruit is modelled with geometric primitives such
as cylinders and ellipsoids [20]. In this work, such a high precision is not
necessary, hence approximate 3D positions of the detected peppers are obtained
applying the 2D detection bounding boxes on the organised pointcloud output
from the camera depth channel. From among the filtered points representing a
single detected pepper, a centroid point can be chosen as a reasonable
approximation of the pepper position. It should be noted that this way, a
point on the pepper’s surface is chosen, resulting in varying position
estimates for the same pepper depending on the camera perspective with respect
to the pepper.

The harvesting mission can be planned based on the ripe fruit count. The
aerial manipulator executes a trajectory that enables recording the plants
from multiple perspectives. Most of the fruit, unless heavily occluded, is
recorded, and detected, from several perspectives during the UAV motion. As
stated, these detections do not match perfectly, since the positions of the
pepper surface are not at the same global position. A counting method based on
unsupervised learning is devised, that provides an estimate of the yield.

Upon trajectory execution, a set of detections is collected. In these
detections, subsets of data are separated for each greenhouse row, using known
layout of the greenhouse. Furthermore, with a predefined UAV trajectory, a
known camera frame rate, and an estimate of the possible pepper yield, random
subsampling is conducted on the detection dataset, filtering out most of false
positive detections. For an expected yield of up to 10 peppers uniformly
distributed across the plant body, this filtering step retains approximately
5-10 detections of each fruit. The remaining detections are augmented using
random noise, to a fixed size set that enables unsupervised learning methods
to properly separate the search space. We have found that 500 points are
sufficient to properly separate the data, ranging in fruit count to up to 40
peppers per row (10 fruit per plant), with detection variation at centimeter
level, and spread over approx. $2\times 0.5\times 0.5\textup{m}$. From among
the 500 points, the remaining false positive detections (e.g. pot, table, and
other detected outside the plant bounding volumes) are filtered using the
known greenhouse layout. In case the greenhouse layout is not a-priori known,
outlier removal methods can be used for the false positive filtering [21].

###  III-B Ripe fruit count

On the augmented dataset of detections, the OPTICS algorithm is deployed as a
clustering method [22] that separates the dataset by defining core points.
These core points, i.e. cluster centres, are considered as individual peppers.
The inputs to the algorithm are maximum distance $\epsilon$, and minimum
cluster size $MinPts$. The parameter $\epsilon$ represents the maximum
distance of the cluster points to the core point, in order to be considered as
a member of the cluster. In our case, this is set to $4$ cm, given the
expected size of the pepper fruit. The second parameter, $MinPts$, is a
requirement on the cluster size, i.e. on the minimum number of detections to
be considered a pepper. Empirical results in manipulator harvesting
experiments showed that at least 2 detections from various viewpoints were
needed to estimate the pepper position reliably. Due to the dataset
augmentation to the fixed size of 500 points, the minimum number of detections
is a function of the initial detections dataset size $n_{init}$, i.e. the
original requirement of 2 detections increases by the factor of $\left\lceil
500/n_{init}\right\rceil$. The method clusters the points satisfying the
provided conditions, and leaves the remaining points undefined.

Figure 4: Intermediate results of clustering based pepper fruit counting on a
subset of the greenhouse scene. The initial set of detections generated by
applying the CNN detector results to the scene point cloud are augmented with
a Gaussian random noise until a dataset of fixed size 500 is generated. The
OPTICS algorithm separates the augmented search space, identifying 7 clusters
(7 peppers). The remaining undefined points are removed for clarity.

Thanks to the organised structure of the greenhouse, the positions of pepper
growth containers (pots) are known (or can be known if a similar detection was
deployed for pot detection). Although this is not crucial in this paper, it is
an important feature for the ultimately envisioned case of autonomous robotic
team in a greenhouse environment. This information can be used to select a
growth unit to be harvested and the detected pepper positions can be used as
a-priori information for harvesting with a manipulator.

##  IV Simulation setup

Figure 5: The functional diagram representing the system pipeline. At the
start, a request for scanning the plants is generated, upon which the aerial
manipulator executes a scan trajectory in the greenhouse. The gathered camera
position and orientation is used through Blender to generate a realistic
dataset for fruit counting. The dataset is afterwards processed in ROS to
obtain the synchronized point cloud and RGB image for each camera
configuration. Using that data, the fruit count is estimated.

One of the key features which make Gazebo one of the most popular physics
simulators in robotics is its simple and straight-forward interface with ROS.
However, due to its limitations in photorealistic visualisation and the
complexity of modeling a world containing a large number of detailed models,
we propose extending the simulation pipeline with Blender production rendering
software, Cycles. The complete pipeline is described with a schematic in Fig.
5.

###  IV-A ROS/Gazebo

As depicted with the schematic in Fig. 5, for the most part, our solution to
the proposed yield estimation is implemented within the ROS environment. This
setup enables us to simply switch between the simulation and the real hardware
when it is ready to be deployed. The aerial manipulator performs a scanning
trajectory, so that the end-effector tip position is used to track the camera
motion. The tracked pose of the camera is later imported in Blender, and a
photorealistic image of the environment is rendered using Cycles on the exact
pose the camera had in the simulation. Based on the rendered camera image,
fruit yield is estimated per plant, and provided to the mission planner. The
mission planner generates a schedule for the robot(s) involved in the
greenhouse cultivation. The speed of the mission execution is a function of
available resources and the setup. Mission planning goes beyond the scope of
this paper, but interested reader can find more details in [23].

The rotorcraft UAV is simulated in Gazebo using RotorS simulator [1] to
actuate the body. The UAV is modeled as a rigid body with four propeller units
mounted at each of the UAV’s arms. The aforementioned three DoF manipulator is
attached to the UAV body. More information about this setup and the
corresponding GitHub repository is available at [24]. Low level position
control both for the UAV body and the manipulator joints are developed in the
form of classic cascade PID controllers.

A simplified 3D model of the greenhouse is used during the Gazebo simulation.
Such a setup allows us to simulate the trajectory in real time, while keeping
track of potential collisions. Different plugins enable emulating various
sensors needed for localization and control. This also allows one to simulate
the sensors with variable degree of noise and accuracy, which is important for
troubleshooting the solution and it makes the simulation more realistic.

While performing the simulation, the end-effector configuration is recorded
and saved to a file. This is necessary to provide the Blender environment with
positions and corresponding orientations of the RGB-D camera for rendering,
which is explained in detail in the following section.

###  IV-B Blender

The model of the greenhouse consists of rows of tables carrying pepper growth
containers. Each row consists of 4 tables carrying a single pepper plant. The
plants are modeled using the basic building blocks developed for synthetic
dataset generation, namely realistic pepper models, planar leaf models with
realistic textures, flower models, plant stem, and a plastic pot. Additional
advantage Blender has over Gazebo at this point is that the plants are
generated procedurally with a random particle generator. This means that in
different trials one can vary the pepper, the flower and the leaf number,
position, and orientation with respect to the plant stem. More info on the
setup we use to train the AI to detect peppers can be found in [25], where we
successfully demonstrated pepper detection method trained using both synthetic
and actual pepper images. The 3D assets used to build this synthetic dataset
were used to create random 3D models of the greenhouse. Example of these 3D
models are shown in Fig 6.

Figure 6: Example models used to build the 3D models of the plants in the
simulated greenhouse.

Once the trajectory is simulated using the simplified model in Gazebo, we
record the position and orientation of the RGB-D camera. Camera poses are
imported in Blender, where each sample from the Gazebo simulation is rendered
using Cycles. Since we can vary the number of peppers and the shape of the
plants, this setup allows us to obtain multiple variable attempts using the
same trajectory generated in the Gazebo.

###  IV-C Simulated RGB-D camera

Another important role Blender plays in this paper is the animated camera
system, simulating an RGB-D sensor for the Gazebo robot simulation. The
simulation pipeline links the motion of the end-effector of the aerial
manipulator, with an eye-in-hand RGB-D camera, to camera animation in the
Blender environment model. The RGB image of the scene is rendered using built-
in Blender renderer, Cycles. The depth image is generated from the depth map
calculated within Blender Cycles rendering engine. The depth map is a part of
the depth buffer, also known as a z-buffer, usually used in computer graphics
to represent depth information of objects from a particular perspective [26].
Cycles is a ray tracing computer graphics algorithm that uses depth buffers to
render a scene ensuring that the correct polygons are shown from the provided
camera perspective.

Modifying the algorithm in [27], we were able to render depth and RGB images
from camera viewpoints and corresponding camera positions. This data is
afterwards processed offline in ROS, transforming the recorded depth image
into an organized point cloud [28]. This framework loads RGB and depth image
and publishes them sequentially, together with the intrinsic camera
parameters, camera position and orientation. We rely on the standard
depth_image_proc ROS package that generates point cloud from the depth image,
and publishes it accordingly. The point cloud is published along with the
corresponding RGB image and global camera transformation. This format
corresponds to the output of a real RGB-D camera, and can be used by the
detection and yield estimation software solution. Extrinsic parameters of the
camera are required so that the pepper detection and clustering algorithm can
estimate the yield of the simulated greenhouse. The detection and counting is
implemented in Python, relying on the TensorFlow and Pythons Sklearn module.

##  V Simulation results and conclusion

The greenhouse layout we used for simulation is shown in Fig. 7. The structure
consists of eight tables, each with four pepper plant containers. The
stationary manipulator workstation is located in the center of the structure.
For this layout, we ran several simulations with a random number of peppers
per plant. For every run the plants were generated randomly. We tested the
entire proposed pipeline generating the UAV trajectory for greenhouse
inspection, counting the fruits from the video rendered in Blender, and
planning the mission based on the inputs from the pepper counting.

Figure 7: The layout of the proposed greenhouse structure consisting of eight
tables, each with four pepper plant containers. The stationary manipulator
workstation is located in the center of the structure.

To test the planner with a random number of fruits per pepper plant, we
selected $20$ of the $40$ available plants that bore fruit. For these plants,
a random number of peppers between $1$ and $5$ was generated from a uniform
distribution for the left and right sides of the plant. The simulation setups
for this use case in terms of number of peppers ranged from 3 to 6 peepers per
plant. For every scenario the yield estimation was correct with up to 10%
margin of error.

The results show that _Render-in-the-loop_ setup can be used to simulate both
the dynamics of the system and vision algorithms for the proposed yield
estimation scenario. The complexity of the proposed scenario allows us to
assume this simulator pipeline can be used in different scenarios in various
indoor and outdoor environments. Future work includes simulating various
infrastructure inspection missions as well as outdoor agriculture
applications.

## References

  * [1] F. Furrer, M. Burri, M. Achtelik, and R. Siegwart, Robot Operating System (ROS): The Complete Reference (Volume 1), ch. RotorS—A Modular Gazebo MAV Simulator Framework, pp. 595–625.  Cham: Springer International Publishing, 2016.
  * [2] R. Garcia and L. Barnes, “Multi-uav simulator utilizing x-plane,” in Selected papers from the 2nd International Symposium on UAVs, Reno, Nevada, USA June 8–10, 2009, pp. 393–406, Springer, 2009.
  * [3] P. Pueyo, E. Cristofalo, E. Montijano, and M. Schwager, “Cinemairsim: A camera-realistic robotics simulator for cinematographic purposes,” in 2020 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), pp. 1186–1191, IEEE, 2020.
  * [4] B. Iraci, Blender cycles: lighting and rendering cookbook.  Packt Publishing Ltd, 2013.
  * [5] M. Polic, A. Ivanovic, B. Maric, B. Arbanas, J. Tabak, and M. Orsag, “Structured ecological cultivation with autonomous robots in indoor agriculture,” in 2021 16th International Conference on Telecommunications (ConTEL), pp. 189–195, IEEE, 2021.
  * [6] R. Milijas, L. Markovic, A. Ivanovic, F. Petric, and S. Bogdan, “A comparison of lidar-based slam systems for control of unmanned aerial vehicles,” in 2021 International Conference on Unmanned Aircraft Systems (ICUAS), pp. 1148–1154, IEEE, 2021.
  * [7] J. Macdonald, R. Leishman, R. Beard, and T. McLain, “Analysis of an improved IMU-based observer for multirotor helicopters,” Journal of Intelligent & Robotic Systems, vol. 74, pp. 1049–1061, May 2013.
  * [8] H. Pham and Q. Pham, “A new approach to time-optimal path parameterization based on reachability analysis,” IEEE Transactions on Robotics, vol. 34, pp. 645–659, June 2018.
  * [9] L. Fu, F. Gao, J. Wu, R. Li, M. Karkee, and Q. Zhang, “Application of consumer rgb-d cameras for fruit detection and localization in field: A critical review,” Computers and Electronics in Agriculture, vol. 177, p. 105687, 2020.
  * [10] T. Lin, M. Maire, S. J. Belongie, L. D. Bourdev, R. B. Girshick, J. Hays, P. Perona, D. Ramanan, P. Dollár, and C. L. Zitnick, “Microsoft COCO: common objects in context,” CoRR, vol. abs/1405.0312, 2014.
  * [11] S. Hinterstoisser, O. Pauly, H. Heibel, M. Martina, and M. Bokeloh, “An annotation saved is an annotation earned: Using fully synthetic training for object detection,” in Proceedings of the IEEE/CVF International Conference on Computer Vision Workshops, pp. 0–0, 2019.
  * [12] S. Khan, B. Phan, R. Salay, and K. Czarnecki, “Procsy: Procedural synthetic dataset generation towards influence factor studies of semantic segmentation networks.,” in CVPR Workshops, pp. 88–96, 2019.
  * [13] A. Barisic, F. Petric, and S. Bogdan, “Sim2air-synthetic aerial dataset for uav monitoring,” IEEE Robotics and Automation Letters, 2022.
  * [14] M. Di Cicco, C. Potena, G. Grisetti, and A. Pretto, “Automatic model based dataset generation for fast and accurate crop and weeds detection,” in 2017 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), pp. 5188–5195, IEEE, 2017.
  * [15] J. Olatunji, G. Redding, C. Rowe, and A. East, “Reconstruction of kiwifruit fruit geometry using a cgan trained on a synthetic dataset,” Computers and Electronics in Agriculture, vol. 177, p. 105699, 2020.
  * [16] K. Zhang, Q. Wu, and Y. Chen, “Detecting soybean leaf disease from synthetic image using multi-feature fusion faster r-cnn,” Computers and Electronics in Agriculture, vol. 183, p. 106064, 2021.
  * [17] R. Barth, J. IJsselmuiden, J. Hemming, and E. J. Van Henten, “Data synthesis methods for semantic segmentation in agriculture: A capsicum annuum dataset,” Computers and electronics in agriculture, vol. 144, pp. 284–296, 2018.
  * [18] C. Lehnert, A. English, C. McCool, A. W. Tow, and T. Perez, “Autonomous sweet pepper harvesting for protected cropping systems,” IEEE Robotics and Automation Letters, vol. 2, no. 2, pp. 872–879, 2017.
  * [19] I. Sa, C. Lehnert, A. English, C. McCool, F. Dayoub, B. Upcroft, and T. Perez, “Peduncle detection of sweet pepper for autonomous crop harvesting—combined color and 3-d information,” IEEE Robotics and Automation Letters, vol. 2, no. 2, pp. 765–772, 2017.
  * [20] C. Lehnert, I. Sa, C. McCool, B. Upcroft, and T. Perez, “Sweet pepper pose detection and grasping for automated crop harvesting,” in 2016 IEEE International Conference on Robotics and Automation (ICRA), pp. 2428–2434, IEEE, 2016.
  * [21] M. M. Breunig, H.-P. Kriegel, R. T. Ng, and J. Sander, “Lof: identifying density-based local outliers,” in Proceedings of the 2000 ACM SIGMOD international conference on Management of data, pp. 93–104, 2000.
  * [22] M. Ankerst, M. M. Breunig, H.-P. Kriegel, and J. Sander, “Optics: Ordering points to identify the clustering structure,” ACM Sigmod record, vol. 28, no. 2, pp. 49–60, 1999.
  * [23] B. Arbanas Ferreira, T. Petrović, M. Orsag, J. Ramiro Martínez-de Dios, and S. Bogdan, “Distributed allocation and scheduling of tasks with cross-schedule dependencies for heterogeneous multi-robot teams,” arXiv e-prints, pp. arXiv–2109, 2021.
  * [24] B. Arbanas, A. Ivanovic, M. Car, M. Orsag, T. Petrovic, and S. Bogdan, “Decentralized planning and control for uav–ugv cooperative teams,” Autonomous Robots, Feb 2018.
  * [25] M. Polić, J. Tabak, and M. Orsag, “Pepper to fall: a perception method for sweet pepper robotic harvesting,” Intelligent Service Robotics, 2021.
  * [26] T. Theoharis, G. Papaioannou, and E.-A. Karabassi, “The magic of the z-buffer: A survey,” 2001.
  * [27] J. Martin, “Blendertorgbd.” https://github.com/JavonneM/BlenderToRGBD, 2019\.
  * [28] A. Ivanovic, “Blender rgbd ros.” https://github.com/larics/blender_rgbd_ros, 2019.