# Squeeze flow of micro-droplets: convolutional neural network with trainable
and tunable refinement

Aryan Mehboudi  Shrawan Singhal  S.V. Sreenivasan

(Nov. 15, 2022)

###### Abstract

We propose a platform based on neural networks to solve the image-to-image
translation problem in the context of squeeze flow of micro-droplets. In the
first part of this paper, we present the governing partial differential
equations to lay out the underlying physics of the problem. We also discuss
our developed Python package, sqflow, which can potentially serve as free,
flexible, and scalable standardized benchmarks in the fields of machine
learning and computer vision. In the second part of this paper, we introduce a
residual convolutional neural network to solve the corresponding inverse
problem: to translate a high-resolution (HR) imprint image with a specific
liquid film thickness to a low-resolution (LR) droplet pattern image capable
of producing the given imprint image for an appropriate spread time of
droplets. We propose a neural network architecture that learns to
systematically tune the refinement level of its residual convolutional blocks
by using the function approximators that are trained to map a given input
parameter (film thickness) to an appropriate refinement level indicator. We
use multiple stacks of convolutional layers the output of which is translated
according to the refinement level indicators provided by the directly-
connected function approximators. Together with a non-linear activation
function, such a translation mechanism enables the HR imprint image to be
refined sequentially in multiple steps until the target LR droplet pattern
image is revealed. The proposed platform can be potentially applied to data
compression and data encryption. The developed package and datasets are
publicly available on GitHub at https://github.com/sqflow/sqflow.

Keywords: Machine learning, Squeeze flow, Lubrication theory, Fluid dynamics,
Data compression, Data encryption, Image-to-image translation, Multi-class
multi-label classification

###### Contents

  1. 1 Introduction
    1. 1.1 Imprint lithography
    2. 1.2 Droplet pattern image & Imprint image
  2. 2 Motivation
  3. 3 Physics-based solver
    1. 3.1 Governing equations
    2. 3.2 Parameters
  4. 4 Dataset
    1. 4.1 Data generation
    2. 4.2 Breakdown of the compiled dataset
    3. 4.3 Structure of each dataset
  5. 5 A gentle dive into the fluid dynamics of the system
  6. 6 Network structure
    1. 6.1 Function approximation
    2. 6.2 Convolutional blocks with tunable refinement with (CNN-TR-DS) and without (CNN-TR) down-scaling
    3. 6.3 Output preparation block
  7. 7 Results and discussion
    1. 7.1 Pre-processing
    2. 7.2 Training
    3. 7.3 Function approximators
    4. 7.4 Sequential refinement
    5. 7.5 Kernels
    6. 7.6 Custom image
  8. 8 Future works

##  1 Introduction

###  1.1 Imprint lithography

The imprint lithography and particularly a subclass of it called step and
flash imprint lithography (SFIL) [1], are considered as promising alternatives
to photolithography process. In this method, the imprint resist droplets are
dispensed on a substrate, followed by imprinting the droplets using a
superstrate (template), and curing the resulted liquid film into a desired
pattern [2, 1].

Inspired by SFIL, we aim at developing a methodology for producing a target
imprint image from an appropriate droplet pattern map by using a blank
superstrate. Let us consider a pair of parallel plates (substrate and
superstrate) initially separated from each other as shown in Fig. 1.

Figure 1: A schematic representation of the imprint lithography with blank
superstrate: (A) dispensing imprint resist droplets according to a droplet
pattern, (B) alignment of superstrate (template) and substrate as needed
followed by squeezing the droplets, (C) spread and merge of imprint resist
droplets, (D) polymerization of imprint resist using the UV exposure, and (E)
separation of template from substrate.

First, one or multiple droplets of imprint resist are dispensed on the
substrate using an inkjet dispenser according to a droplet pattern map. The
superstrate is then brought in contact with the liquid droplets. Because of
the capillary force, the droplets start to spread and merge. This flow of
liquid continues for a period of time, the length of which is called spread
time. The resulted liquid film is then polymerized by using an ultraviolet
(UV) exposure. The superstrate is debonded from the substrate at the end. A
schematic representation of these steps can be found in Fig. 1. The main
difference of this methodology from SFIL is that in this work, we consider a
blank superstrate aiming at finding the appropriate initial position of
droplets that can produce a given pattern with a specific film thickness.

###  1.2 Droplet pattern image & Imprint image

In practice, an inkjet printhead has an array of nozzles separated by a finite
distance, which is typically of the order of tens of microns. The pitch of the
nozzle array dictates the pixel size of the droplet pattern image. The droplet
pattern image is simply a representation of the firing nozzles. For each
nozzle that dispenses a droplet, the corresponding pixel is considered to have
the value of ‘1’ (On). Otherwise, a value of ‘0’ (Off) is assigned to the
corresponding pixel.

On the other hand, the imprint image is a top-down view of the liquid film.
The areas with liquid are assigned a value of ‘1’ (Wet), while a value of ‘0’
is considered for areas without liquid (Dry). We have developed a physics-
based model to simulate the squeeze flow of droplets for any given droplet
pattern (forward problem). The physical domain of interest is discretized into
a finite number of computational cells. Therefore, the pixel size of the
imprint image is dictated by the size of the computational cells.

##  2 Motivation

The schematic representation depicted in Fig. 2 illustrates the forward and
inverse problems of the squeeze flow of droplets.

Figure 2:  A general schematic representation of the forward (solid black
arrows) and inverse (orange dashed arrows) problems for the squeeze flow of
droplets.

Although the forward problem, i.e. prediction of the 1. imprint image and 2.
film thickness from a given tuple of 1. drop pattern image and 2. spread time,
is achievable using a physics-based solver, the inverse problem may not be
straightforward. That is, given an imprint image with a specific film
thickness, it may not be trivial to find the appropriate drop pattern image
and the required spread time.

As an example shown in Fig. 3, two imprint images are produced with two
different film thicknesses for a given drop pattern image under two different
spread times.

Figure 3: A droplet pattern image together with the (numerically) produced
imprint images for two different spread times as well as the schematic
representations of the cross section view of the liquid film. The solid black
arrows show the forward problem, while the orange dashed arrows denote the
inverse problem. Here, the droplet pattern image has $20\times 20$ pixels,
while the imprint images consist of $160\times 160$ pixels. The liquid’s
thermophysical properties and other related parameters are discussed in the
text.

It can be discerned that the droplet pattern map itself can be a good model
representation of the imprint image at the beginning of the spread of droplets
(after applying an appropriate transformation to match the dimensions of
images). However, as time elapses, the error associated with such a crude
model that ignores the film thickness can become noticeably larger.

The goal of this work is primarily to address the inverse problem in the
context of squeeze flow of droplets. However, we believe that providing our
compiled dataset as a flexible standardized benchmark may potentially benefit
the community of machine learning and computer vision for testing algorithms
on time series images, image-to-image translation [3], multi-class multi-label
classifications [4, 5, 6], image super-resolution [7, 8, 9, 10, 11], image
deblurring [12, 13, 14], etc. Therefore, we decided to publicly distribute the
developed package and the compiled datasets on GitHub at
https://github.com/sqflow/sqflow.

##  3 Physics-based solver

In this section, we discuss the governing equations for the squeeze flow of
droplets. We also provide a high-level description of the physics-based solver
that we have developed for the forward problem.

###  3.1 Governing equations

From the lubrication theory [15], the velocity ($V_{l}$) and pressure ($p$)
fields within a liquid film are coupled as

$\vec{V_{l}}=-\frac{h^{2}}{12\mu_{l}}\vec{\nabla}p,$ (1)

where $\mu_{l}$ denotes the liquid’s viscosity, and $h$ shows the
instantaneous local liquid film thickness. After change of variable
$\hat{p}\equiv p-p_{amb}+(\cos\theta_{1}+\cos\theta_{2})\sigma/h$, wherein
$\theta_{1}$ and $\theta_{2}$ denote the contact angles on the top and bottom
plates, $\sigma$ shows the surface tension, and $p_{amb}$ refers to the
ambient pressure, the governing equation describing the spread of liquid film
confined between two plates can be written as

$\begin{split}\vec{\nabla}.\Big{(}h^{3}\vec{\nabla}\hat{p}\Big{)}&=12\mu_{l}\frac{\partial
h}{\partial t},\\\ \text{at interface:}\quad\hat{p}&=0.\end{split}$ (2)

The boundary value problem in Eqs. 2 needs to be solved for an appropriate
instantaneous gap reduction rate $\partial h/\partial t$. The instantaneous
gap reduction rate is determined so that the force balance is satisfied as
described in [16].

The solver we have developed for the forward problem is based on the finite
volume method (FVM). The volume fraction $f$ is defined for each computational
cell as the ratio of liquid volume in the cell to the total volume of that
cell. After finding the instantaneous pressure and velocity fields, the volume
fraction field needs to be modified according to the velocity field to reflect
the evolution of the liquid film over time. In order to achieve that, the
following partial differential equation is solved:

$\begin{split}\frac{\partial f^{\ast}}{\partial
t}+\vec{\nabla}.(f^{\ast}\vec{V_{l}})&=0,\end{split}$ (3)

wherein $f^{\ast}$ is the modified volume fraction parameter defined as
$f^{\ast}\equiv fh/h_{\text{ref}}$, and $h_{\text{ref}}$ is a reference gap
distance. Given the instantaneous velocity and modified volume fraction fields
at time $t$, the solver then solves Eq. 3 to determine the modified volume
fraction field at time $t+\Delta t$. The volume fraction ($f$) field can then
be obtained for the next time step. The same process is repeated for each time
step until at least one of the termination criteria is met. A more detailed
description of this implementation can be found in [16].

###  3.2 Parameters

The dataset has been compiled for an array of $20\times 20$ nozzles with a
pitch of $84.5\leavevmode\nobreak\ \mu m$. We have also set the size of each
computational cell to be $1/8$ of the pitch of the nozzle array, i.e.
$10.5625\leavevmode\nobreak\ \mu m$, which results in imprint images with
$160\times 160$ pixels. In addition, we have considered a uniform droplet
volume (6 pl) when producing the dataset. Some of the other important
parameters are shown below:

  * •

Surface tension: $0.032$ N/m.

  * •

Viscosity: $0.001$ Pa.s.

  * •

Initial thickness of film: $1\leavevmode\nobreak\ \mu m$.

##  4 Dataset

###  4.1 Data generation

The physics-based solver has been used to generate data from different droplet
pattern images. The criteria for terminating each simulation includes:

  * •

The liquid film covering more than $90\%$ of the field, i.e. any imprint image
has less than $90\%$ of its pixels ‘On’.

  * •

The spread time exceeding 1 second, i.e. all the examples in the dataset have
a spread time ranging from 0 to 1 second.

  * •

The film thickness decreasing below $5\leavevmode\nobreak\ nm$, i.e. all the
examples in the dataset have a film thickness ranging from
$5\leavevmode\nobreak\ nm$ to $1\leavevmode\nobreak\ \mu m$.

Based on the aforementioned termination criteria, each simulation for a given
droplet pattern image ($20\times 20$ pixels) produces approximately about
40–60 examples for different spread times. In addition, the liquid can flow
out of the system partially. In particular, this leads to fewer examples for
the cases of having one or a few droplets close to the corner with the
coordinates of $(20,20)$.

Some of the examples generated for a random droplet pattern image are shown in
Fig. 4.

Figure 4:  Some of the examples from the dataset: a drop pattern image with
$20\times 20$ pixels (red border) together with the evolution of the liquid
film (imprint images with $160\times 160$ pixels) for various spread times:
$0$, $0.6\leavevmode\nobreak\ ms$, $1.8\leavevmode\nobreak\ ms$,
$4.2\leavevmode\nobreak\ ms$, and $8.0\leavevmode\nobreak\ ms$, increasing
from left to right and top to bottom. The corresponding liquid film
thicknesses are $1\leavevmode\nobreak\ \mu m$, $263\leavevmode\nobreak\ nm$,
$179\leavevmode\nobreak\ nm$, $140\leavevmode\nobreak\ nm$, and
$116\leavevmode\nobreak\ nm$, decreasing from left to right and top to bottom.

We can see that the droplets spread and the film thickness reduces as time
elapses. This phenomenon happens fairly quickly at the beginning, but it slows
down over time.

###  4.2 Breakdown of the compiled dataset

We have run the physics-based solver to generate data for several categories
based on the number of ‘On’ pixels in the droplet pattern image, i.e. number
of firing nozzles of the printhead to dispense droplet. More information about
the breakdown of these datasets can be found in Table 1.

Category | 1 | 2 | 3 | 4 | 5 | 7 | 10 | 15 | 20 | 25 | 30 | 40
---|---|---|---|---|---|---|---|---|---|---|---|---
#Simulations | 400 | 100 | 100 | 100 | 100 | 100 | 100 | 100 | 100 | 100 | 50 | 50
#Examples | 17,047 | 4,912 | 5,627 | 5,882 | 6,189 | 6,441 | 6,203 | 5,745 | 5,101 | 5,028 | 2,295 | 2,007
$\overline{|DP|}$ | 1.00 | 2.00 | 3.00 | 4.00 | 4.97 | 6.94 | 9.84 | 14.80 | 19.54 | 24.21 | 29.03 | 38.27
$\overline{|DP|}$: Average count number of ‘On’ pixels in the drop pattern
images.
Table 1: Breakdown of the compiled datasets.

The dataset for each category is provided independently from others. Some of
the advantages include:

  * •

Enabling researchers to select the categories that may best suit their
application in future.

  * •

Providing a convenient way to properly prepare the training, validation, and
test datasets from the data related to different categories not seen by each
other.

  * •

Reducing the size of dataset by splitting it into smaller partitions.

When preparing the training, validation, and test datasets for a supervised
learning task, one or multiple categories can be used partially or in its
entirety. However, in order to mitigate the data leakage issue, one should be
careful to NOT include examples with the same droplet pattern image in more
than one dataset. This can be achieved by using data from different categories
when preparing the training, validation, and test datasets. Alternatively, one
can preprocess the data to ensure that the examples associated with the same
droplet pattern image are not included in more than one of the training,
validation, and test datasets.

Fig. 5 presents the breakdown of the compiled datasets for different
categories in terms of the count number of occurrences for each pixel in the
droplet pattern images.

(a) Category: 1

(b) Category: 2

(c) Category: 3

(d) Category: 4

Figure 5:  Breakdown of the compiled datasets: the count number of occurrences
for each pixel in the droplet pattern images with $20\times 20$ pixels.

(a) Category: 5

(b) Category: 7

(c) Category: 10

(d) Category: 15

Figure 6: (continued) Breakdown of the compiled datasets: the count number of
occurrences for each pixel in the droplet pattern images with $20\times 20$
pixels.

(a) Category: 20

(b) Category: 25

(c) Category: 30

(d) Category: 40

Figure 7: (continued) Breakdown of the compiled datasets: the count number of
occurrences for each pixel in the droplet pattern images with $20\times 20$
pixels.

###  4.3 Structure of each dataset

Each dataset includes four csv files:

  * •

“t”–consisting of one column of data. The $i$th row holds the value of spread
time (unit: second) for the $i$th example in that dataset partition.

  * •

“h”–consisting of one column of data. The $i$th row holds the value of film
thickness (unit: meter) for the $i$th example in that dataset partition.

  * •

“dp”–consisting of multiple columns. The $i$th row holds the index of all ‘On’
pixels (value of 1) in the droplet pattern image for the $i$th example in that
dataset partition. These pixels denote the firing nozzles of the printhead to
dispense droplet. Other pixels are ‘Off’ (value of 0).

  * •

“vof”–consisting of multiple columns. The $i$th row holds the index of all
‘On’ pixels (value of 1) in the imprint image for the $i$th example in that
dataset partition. These pixels denote the wet area of the field. Other pixels
are ‘Off’ (value of 0).

The developed Python package comes with some simple functions to read the
dataset partition(s) from a given directory, or by walking through a root
directory, and returning the compiled t, h, dp, and vof as NumPy arrays. The
dp array has a shape of (dataset size, 400) since each drop pattern image has
400 pixels. The vof array has a shape of (dataset size, 25600) since each
imprint image has $160\times 160$ pixels.

##  5 A gentle dive into the fluid dynamics of the system

Figure 8: The temporal evolution of the film thickness for three different
drop pattern images presented in Fig. 9.

Figure 9: The liquid film map (imprint image vof) for different spread times
given three different drop pattern images (dp): Example 1 (Top), Example 2
(Middle), and Example 3 (Bottom). The solid black arrows show the translation
of drop pattern image to an imprint image for different spread times (forward
problem). The orange dashed arrows denote the translation of an imprint image
with a given film thickness to an appropriate drop pattern image and the
required spread time (inverse problem). Here, the droplet pattern images have
$20\times 20$ pixels, while the imprint images contain $160\times 160$ pixels.

Before using the dataset, it is good to have some level of familiarity with
the underlying physics of the system the dataset is derived from. For example,
an accurate and robust image-to-image mapping is not achievable by considering
solely the dp and vof images in its entirety. Any reasonable mapping between
these images in its entirety should capture the strongly-nonlinear correlation
of these images with each other and also with the spread time and liquid film
thickness. Alternatively, one may omit the effects of spread time (or liquid
film thickness) by confining the dataset to a narrow range of spread time (or
liquid film thickness).

In this section, we aim at explaining, from a high-level perspective, how the
parameters of an example are correlated qualitatively: t (scalar), h (scalar),
dp (drop pattern image), and vof (imprint image).

Let us start with a simple system with only one droplet (Example 1) as shown
in Fig. 9 (Top). As we can see, only one pixel of the dp is ‘On’. At the
beginning of the imprinting process, the imprint image resembles reasonably
the drop pattern image, except that it has $160\times 160$ pixels compared to
$20\times 20$ pixels of the drop pattern image. As time elapses, the droplet
spreads further and the liquid film thickness decreases (conservation of mass
principle).

The dynamics of the system is relatively fast at the beginning, but it becomes
slower with time. For example, it takes about $1.27\leavevmode\nobreak\ ms$
for the liquid film thickness to reduce from $1\leavevmode\nobreak\ \mu m$ to
$140\leavevmode\nobreak\ nm$, while it needs another $7.13\leavevmode\nobreak\
ms$, i.e. total spread time of $8.4\leavevmode\nobreak\ ms$, to get to
$54\leavevmode\nobreak\ nm$ film thickness.

Now, let us consider a dp image with only five ‘On’ pixels (Example 2) as
shown in Fig. 9 (Middle). The dynamics of the system is the same as the
previous one until the droplets start to merge, i.e. it takes about
$1.27\leavevmode\nobreak\ ms$ for the liquid film thickness to reduce from
$1\leavevmode\nobreak\ \mu m$ to $140\leavevmode\nobreak\ nm$. However, the
current system becomes slower than the previous one as a result of the merge
of droplets. This behavior can be discerned more easily from Fig. 8 which
shows the dynamic evolution of the liquid film thickness.

Similarly, the same fluid flow characteristics can be observed for a dp image
with 29 ‘On’ pixels (Example 3) as shown in Fig. 9 (Bottom), except that the
merge of droplets that are, now, more in quantity and closer to each other in
some areas, causes the dynamics of the system to become slower than in Example
2. Fig. 8 should help more easily observe this behavior.

In summary, it can be discerned from Figs. 9 and 8 that:

  * •

The dynamics of the system becomes slower with time.

  * •

The dynamics of the system becomes slower as a result of the merge of
droplets.

##  6 Network structure

The basic principle of our proposed model is to train and tune the refinement
level of convolutional layers. At the core of this convolutional neural
network with trainable and tunable refinement level (CNN-TTR) is residual
convolutional blocks with tunable refinement (CNN-TR) connected to function
approximators that are trained to find the optimal refinement level. The
refinement level of CNN-TR blocks is then tuned by the trained function
approximator. A schematic representation of our proposed CNN-TTR model is
depicted in Fig. 10.

Figure 10: (A) The architecture of our residual convolutional neural network
with trainable and tunable refinement level (CNN-TTR) to map a high-resolution
imprint image (VOF), in conjunction with a normalized film liquid thickness
$h^{\ast}$, to its target low resolution droplet pattern image (DP), as well
as (B) the building blocks used in our model.

Our convolutional neural network consists of components to accomplish the
following tasks:

  * •

Function approximation: dense blocks of $f_{1},f_{2},...$ aiming at predicting
the appropriate refinement level in each stage.

  * •

Down-scaling: a stack of convolutional layers with tunable refinement and
down-scaling (CNN-TR-DS). An input image with $160\times 160$ pixels is down-
scaled to $80\times 80$ (Block #1), $40\times 40$ (Block #2), and $20\times
20$ pixels (Block #3) sequentially.

  * •

Widening of field of view: a stack of convolutional layers with tunable
refinement and no down-scaling (CNN-TR).

  * •

Output preparation: a convolutional layer with tunable refinement and sigmoid
activation function.

###  6.1 Function approximation

The similarity between an imprint image vof and its target drop pattern image
dp decreases as the droplets spread and merge further with time. As a result,
in order to predict the target dp image, more refinement (removing ‘On’ pixels
properly) of the vof image is required. Finding the appropriate refinement
level in each stage is accomplished by a function approximator. Such function
approximators aim at finding the appropriate refinement levels to be used by
their corresponding residual convolutional blocks with tunable refinement
(CNN-TR and CNN-TR-DS). Here, a function approximator [17] is a dense neural
network with one hidden layer of 128 units and activation function of ReLU
[18].

###  6.2 Convolutional blocks with tunable refinement with (CNN-TR-DS) and
without (CNN-TR) down-scaling

Each CNN-TR and CNN-TR-DS block has three inputs: 1\. feature maps from the
preceding block, 2\. a single map consisting of the element-wise average of
the feature maps from the preceding block, and 3. a refinement level indicator
provided by the corresponding function approximator. Providing the average of
the feature maps of the last preceding block enables access to the
corresponding hierarchical features and mitigate the loss of information that
should be preserved.

Each CNN-TR and CNN-TR-DS block consists of a stack of three connected Conv.
Norm. blocks. The output of this stack is added (element-wise) to inputs #2
and 3: the refinement level indicator and the average feature map from the
last preceding block. The difference of CNN-TR and CNN-TR-DS is that the
latter uses a max pooling layer for down-scaling, while CNN-TR blocks are used
to widen the field of view.

Each Conv. Norm. block comprises a convolutional layer with the kernel size of
5, followed by an instance normalization [19], and ReLU activation function
[18]. For the CNN-TR-DS blocks, each convolutional layer has 32, 64, and 128
filters in Blocks #1, 2, and 3, respectively. For the CNN-TR blocks, each
convolutional layer has 128, 64, and 32 filters in Blocks #4, 5, and 6,
respectively.

###  6.3 Output preparation block

Block #7 is similar to a CNN-TR block with the difference that it has only one
convolutional layer with no normalization or ReLU activation. It has a sigmoid
activation function at the end of the block. The convolutional layer in this
block has 8 filters with a kernel size of 5.

##  7 Results and discussion

###  7.1 Pre-processing

When preparing the training, validation, and test datasets, we selected data
from different categories according to Table 2.

Category | 1 | 2 | 3 | 4 | 5 | 7 | 10 | 15 | 20 | 25 | 30 | 40
---|---|---|---|---|---|---|---|---|---|---|---|---
Training | 0.125 | - | - | 1.0 | - | 1.0 | 1.0 | - | 1.0 | 1.0 | - | 1.0
Validation | 0.125 | - | - | - | 1.0 | - | - | 1.0 | - | - | - | -
Test | 0.125 | - | - | - | - | - | - | - | - | - | 1.0 | -
Table 2: Fraction of each category selected for pre-processing when preparing
training, validation, and test datasets.

To mitigate the issue with examples that can potentially have multiple
solutions, we omitted the examples with too large continuous liquid film in
their vof image. We displaced an interrogation window with $72\times 72$
pixels across each vof image. Any example with a local coverage greater than
$90\%$ was removed from the datasets. The final number of examples in
training, validation, and test datasets are shown in Table. 3.

Training | Validation | Test
---|---|---
12,476 | 5,703 | 2193
Table 3: Number of examples in training, validation, and test datasets after
pre-processing.

We use $t^{\ast}$ and $h^{\ast}$ as the normalized scalars for $t$ and $h$,
respectively, as defined below:

$\begin{split}t^{\ast}&=\frac{\ln{t}-\mu_{t}}{\sigma_{t}}\\\\[5.0pt]
h^{\ast}&=\frac{\ln{h}-\mu_{h}}{\sigma_{h}},\end{split}$ (4)

wherein $\mu_{t}$ and $\sigma_{t}$ refer to the arithmetic mean and standard
deviation of $\ln{t}$ for the examples within the training dataset,
respectively. Similarly, $\mu_{h}$ and $\sigma_{h}$ refer to the arithmetic
mean and standard deviation of $\ln{h}$ for the examples within the training
dataset, respectively. The distribution of $h^{\ast}$ and $t^{\ast}$ is shown
in Fig. 11 for the training, validation, and test datasets.

(a) Training dataset

(b) Validation dataset

(c) Test dataset

Figure 11:  The distribution of $h^{\ast}$ and $t^{\ast}$ for the training,
validation, and test datasets.

###  7.2 Training

The model has 3,085,935 trainable parameters in total. We used TensorFlow [20]
and Keras [21] to train our model. We considered the binary cross entropy loss
while applying a $L_{2}$ regularization penalty to bias and kernel of the
convolutional layers.

We use the precision-recall (PR) plot and the associated area under curve
(AUC) metric to evaluate the performance of our model. In average, only about
$3\%$ of pixels of dp images are ‘On’. We noticed that in our case, because of
having an imbalanced dataset, the AUC of receiver operating characteristics
(ROC) misleadingly leads to a much higher value than the AUC-PR. Therefore,
the AUC-ROC could be deceptive when drawing conclusions about the performance
of the model on the current dataset wherein the number of negatives outweighs
the number of positives significantly [22].

The validation dataset was used to examine the performance of the trained
model for various classification thresholds. The results are shown in Fig. 12
(a) in terms of precision, recall, and F-1 Score.

(a)

(b)

Figure 12:  (a) The variations of precision, recall, and F-1 score with the
classification threshold when running the model on the validation dataset. (b)
The precision-recall plot as well as (Inset) the F-1 score variations for the
magnified region of interest.

The PR curve is also presented in Fig. 12 (b). The obtained AUC-PR for the
validation dataset equals $0.9956$. The PR curve shows an acceptable
performance over a wide range of thresholds, but the highest F-1 score related
to the validation dataset is found when the threshold is about $0.505$. The
corresponding F-1 score is $0.9735$. Finally, by using this optimized model
and classification threshold, we obtain the F-1 score of $0.9647$ on the test
dataset.

We used the test dataset to compare the trained model’s performance with that
of a crude model consisting of three max pooling layers with no trainable
parameters. The results are shown in Table. 4.

Model | AUC-PR | Precision | Recall | F-1
---|---|---|---|---
Baseline crude model | 0.1159 | 0.1159 | 1.0000 | 0.0886
CNN-TTR | 0.9951 | 0.9734 | 0.9562 | 0.9647
Table 4: Comparison of the performances of the model CNN-TTR and a crude model
consisting of three max pooling layers without any trainable parameter when
acting on the test dataset. Here, the threshold equals 0.5 when calculating
the precision, recall, and F-1 score.

As expected, the crude model can predict all the ‘ON’ pixels of the dp images
(Recall=1.0) at the cost of a relatively low precision of about 11.59%.

###  7.3 Function approximators

The output of function approximators $f_{1},f_{2},...,f_{7}$ is shown in Fig.
13 for various $h^{\ast}$ values.

(a)

(b)

Figure 13:  (a) The output of function approximators for various $h^{\ast}$
values using the trained model CNN-TTR and (b) its counterpart with a linear
activation function within the CNN-TR and CNN-TR-DS blocks and no
normalization layer in the Conv. Norm. blocks.

Fig. 13 (a) shows the results for the trained model CNN-TTR the schematic
representation of which was shown in Fig. 10. When drawing conclusions about
these results, we need to keep in mind that the model CNN-TTR consists of an
instance normalization layer within the Conv. Norm. blocks and the nonlinear
activation function of ReLU within the CNN-TR and CNN-TR-DS blocks, as a
result of which the observed trends my not be easily interpretable.
Nevertheless, the plot shows that, in general, the outputs are negative.
Together with the immediate ReLU activation function within the CNN-TR and
CNN-TR-DS blocks, such negative translations lead to a sequential refinement,
i.e removal of ‘On’ pixels of image in a sequential fashion.

Fig. 13 (b) shows the results for a less effective (slower training and lower
AUC-PR) variant of CNN-TTR with a linear activation function within the CNN-TR
and CNN-TR-DS blocks and no normalization layer within the Conv. Norm. blocks.
Because of omitting the non-linearity associated with the aforementioned
components, the behavior of function approximators is, now, more
interpretable: a vof image with a thinner liquid film (smaller $h^{\ast}$)
requires more refinement to obtain its target dp image. Hence, the function
approximators output a smaller value (negative number with a larger magnitude)
as $h^{\ast}$ decreases to lead to a larger refinement level (removal of more
‘On’ pixels).

###  7.4 Sequential refinement

The concept of sequential refinement is demonstrated in Fig. 14 for a given
example from the test dataset.

Figure 14:  The average response maps of Blocks #1–6 as well as the logits and
the sigmoid activation map related to Block # 7, together with the ground-
truth dp image for a given example from the test dataset with
$h^{\ast}=-0.001061$.

The model refines the vof image in multiple steps based on the normalized
thickness of the liquid film $h^{\ast}$ until the predicted dp image is
obtained. It can be discerned that a region with a relatively small puddle of
liquid is more easily refined, while an area with a large continuous liquid
film is more challenging to be refined properly.

There are two potential culprits for the problem associated with the large-
coverage regions: 1\. insufficient examples during the training of model, and
2\. non-uniqueness of solution for such regions. In particular, a large-
coverage region my be obtainable by multiple droplet patterns that differ only
slightly in one or a few pixels (usually close to the center of such large-
coverage regions).

Figure 15:  Prediction of the model CNN-TTR for several examples from the test
dataset. Here, the classification threshold equals 0.4 when obtaining the
“Predicted DP” from the sigmoid layer’s activation. From top to bottom,
$h^{\ast}$ is -2.641, -2.574, -2.435, -0.173, -0.136, 0.317, 0.169, and 0.123.

One of the consequences of the model getting less confident in these
situations is that some of the ‘On’ pixels can be lost when using the
classification threshold optimized for the highest F-1 score. In such
situations, reducing the classification threshold can lead to a higher recall
value at the cost of reducing the precision and F-1 score. However, as shown
in Fig. 12, the model performs acceptably well on a wide range of the
classification threshold. Hence, one can reduce the classification threshold
to improve the recall metric. We have noticed that the subsequent reduction of
precision is partly due to the slight misplacement of one or few ‘On’ pixels
close to the center of large-coverage regions, which is still much better than
completely loosing those ‘On’ pixels. The results obtained using a
classification threshold of 0.4 are shown in Fig. 15 for some of the examples
from the test dataset.

###  7.5 Kernels

Some of the kernels of the first convolutional layer in Blocks #1–7 are shown
in Fig. 16.

(a)

(b)

(c)

Figure 16:  Some of the kernels of the first convolutional layer for different
Blocks (normalized to vary between 0 and 1).

(a)

(b)

(c)

Figure 17: (continued) Some of the kernels of the first convolutional layer
for different Blocks (normalized to vary between 0 and 1).  (a)  Figure 18:
(continued) Some of the kernels of the first convolutional layer for different
Blocks (normalized to vary between 0 and 1).

The kernels of the earlier blocks act as a collection of various liquid-gas
interface detectors. As moving higher, the kernels act more like a set of
detectors for various local droplet patterns.

###  7.6 Custom image

In order to examine how the model performs on a custom image, we used an image
of the letter ‘A’ as shown in Fig. 19 (Left).

Figure 19:  (Left) A custom image of the letter ‘A’ together with (Right) the
droplet pattern produced by a crude model consisting of three max pooling
layers with no trainable parameters.

A crude model consisting of three max pooling layers with no trainable
parameters predicts a droplet pattern with 97 ‘On’ pixels (number of droplets
to be dispensed) as shown in Fig. 19 (Right).

We used the CNN-TTR model to predict the droplet pattern images for different
target liquid film thicknesses (classification threshold of 0.3). Some of the
results are shown in Fig. 20.

(a)

(b)

(c)

(d)

Figure 20:  (Left) The droplet pattern image (DP) predicted by the CNN-TTR
model for the custom image of the letter ‘A’ shown in Fig. 19 together with
the (Right) corresponding imprint image obtained from the physics-based solver
using the predicted droplet pattern image for different cases wherein the
target liquid film thickness equals (a) 60 nm, i.e. $h^{\ast}=-0.641$, (b) 70
nm, i.e. $h^{\ast}=-0.469$, (c) 90 nm, i.e. $h^{\ast}=-0.189$, and (d) 200 nm,
i.e. $h^{\ast}=0.701$.

It can be perceived that the model generally predicts a fewer number of
droplets (‘On’ pixels in dp image) to be dispensed for thinner liquid films
(10–20 droplets) since the droplets spread more while the liquid film
thickness reduces to a smaller value.

In the next step, we examined the droplet patterns predicted for these inverse
problems, i.e. left images in Fig. 20 (a–d). We used the physics-based solver
to find the imprint vof image for the predicted droplet pattern and the
corresponding target liquid film thickness: right images in Fig. 20 (a–d).

We notice that the obtainable imprint images do not have as high resolution as
the input image of the letter ‘A’ shown in Fig. 19. This can be attributed to
the finite size of droplets (6 pl in this work) and the fact that they become
even larger (lower obtainable resolution) as the liquid film is squeezed
further. As a result, the best obtainable resolution is dictated by the size
of droplets dispensed by the printhead.

The time series imprint vof images obtained by using the physics-based solver
for the case of the target liquid film thickness of 90 nm are shown in Fig.
21.

Figure 21:  The droplet pattern produced by the CNN-TTR model for a custom
image of the letter ‘A’ shown in Fig. 19 and the target liquid film thickness
of 90 nm, i.e. $h^{\ast}=-0.189$, together with the time series imprint vof
images obtained from the physics-based solver for the predicted droplet
pattern image.

The corresponding results related to the droplet pattern predicted by a crude
model consisting of three max pooling layers with no trainable parameters are
also shown in Fig. 22.

Figure 22:  The droplet pattern produced by a crude model consisting of three
max pooling layers with no trainable parameters for a custom image of the
letter ‘A’ shown in Fig. 19 together with the time series imprint vof images
obtained from the physics-based solver for the predicted droplet pattern
image.

We can see that the crude model’s performance deteriorates as the liquid film
thickness decreases since it does not take into account the spread and merge
of droplets.

##  8 Future works

The results of the inverse problems shown in Fig. 20 can be improved by
human/physics-based solver intervention, i.e. using the model’s predicted
droplet pattern image as a starting point to add, remove, or displace droplets
for a more desired outcome. Hence, the future efforts can be devoted to two
purposes: 1. improving the neural network model, and 2. enhancing the neural
network model by including the physics-based solver in the decision making
process. Therefore, some of the potential future works could include:

  * •

including both neural network (inverse problem) and physics-based solver
(forward problem) in the decision making process.

  * •

extending the dataset to include smaller droplets that can potentially lead to
higher resolution patterns.

  * •

extending the dataset to include various droplet diameters to increase the
flexibility of generating droplet patterns for higher resolution patterns.

## Conclusion

In this paper, we proposed a residual convolutional neural network with a new
structure suitable for refining a high resolution image to a low resolution
one. We demonstrated the concept in the context of squeeze flow of micro-
droplets. We showed that the CNN-TTR model learns to systematically tune the
refinement level of its residual convolutional blocks by using multiple
function approximators. In future, the datasets can be extended to include
various droplet volumes to enhance the flexibility of generating droplet
pattern maps for producing higher resolution imprint patterns. We foresee that
the proposed platform will be used in data compression and data encryption.
Moreover, it may pave the way for a light-weight LR data communication which
can be decoded to its corresponding HR format by using the corresponding
forward model on each side of the communication. In addition, the developed
package and compiled datasets can potentially serve as free, flexible, and
scalable standardized benchmarks in connection with time series images, multi-
class multi-label classifications with many classes/labels, image deblurring,
etc.

## Data Availability Statement

The developed package and datasets are publicly available on GitHub at
https://github.com/sqflow/sqflow.

## Acknowledgments

The authors acknowledge the Texas Advanced Computing Center (TACC) at The
University of Texas at Austin for providing HPC resources that have
contributed to the research results reported within this paper. URL:
http://www.tacc.utexas.edu

## Conflict of Interest

The authors have no conflicts to disclose.

## References

  * [1] S. Sreenivasan, “Nanoimprint lithography steppers for volume fabrication of leading-edge semiconductor integrated circuits,” Microsystems & Nanoengineering 3 no. 1, (Dec., 2017) . http://www.nature.com/articles/micronano201775.
  * [2] S. Reddy and R. T. Bonnecaze, “Simulation of fluid flow in the step and flash imprint lithography process,” p. 200.  San Jose, CA, May, 2005.  http://proceedings.spiedigitallibrary.org/proceeding.aspx?doi=10.1117/12.601413.
  * [3] Y. Pang, J. Lin, T. Qin, and Z. Chen, “Image-to-Image Translation: Methods and Applications,” July, 2021.  http://arxiv.org/abs/2101.08629. arXiv:2101.08629 [cs].
  * [4] J. Lanchantin, T. Wang, V. Ordonez, and Y. Qi, “General Multi-label Image Classification with Transformers,” Nov., 2020.  http://arxiv.org/abs/2011.14027. arXiv:2011.14027 [cs].
  * [5] S. Yun, S. J. Oh, B. Heo, D. Han, J. Choe, and S. Chun, “Re-labeling ImageNet: from Single to Multi-Labels, from Global to Localized Labels,” July, 2021.  http://arxiv.org/abs/2101.05022. arXiv:2101.05022 [cs].
  * [6] W. Liu, H. Wang, X. Shen, and I. W. Tsang, “The Emerging Trends of Multi-Label Learning,” IEEE Transactions on Pattern Analysis and Machine Intelligence 44 no. 11, (Nov., 2022) 7955–7974. http://arxiv.org/abs/2011.11197. arXiv:2011.11197 [cs].
  * [7] Z. Wang, J. Chen, and S. C. H. Hoi, “Deep Learning for Image Super-resolution: A Survey,” Feb., 2020.  http://arxiv.org/abs/1902.06068. arXiv:1902.06068 [cs].
  * [8] J. Li, Z. Pei, and T. Zeng, “From Beginner to Master: A Survey for Deep Learning-based Single-Image Super-Resolution,” Sept., 2021.  http://arxiv.org/abs/2109.14335. arXiv:2109.14335 [cs, eess].
  * [9] A. Liu, Y. Liu, J. Gu, Y. Qiao, and C. Dong, “Blind Image Super-Resolution: A Survey and Beyond,” July, 2021.  http://arxiv.org/abs/2107.03055. arXiv:2107.03055 [cs].
  * [10] C. Saharia, J. Ho, W. Chan, T. Salimans, D. J. Fleet, and M. Norouzi, “Image Super-Resolution via Iterative Refinement,” June, 2021.  http://arxiv.org/abs/2104.07636. arXiv:2104.07636 [cs, eess].
  * [11] B. C. Maral, “Single Image Super-Resolution Methods: A Survey,” Feb., 2022.  http://arxiv.org/abs/2202.11763. arXiv:2202.11763 [cs, eess].
  * [12] H. Zhao, Z. Ke, N. Chen, S. Wang, K. Li, L. Wang, X. Gong, W. Zheng, L. Song, Z. Liu, D. Liang, and C. Liu, “A new deep learning method for image deblurring in optical microscopic systems,” Journal of Biophotonics 13 no. 3, (Mar., 2020) . https://onlinelibrary.wiley.com/doi/10.1002/jbio.201960147.
  * [13] X. Tao, H. Gao, X. Shen, J. Wang, and J. Jia, “Scale-Recurrent Network for Deep Image Deblurring,”.
  * [14] K. Zhang, W. Ren, W. Luo, W.-S. Lai, B. Stenger, M.-H. Yang, and H. Li, “Deep Image Deblurring: A Survey,” May, 2022.  http://arxiv.org/abs/2201.10700. arXiv:2201.10700 [cs].
  * [15] B. Hamrock, Fundamentals of Fluid Film Lubrication.  NASA Reference Publication 1255, 1991.
  * [16] A. Mehboudi, S. Singhal, and S. V. Sreenivasan, “Modeling the squeeze flow of droplet over a step,” Physics of Fluids 34 no. 8, (Aug., 2022) 082005. https://aip.scitation.org/doi/10.1063/5.0098597. Publisher: American Institute of Physics.
  * [17] K. Hornik, M. Stinchcombe, and H. White, “Multilayer feedforward networks are universal approximators,” Neural Networks 2 no. 5, (Jan., 1989) 359–366. https://linkinghub.elsevier.com/retrieve/pii/0893608089900208.
  * [18] X. Glorot, A. Bordes, and Y. Bengio, “Deep Sparse Rectiﬁer Neural Networks,” Artificial Intelligence and Statistics 15 (2011) 9.
  * [19] D. Ulyanov, A. Vedaldi, and V. Lempitsky, “Instance Normalization: The Missing Ingredient for Fast Stylization,” Nov., 2017.  http://arxiv.org/abs/1607.08022. arXiv:1607.08022 [cs].
  * [20] M. Abadi, A. Agarwal, P. Barham, E. Brevdo, Z. Chen, C. Citro, G. S. Corrado, A. Davis, J. Dean, M. Devin, S. Ghemawat, I. Goodfellow, A. Harp, G. Irving, M. Isard, Y. Jia, R. Jozefowicz, L. Kaiser, M. Kudlur, J. Levenberg, D. Mane, R. Monga, S. Moore, D. Murray, C. Olah, M. Schuster, J. Shlens, B. Steiner, I. Sutskever, K. Talwar, P. Tucker, V. Vanhoucke, V. Vasudevan, F. Viegas, O. Vinyals, P. Warden, M. Wattenberg, M. Wicke, Y. Yu, and X. Zheng, “TensorFlow: Large-Scale Machine Learning on Heterogeneous Distributed Systems,” Mar., 2016.  http://arxiv.org/abs/1603.04467. arXiv:1603.04467 [cs].
  * [21] F. Chollet and others, “Keras,” 2015.  https://keras.io.
  * [22] T. Saito and M. Rehmsmeier, “The Precision-Recall Plot Is More Informative than the ROC Plot When Evaluating Binary Classifiers on Imbalanced Datasets,” PLOS ONE 10 no. 3, (Mar., 2015) e0118432. https://dx.plos.org/10.1371/journal.pone.0118432.


# Learning Emergent Behavior in Robot Swarms with NEAT

Pranav Rajbhandari  Computer Science Department
Carnegie Mellon University
Pittsburgh, PA
<EMAIL_ADDRESS>
ORCID: 0009-0004-4933-5204  Donald Sofge  Head, Distributed Autonomous Systems
Section
U.S. Naval Research Laboratory
Washington, DC
<EMAIL_ADDRESS>
ORCID: 0000-0003-0153-3581

###### Abstract

When researching robot swarms, many studies observe complex group behavior
emerging from the individual agents’ simple local actions. However, the task
of learning an individual policy to produce a desired emergent behavior
remains a challenging and largely unsolved problem. We present a method of
training distributed robotic swarm algorithms to produce emergent behavior.
Inspired by the biological evolution of emergent behavior in animals, we use
an evolutionary algorithm to train a ‘population’ of individual behaviors to
approximate a desired group behavior. We perform experiments using simulations
of the Georgia Tech Miniature Autonomous Blimps (GT-MABs) aerial robotics
platforms conducted in the CoppeliaSim simulator. Additionally, we test on
simulations of Anki Vector robots to display our algorithm’s effectiveness on
various modes of actuation. We evaluate our algorithm on various tasks where a
somewhat complex group behavior is required for success. These tasks include
an Area Coverage task, a Surround Target task, and a Wall Climb task. We
compare behaviors evolved using our algorithm against designed policies, which
we create in order to exhibit the emergent behaviors we desire.

###### Index Terms:

Swarm robotics, Emergent behavior, Evolutionary algorithm, Neuroevolution

##  I Introduction

###  I-A Emergent Behavior

Emergent behavior is a phenomenon observed in swarms of agents. It is
generally defined as a complex swarm behavior which occurs as a consequence of
each individual agent following a relatively simple control scheme [13, 15].
Examples of this can be found in the behavior of groups of animals, such as
how some species of fire ants have been observed to create rafts out of their
bodies to survive flooding [11]. In this case, the individual behavior of
‘grab onto neighbors’ creates the emergent behavior of ‘build a raft’.
Emergent behavior is also exhibited in migrating swarms of some species of
caterpillars, which walk on top of each other in order to create a ‘rolling
swarm’ [19]. This allows the caterpillars to migrate quicker than simply
walking. Similar instances of complex emergent behavior have been observed in
the field of swarm robotics.

###  I-B Swarm Robotics

Though there are various definitions of swarm robotics, it is generally agreed
that robot swarms are multi-agent systems where each agent is autonomous and
has low complexity [3, 8, 12]. They are used for various tasks, including
exploration and surveillance. They are characterized by using locally
communicating distributed systems as opposed to a central controller. Due to
this, they remain functional as the swarm size is increased. However, they do
not have a central controller, so careful fine tuning of the individual
policies is required to achieve a desired swarm behavior.

It is well known that emergent behavior can arise in robot swarms, and past
studies have explored this. Pagello et al. [16] (2001) study how to make a
robot swarm perform a cooperative task by creating emergent behaviors. They
implement this by dynamically assigning predefined roles to each robot. A
function $Q$ is learned for each agent, which takes in local information and
outputs the best role for the agent to adopt. After refining their $Q$
functions, they observe cooperative emergent behavior in their chosen task,
robot soccer.

In a more recent study, Oliveri et al. [14] use a Monte Carlo scheme to
continuously refine the behavior of individual agents in a swarm. The robots
they choose to use are arranged in a line, and connected with soft tubes.
Their actuation is in the form of a pump, which increases the distance between
a robot and the neighbor behind it. In training, each robot attempts to learn
its optimal phase time (the time that it takes for a cycle of the pump filling
then unfilling the tube). The robots randomly sample nearby phase times, and
update depending on the change in velocity they experience. The study
concludes that the result of each agent attempting to optimize its velocity is
an emergent behavior where the group of connected robots crawls forward.

Similar to the research of Oliveri and Pagello, we focus on optimizing an
individual policy to best create a desired emergent swarm behavior. Inspired
by how biological evolution was effective in creating various emergent
behaviors in animals, we propose using an evolutionary algorithm to solve this
task.

###  I-C Evolutionary Algorithms/NEAT

Evolutionary algorithms are algorithms inspired by biological evolution, where
a fitness function is optimized by repeatedly mutating and evaluating members
of a ‘population’ of solutions. Importantly, these methods are generic, as
they do not make any assumptions about the format of the solutions or the
topography of the fitness function. Since we want to solve a problem with
general inputs/outputs, we use neuroevolution (NE), a type of evolutionary
algorithm that evolves neural networks. In early NE algorithms, the
topology111The topology of a neural network refers to the graph structure of
the nodes and the connections between them of the neural network is initially
set, and the mutations take place in the weights of the connections [22].
However, using the Neuroevolution of Augmenting Topologies (NEAT) algorithm,
we are able to evolve the structure of the network as well [23]. This is
useful since it allows us to begin with a simple neural network that evolves
in complexity according to the requirements of the task.

Thus, we evaluate NEAT as a candidate for evolving emergent behavior in robot
swarms. The method we propose is using a homogeneous swarm (where each agent
is using the same neural network) to act in an experiment. At the end of the
experiment, the ‘performance’ of the swarm is evaluated with some user-defined
metric, and we use this as the fitness function for the neural network. This
allows us to use NEAT to evolve a population of neural networks to better
exhibit a desired emergent behavior. A trade off with this is that after
running a simulation, we are left with a single fitness score for the network,
in contrast to the many times we actually call the network.

However, this seems to be an issue inherent to our problem, since the
performance of a swarm can only be evaluated when considering the behavior of
the swarm as a whole. Thus, the actions taken by each individual agent in a
swarm are difficult to assign a reward to. This makes alternate learning
algorithms like Reinforcement Learning (RL) difficult to use, since they
require rewards to be assigned to every action taken. In contrast, this is
coherent with the evolutionary structure we use in our algorithm, since each
‘evaluation’ only needs to yield one fitness value.

###  I-D Organization

The rest of the paper is organized as follows: Section II discusses related
work. Section III details our proposed learning algorithm, as well as common
sensing and output schemes we use. We present our experiments in Section IV,
and discuss the results in Section V. The conclusions we draw are in Section
VI.

##  II Related Work

Much research has focused on the goal of learning individual behaviors for a
robot swarm [2, 5, 6, 3, 14, 16]. Behjat et al. [2] explore how to learn
tactical swarm behavior through a combination of various techniques. These
techniques include learning neural network based robot policies, dynamically
organizing the swarm into groups, and Pareto filtering of points of interest
to reduce the problem dimensionality.

Fan et al. [5] compare different swarm intelligence algorithms against each
other to evaluate their relative performance in obstacle performance under
different circumstances. The algorithms they evaluate are the bat algorithm
(BA), particle swarm optimization (PSO), and the grey wolf optimizer (GWO).
They find that PSO outperforms BA, which outperforms GWO in general. However,
GWO performs better than the other two algorithms in the case of large swarms
and large communication ranges.

Dorigo et al. [6] explore how to create communication protocols between agents
in a swarm that best help methods such as deep reinforcement learning create
good decentralized control policies .

In addition to the studies by Oliveri and Pagello which focus on the learning
aspect (described in Section I-B), there has also been research on helping the
stability of emergent behaviors in swarms. In 2022, Chen and Ng explore secure
communications between robots in a swarm [3]. They model these communications
as a series of random graphs, and use a method involving hash chains to
identify ‘rogue robots’ with high probability. This allows them to ensure that
the emergent behavior exhibited by the swarm is protected. In their paper,
they also create a system to distinguish different classes of robot swarms by
identifying differences in the robot homogeneity, the interactions between
robots, and the interactions with a central control.

Our work is closely related to research done by Trianni et al. in 2003 [26].
They use an evolutionary algorithm to evolve an ‘aggregation’ behavior in a
swarm of robots, inspired by the self-organized aggregation behavior observed
in the cellular slime mold Dictyostelium discoideum [27]. Trianni’s
evolutionary algorithm uses a neural network with no hidden layers, and
evolves by mutating the weights. They use a fitness function that seeks to
minimize the mean distance of each agent from the swarm’s center of mass.
Later research by Bahceci explores this further by varying the parameters of
the evolutionary algorithm to see how it affects the evolution of aggregation
behaviors [1].

We expand on Trianni’s idea by generalizing the algorithm to allow user-
defined fitness functions. Additionally, we use the NEAT algorithm as opposed
to simply evolving the weights of a set network. These extensions allow our
algorithm to theoretically evolve a network to arbitrary complexity in order
to create a desired emergent behavior.

We use the GT-MABs as the robots for our experiments [25]. The control system
for their actuation was designed based on the dynamic model created in [24,
4].

##  III Methodology

###  III-A Evolutionary Algorithm

Algorithm 1 Swarm Evolutionary Algorithm

Input $E$: Experiment setup
$F$: Swarm evaluation function

1:$G\leftarrow$ Population of randomly initialized neural networks

2:loop For desired # of generations

3: for network $g$ in $G$ do

4: Initialize robot swarm experiment $\triangleright$ Using $E$

5: Give each robot $g$ for control

6: Run experiment to termination

7: $f_{g}\leftarrow$ swarm fitness $\triangleright$ evaluated with $F$

8: end for

9: $G^{\prime}\leftarrow$ NEAT reproduction on $G$ $\triangleright$ Fitness
$f_{g}:g\in G$

10: $G\leftarrow G^{\prime}$ $\triangleright$ Update next generation

11:end loop$\triangleright$ Note: parallelization is possible since each

12: genome evaluation is independent

In Algorithm 1, we apply the NEAT algorithm to evolving a swarm of robots to
exhibit some desired emergent behavior. The experiment setup $E$ includes the
objects/agents involved, as well as the definition and implementation of the
network inputs and outputs. In our case, we define $E$ through a CoppeliaSim
simulation [18]. The evaluation function $F$ evaluates the overall behavior of
the swarm throughout the experiment and returns a real number representing the
‘fitness’ of the swarm. $F$ can either directly evaluate the desired emergent
behavior, or evaluate some aspect of the swarm correlated with it. To
implement this, we use the NEAT-Python package [10].

###  III-B Network Inputs

The inputs into the neural network vary for each task, but usually include the
following method of sensing neighbors.

Figure 1: $k$-tant distance sensing with $k=8$

We consider the polar coordinates of each neighbor and split the angular
($\theta$) parameter into $k$ regions. We allow each agent to sense the
inverse distance to the nearest robot on each k-tant, as shown in Fig. 1. We
use the inverse distance since having the sensor read 0 when not detecting
neighbors made more sense with this mapping. It is continuous in the sense
that the limit of the sensor output as the closest neighbor gets arbitrarily
far away is the same as the value when not sensing a neighbor.

Figure 2: $k$-tant neighbor sensing with $k=8$

We also experiment with allowing each agent to sense the number of neighbors
within a range on each $k$-tant instead, as shown in Fig. 2. We introduce this
variant in certain tasks to explore if it is more helpful for an agent to know
the proximity or the count of its neighbors in a given direction.

For 3D cases, we can extend this sensing method into ‘$\ell$-$k$-tant’
sensing. We consider the spherical coordinates of each neighbor and split the
$\phi$ and $\theta$ parameters into $\ell$ parts and $k$ parts respectively,
for a total input dimension of $\ell\cdot k$. However, for our experiments we
only use $k$-tant sensing.

###  III-C Network Outputs/Agent Actuation

Figure 3: GT-MAB (left) and Anki Vector (right) actuation scheme

To actuate the GT-MABs, we use a velocity controller that takes input in
$\mathbb{R}^{3}$. Thus, we must transform the output of the neural network
into a velocity vector, as shown in Fig. 3. For our experiments, we want the
GT-MABs to keep a certain height, so we directly code the $z$ component to
keep the blimp at 1 meter. We set the output dimension of the evolved neural
networks to 2, and assign one output to the $x$ component and one to the $y$
component. Additionally, since the neural network uses a sigmoid activation,
the output is on the interval $[0,1]$. We transform this to allow a command
vector on $[-1,1]^{3}$.

We also run experiments using a simulation of the Anki Vector robots, which
use a wheel velocity controller. We believe this is a significant variation
from a global velocity controller since the motion of the agent is dependent
on the heading angle as well as the network output. To go in a particular
direction, an agent must rotate to face that direction before moving forward.
Thus, we use the Anki Vector experiments to evaluate if our algorithm is
robust to variations in the actuation scheme.

As shown in Fig. 3, we use similar transformations to the GT-MABs to result in
a command vector in $[-1,1]^{2}$. We use this to actuate the Anki Vectors by
assigning one dimension to the left wheel velocity and the other to the right
wheel velocity.

Figure 4: CoppeliaSim model of Anki Vector and GT-MAB

##  IV Experiments

To test our method, we create various tasks for the robot swarms to complete.
We aim to make the tasks achievable through a swarm behavior. Since each
agent’s neural network is only given local information, the algorithm must
evolve an individual policy that produces emergent behavior for the swarm to
perform well.

For each experiment, we train our evolutionary algorithm for 50 generations in
a CoppeliaSim physics simulation. The agents we use are simulated versions of
the GT-MAB and the Anki Vector robots (Fig. 4). We compare the performance of
the evolved policy in each experiment to the performance of a designed policy
that we hard-code to solve the task.

We use ROS2 Foxy to handle transferring sensing and actuation messages between
CoppeliaSim and Python [9]. Our full implementation is available on Github
[17].

###  IV-A Task: Area Coverage

In this task, we simulate a ‘search and rescue’ scenario. The environment we
use is an empty square arena with the agents spawned randomly near the center
(Fig. 6 & Fig. 7). For the fitness function, we adapt deployment entropy, a
measure of how well distributed the agents become in the environment [7].
Deployment entropy is defined by discretizing the environment into a grid, and
measuring the entropy of the distribution of agents in that grid, as shown in
Fig. 5. Since entropy is maximized by a uniform distribution, we believe using
this as a fitness function would encourage the agents to spread out, which is
the desired behavior for a ‘search and rescue’ mission. We implement this task
for the both the Anki Vector robots and the GT-MABs.

Figure 5: Calculation of deployment entropy

####  IV-A1 GT-MABs

Figure 6: GT-MAB Area Coverage experiment initial position

For the GT-MABs, we use 20 agents and divide the environment into a grid with
16 units for calculating deployment entropy. We terminate each experiment
after 60 seconds. We use a variation of Octant Distance Sensing222$k$-tant
distance sensing with $k=8$ shown in Fig. 1 as the input. We additionally
allow the agents the ability to sense the border walls, implemented by
counting the closest wall point on each octant as an additional agent. We add
this capability since the agents would otherwise have no way of sensing the
boundaries. For output, we use our GT-MAB actuation scheme (Fig. 3).

We define a designed policy where each agent moves away from the closest
neighbors. We expect this to cause the agents to distribute themselves evenly,
producing the desired swarm behavior.

####  IV-A2 Anki Vectors

Figure 7: Anki Area Coverage experiment initial position

For the Anki Vectors, we use 10 agents and divide the environment into a grid
with 9 units for calculating deployment entropy. As before, we terminate each
experiment after 60 seconds.

We use Octant Distance Sensing (Fig. 1) along with the Anki’s proximity sensor
as the input. In this case, since the proximity sensor is able to sense the
boundary walls, we do not allow the agent to sense the border with Octant
Distance Sensing. For output, we use the Anki’s wheel velocity actuation
scheme (Fig. 3).

We implement a similar designed policy for each Anki Vector to move away from
its closest neighbors. We notice that this policy is more difficult to design,
since we must include logic to turn the agent to face a more desirable
direction. This supports our assumption that the control scheme of the Anki
Vector robots is different enough from the GT-MAB’s to evaluate our
algorithm’s robustness to control schemes.

###  IV-B Task: Wall Climb

Figure 8: Wall Climb experiment initial position

In previous research we established that a swarm of GT-MABs was able to climb
a wall that was much taller than the altitude they were assigned to hold [20].
Upon closer inspection, we realized that that collisions were resulting in the
agents stacking on top of one another. Since they read their altitude using
ultrasonic sensors, the top GT-MAB would register the other one as a floor,
effectively adding their desired heights. A stack of three agents would allow
the top one to pass over the wall.

We also note that the ‘wall climbing’ behavior is susceptible to the sensing
angle of the GT-MAB’s ultrasound-based range sensor. With a large sensing
angle, the GT-MABs are able to sense the wall itself when they are next to it,
which results in a single agent being able to climb the wall on its own. To
avoid this, we use a narrow ultrasound angle so that the agents must stack in
order to climb the wall.

####  IV-B1 Octant Distance Sensing

We design the Wall Climb task to replicate this scenario. The environment we
use is an arena with a 3 meter tall wall on the $y$ axis. The GT-MABs spawn
randomly on the right side, as shown in Fig. 8. For the fitness function, we
use the number of agents that end the experiment on the other side of the
wall. If no agents succeed, we subtract a penalty of the distance of the
closest blimp to the wall. We believe using this as a fitness function will
encourage the “aggregation” behavior, since this is the only way the agents
can climb the wall. We use 20 agents, and terminate each experiment after 60
seconds.

We use Octant Distance Sensing as our input (Fig. 1) and our GT-MAB actuation
scheme as output (Fig. 3).

We define a designed policy where each agent moves towards their closest
neighbor, as well as slightly towards the direction of the wall. We expect
this to cause the GT-MABs to flock and stack on top of one another. Once they
are stacked, we predict they will slowly make their way over the wall,
producing the desired swarm behavior.

####  IV-B2 Octant Neighbor Sensing

Figure 9: Surround Target experiment initial position

We also compare this with using the Octant Neighbor Sensing scheme as our
input, as shown in Fig. 2. Our designed policy for this scheme is similar,
where we move towards the direction with the highest neighbor count. We
predict that sensing the number of neighbors is more useful to create the
‘flocking’ behavior we desire. Thus, we expect that this sensing scheme will
be more successful.

###  IV-C Task: Surround Target

In this task, we simulate a goal of securing a perimeter around a target [21].
The environment we use is an arena containing a cylinder target spawned
randomly, with the GT-MABs also spawned randomly in the arena.

For the fitness function, we consider all the agents close to the target.
Using polar coordinates centered at the target, we find the measure of the
largest uncovered arc $a$. We set the fitness to $-a$ since we want to
minimize this value. We also subtract a penalty of the distance of the closest
agent to the target for experiments where no agents are near the target. We
use 6 agents and terminate each experiment after 60 seconds.

We include Octant Distance Sensing in our input (Fig. 1), as well as the
agent’s displacement vector from the target. For output, we use our GT-MAB
actuation scheme (Fig. 3).

We define a designed policy where each agent moves towards the target, as well
as slightly away from their closest neighbor. We expect this to cause the GT-
MABs to head for the target, then slowly spread out, producing the desired
swarm behavior of surrounding the target.

##  V Results

TABLE I: Comparison of Fitness in Evolveda and Designed Behaviors Experiment | Behavior | Mean | Stdev.
---|---|---|---
GT-MAB Area Coverage | Evolved | 2.53 | 0.10
Designed | 2.62 | 0.12
Anki Area Coverage | Evolved | 2.09 | 0.083
Designed | 1.87 | 0.20
Wall Climb | Evolved | 16.70 | 1.12
(Distance Sense) | Designed | 12.13 | 3.78
Wall Climb | Evolved | 16.72 | 1.32
(Neighbor Sense) | Designed | 15.33 | 1.58
Surround Target | Evolved | -2.21 | .70
Designed | -1.24 | 0.14
a: Using the best genome from the final generation.

Figure 10: GT-MAB Area Coverage evolved behavior (left) and designed policy
behavior (right)

For each experiment, after training for 50 generations, we compare the fitness
of the best genome in the final generation with the fitness of our designed
policy. We run each for 60 trials in their respective experiments and collect
the fitness mean and standard deviation. We report our results in Table I.

###  V-A Area Coverage

####  V-A1 GT-MABs

Figure 11: GT-MAB Area Coverage population fitness across generations

After training the GT-MAB models in CoppeliaSim, we observe that the behavior
they evolved seems to be for each agent to move in a direction away from its
close neighbors. Fig. 10 shows an experiment run with the network with the
best genome in the last generation of training. The result of this behavior
seems to be a well distributed swarm, with the agents arranged very close to a
lattice. From Fig. 11, we can see that after around the 15th generation, the
best genome of each generation achieved close to the theoretic maximal
entropy3332.718 from placing one agent in 12 of the grid units, and two agents
in 4 of the units.

The result of the designed policy, shown in Fig. 10, also had the swarm
distribute in the environment, forming a similar lattice. Comparing these
results, we can see that the designed policy outperforms the evolved behavior
by about one standard deviation (Table I).

Figure 12: Anki Area Coverage evolved behavior (left) and designed policy
behavior (right)

Overall, we show that in this task, our algorithm learns a local behavior that
closely approximated the desired ‘search and rescue’ emergent swarm behavior.
From the slight upwards trend in the mean of Fig. 11, we expect that the
evolved results may perform closer to the designed policy with more training.

####  V-A2 Anki Vectors

Figure 13: Anki Area Coverage population fitness across generations

After training the Anki Vector models in CoppeliaSim, we observe the behavior
appears to be for each agent to spin in place, and move away from close
neighbors directly in front of or behind it. Fig. 12 shows an experiment run
with the best genome in the last generation of training. The result of this
behavior appears to be a well distributed swarm. From Fig. 13, we can see that
after a few generations, the best genome of each generation achieved the
theoretic maximal entropy4442.16 from placing one agent in 8 of the grid
units, and two agents in 1 of the units.

The result of the designed policy, shown in Fig. 12, had the swarm distribute
reasonably well in the environment. Comparing these results, we can see that
the evolved behavior outperforms the designed policy, since the mean is about
one standard deviation higher (Table I). This seems to be due to the designed
policy favoring sending agents to the edges of the environment.

Overall, we show that in this task, our algorithm learns a local behavior that
closely approximated the desired ‘search and rescue’ emergent swarm behavior.
The evolved behavior greatly outperforms the designed policy, which shows that
our evolutionary algorithm performs comparably to designing the behavior with
knowledge of the task.

###  V-B Wall Climb

####  V-B1 Octant Distance Sensing

Figure 14: Distance Sensing Wall Climb evolved behavior (left) and designed
policy behavior (right)

Figure 15: Distance Sensing Wall Climb (left) and Neighbor Sensing Wall Climb
(right) population fitness across generations

After training the GT-MAB models in CoppeliaSim, we observe that the behavior
they evolved seems to be for each agent to move in a direction towards its
closest neighbors, in addition to moving in the direction to climb the wall.
Fig. 14 shows an experiment run with the network with the best genome in the
last generation of training. The result of this behavior does seem to be a
flock of agents stacked on top of each other, climbing the wall. From Fig. 15,
we can see that after a few generations, the best genome of each generation
achieved having about 18 GT-MABs make it over the wall.

####  V-B2 Octant Neighbor Sensing

Figure 16: Neighbor Sensing Wall Climb evolved behavior (left) and designed
policy behavior (right)

As for the Neighbor Sensing experiment, we noticed both sensing methods
achieved similar population fitness values across generations (Fig. 15). The
flocking behavior evolved was also similar in both experiments (Fig. 14, Fig.
16). From Table I, we see there is no statistically significant difference
between the two fitness results555We perform a two sample $t$-test with
unequal variance on the fitnesses obtained from the evolved Neighbor Sensing
and Distance Sensing experiments. Using their respective means and standard
deviations of (16.70,1.12) and (16.72,1.32) with a sample size of 60 for both,
we arrive at a $p$-value of $p>.9$, which is much larger than $.05$, the
accepted value for statistical significance. This shows that with our sample
size, there is no significant difference between the fitnesses obtained from
the different sensing methods..

The result of the designed policy in both of these experiments similarly
causes the agents to flock together and climb the wall. However, in this case
we do notice a difference in the two modes of sensing. In the Distance Sensing
experiments, the designed policy seems to be more likely to cause the agents
to form several smaller groups as opposed to one large one. This results in a
lower fitness due to more agents being left behind, as shown in Table I.
Comparing this with our evolved behavior, we can see in both experiments that
the evolved behavior outperforms the best designed policy by about one
standard deviation.

Overall, we show that in this task, our algorithm learns a local behavior that
closely approximated the desired ‘flocking’ emergent swarm behavior. We can
further see that although the different sensing modes has an effect on our
designed policy, the evolutionary algorithm is robust to these variations.

###  V-C Surrounding a Target

Figure 17: Surround Target evolved behavior (left) and designed policy
behavior (right)

Figure 18: Surround Target population fitness across generations

After training the GT-MAB models in CoppeliaSim, we observe that the behavior
they evolved seems to be for each agent to spiral around the target, in
addition to moving away from its closest neighbor. Fig. 17 shows an experiment
run with the network with the best genome in the last generation of training.
The result of this behavior does seem to be a flock of agents surrounding the
target. From Fig. 18, we can see that after a few generations, the best genome
of each generation achieves having a largest uncovered arc of about 2 radians.

The result of the designed policy, shown in Fig. 17, also caused the agents to
surround the target. Comparing these results, we can see that our designed
policy outperforms the evolved behavior by at least one standard deviation
(Table I).

Overall, we show that in this task, our algorithm learns a local behavior that
closely approximates the desired emergent swarm behavior. However, our
designed policy still does outperform the evolved behavior, which suggests
that for some tasks it is more difficult to evolve the optimal behavior. We
believe fine tuning of the evolution parameters may help this issue, but that
is outside the scope of our research.

##  VI Conclusion

In this paper, we present a novel extension of the NEAT algorithm designed to
learn emergent behaviors in robot swarms. The algorithm we present is able to
be applied on robot swarms with various modes of sensing and actuation.
Results from simulations show that individual agent behaviors evolved using
this method are comparable to hand designed policies at producing desired
complex emergent behaviors.

In future research, we plan to test our evolved policies on the physical GT-
MABs and Anki Vector robots. We also plan to evaluate our algorithm on a more
complex set of tasks. We also may explore the fine tuning of NEAT parameters
to try improving our results. Additionally, the accuracy of the physics
evaluation seems to be dependant on the performance of the computer, as this
determines the speed that each agent can respond to stimuli. In this paper, we
choose to speed up our experiments by running multiple simulations in
parallel. However, using less parallelization or using a more powerful
computer could improve the results of our algorithm.

## References

  * [1] Erkin Bahceci and Erol Sahin.  Evolving aggregation behaviors for swarm robotic systems: a systematic case study.  In Proceedings 2005 IEEE Swarm Intelligence Symposium, 2005. SIS 2005., pages 333 – 340, 07 2005.
  * [2] Amir Behjat, Hemanth Manjunatha, Prajit Krisshna Kumar, Apurv Jani, Leighton Collins, Payam Ghassemi, Joseph Distefano, David Doermann, Karthik Dantu, Ehsan Esfahani, and Souma Chowdhury.  Learning robot swarm tactics over complex adversarial environments.  In 2021 International Symposium on Multi-Robot and Multi-Agent Systems (MRS), pages 83–91, Nov 2021.
  * [3] Liqun Chen and Siaw-Lynn Ng.  Securing emergent behaviour in swarm robotics.  Journal of Information Security and Applications, 64:103047, 2022\.
  * [4] Sungjin Cho, Vivek Mishra, Qiuyang Tao, Paul Vamell, Matt King-Smith, Aneri Muni, Weston Smallwood, and Fumin Zhang.  Autopilot design for a class of miniature autonomous blimps.  In 2017 IEEE Conference on Control Technology and Applications (CCTA), pages 841–846, 2017.
  * [5] Jiaqi Fan, Mengqi Hu, Xianghua Chu, and Dong Yang.  A comparison analysis of swarm intelligence algorithms for robot swarm learning.  In 2017 Winter Simulation Conference (WSC), pages 3042–3053, Dec 2017.
  * [6] Maximilian Hüttenrauch, Adrian Šošić, and Gerhard Neumann.  Local communication protocols for learning complex swarm behaviors with deep reinforcement learning.  In Swarm Intelligence, pages 71–83, Cham, 2018. Springer International Publishing.
  * [7] John D. Kelly, Daniel M. Lofaro, and Donald Sofge.  Persistent area coverage for swarms utilizing deployment entropy with potential fields.  In 2020 17th International Conference on Ubiquitous Robots (UR), pages 479–486, 2020.
  * [8] Thalia May Laing, Siaw-Lynn Ng, Allan Tomlinson, and Keith M Martin.  Security in swarm robotics.  Handbook of Research on Design, Control, and Modeling of Swarm Robotics, page 42–66, 2016.
  * [9] Steven Macenski, Tully Foote, Brian Gerkey, Chris Lalancette, and William Woodall.  Robot operating system 2: Design, architecture, and uses in the wild.  Science Robotics, 7(66):eabm6074, 2022.
  * [10] Alan McIntyre, Matt Kallada, Cesar G. Miguel, Carolina Feher de Silva, and Marcio Lobo Netto.  neat-python.  https://github.com/CodeReclaimers/neat-python, 2008.
  * [11] Nathan J. Mlot, Craig A. Tovey, and David L. Hu.  Fire ants self-assemble into waterproof rafts to survive floods.  Proceedings of the National Academy of Sciences, 108(19):7669–7673, 2011.
  * [12] Iñaki Navarro and Fernando Matia.  An introduction to swarm robotics.  ISRN Robotics, 2013, 01 2013.
  * [13] Richard E. Neapolitan and Xia Jiang.  Chapter 5 - decision analysis fundamentals.  In Probabilistic Methods for Financial and Marketing Informatics, pages 177–228. Morgan Kaufmann, Burlington, 2007.
  * [14] Giorgio Oliveri, Lucas C. van Laake, Cesare Carissimo, Clara Miette, and Johannes T. B. Overvelde.  Continuous learning of emergent behavior in robotic matter.  Proceedings of the National Academy of Sciences, 118(21):e2017015118, 2021.
  * [15] Timothy O’Connor.  Emergent Properties.  In The Stanford Encyclopedia of Philosophy. Metaphysics Research Lab, Stanford University, Winter 2021 edition, 2021.
  * [16] E. Pagello, A. D’Angelo, C. Ferrari, R. Polesel, R. Rosati, and A. Speranzon.  Emergent behaviors of a robot team performing cooperative tasks.  Advanced Robotics, 17(1):3–19, 2003.
  * [17] Pranav Rajbhandari.  Swarm Coppeliasim.  https://github.com/pranavraj575/swarm_coppeliasim, July 2023.
  * [18] E. Rohmer, S. P. N. Singh, and M. Freese.  V-rep: A versatile and scalable robot simulation framework.  In 2013 IEEE/RSJ International Conference on Intelligent Robots and Systems, pages 1321–1326, 2013.
  * [19] Destin Sandlin.  Gregarious caterpillar locomotion - raw video for analysis.  https://www.youtube.com/watch?v=YehR0wSUioY, Jun 2013.
  * [20] Tristan Schuler, Cameron Kabacinski, Daniel M. Lofaro, Dhawal Bhanderi, Jennifer Nguyen, and Donald Sofge.  Wall climbing emergent behavior in a swarm of real-world miniature autonomous blimps.  In Proceedings of the 15th International Conference on Agents and Artificial Intelligence, ICAART 2023, Volume 1, Lisbon, Portugal, February 22-24, 2023, pages 225–232. SCITEPRESS, 2023.
  * [21] Daigo Shishika and Vijay Kumar.  A review of multi agent perimeter defense games.  In Decision and Game Theory for Security, pages 472–485, Cham, 2020\. Springer International Publishing.
  * [22] Kenneth Stanley, Jeff Clune, Joel Lehman, and Risto Miikkulainen.  Designing neural networks through neuroevolution.  Nature Machine Intelligence, 1, 01 2019.
  * [23] Kenneth O. Stanley and Risto Miikkulainen.  Evolving neural networks through augmenting topologies.  Evolutionary Computation, 10(2):99–127, 2002.
  * [24] Qiuyang Tao, Jaeseok Cha, Mengxue Hou, and Fumin Zhang.  Parameter identification of blimp dynamics through swinging motion.  In 2018 15th International Conference on Control, Automation, Robotics and Vision (ICARCV), pages 1186–1191, 2018.
  * [25] Qiuyang Tao, Junkai Wang, Zheyuan Xu, Tony X. Lin, Ye Yuan, and Fumin Zhang.  Swing-reducing flight control system for an underactuated indoor miniature autonomous blimp.  IEEE/ASME Transactions on Mechatronics, 26(4):1895–1904, 2021.
  * [26] Vito Trianni, Roderich Groß, Thomas H. Labella, Erol Şahin, and Marco Dorigo.  Evolving aggregation behaviors in a swarm of robots.  In Advances in Artificial Life, pages 865–874, Berlin, Heidelberg, 2003. Springer Berlin Heidelberg.
  * [27] Vanoss, C and Panfilov, Alexander and Hogeweg, P and Siegert, F and Weijer, CJ.  Spatial pattern formation during aggregation of the slime mould Dictyostelium discoideum.  JOURNAL OF THEORETICAL BIOLOGY, 181(3):203–213, 1996\.