# Scalable Diffusion Models with State Space Backbone

Zhengcong Fei Kunlun Inc. Mingyuan Fan Kunlun Inc. Changqian Yu
Junshi Huang*
Kunlun Inc<EMAIL_ADDRESS>
###### Abstract

This paper presents a new exploration into a category of diffusion models
built upon state space architecture. We endeavor to train diffusion models for
image data, wherein the traditional U-Net backbone is supplanted by a state
space backbone, functioning on raw patches or latent space. Given its notable
efficacy in accommodating long-range dependencies, Diffusion State Space
Models (DiS) are distinguished by treating all inputs including time,
condition, and noisy image patches as tokens. Our assessment of DiS
encompasses both unconditional and class-conditional image generation
scenarios, revealing that DiS exhibits comparable, if not superior,
performance to CNN-based or Transformer-based U-Net architectures of
commensurate size. Furthermore, we analyze the scalability of DiS, gauged by
the forward pass complexity quantified in Gflops. DiS models with higher
Gflops, achieved through augmentation of depth/width or augmentation of input
tokens, consistently demonstrate lower FID. In addition to demonstrating
commendable scalability characteristics, DiS-H/2 models in latent space
achieve performance levels akin to prior diffusion models on class-conditional
ImageNet benchmarks at the resolution of 256$\times$256 and 512$\times$512,
while significantly reducing the computational burden. The code and models are
available at: https://github.com/feizc/DiS.

##  1 Introduction

Diffusion models [26, 60, 64, 4] have emerged as potent deep generative models
in recent years [9, 27, 54], for their capacity in high-quality image
generation. Their rapid development has led to widespread application across
various domains, including text-to-image generation [53, 54, 55, 21], image-
to-image generation [6, 70, 69], video generation [25, 45, 45], speech
synthesis [5, 36], and 3D synthesis [52, 71]. Along with the development of
sampling algorithms [10, 32, 41, 42, 62, 66], the revolution of backbones
stands as a pivotal aspect in advancement of diffusion models. A
representative example is U-Net based on a convolutional neural network (CNN)
[26, 63], which has been prominently featured in previous research. The CNN-
based UNet is characterized by a group of down-sampling blocks, a group of up-
sampling blocks, and long skip connections between the two groups [9, 53, 54,
55]. Similarly, Transformer-based architecture replaces sampling block with
self-attention while keeping the remain untouched [51, 67, 1], resulting in
streamlined yet effective performance.

On the other hand, state space models (SSMs) with efficient hardware-aware
designs, have shown great potential in the realm of long sequence modeling
[29, 19, 20, 22, 49, 59]. As self-attention mechanism in Transformer scales
quadratically with the input size, making them resource-intensive when dealing
with long-range visual dependencies, _i.e._ , high resolution images. Recent
efforts, exampified by the work on Mamba [16, 15], have sought to address by
integrating time-varying parameters into the SSM and proposing a hardware-
aware algorithm to enable efficient training and inference. The commendable
scaling performance of Mamba underscores its promise as a viable avenue for
constructing efficient and versatile backbones within the domain of SSMs.
Motivated by the successes observed in language modeling with Mamba, a
pertinent inquiry arises: whether we can build SSM-based U-Net in diffusion
models?

Figure 1: The proposed state space-based diffusion models. It treats all
inputs including the time, condition and noisy image patches as tokens and
employs skip connections between shallow and deep layers. Different from
original Mamba for text sequence modeling, our SSM block process the hidden
states sequence with both forward and backward directions.

In this paper, we aim to design a simple and general state space-based
architecture for diffusion models, referred to as DiS. Following the design
principles of [51], DiS treats all inputs including time, condition and noisy
image patches as discrete tokens. Crucially, we conduct a systematical
ablation study encompassing the incorporation of conditioning factors and the
optimization of model architecture across all components. Noteworthy is DiS’s
adherence to the established best practices of state space models, renowned
for their superior scalability in image generation tasks when compared to CNN
or Transformers, all while maintaining lower computational overhead.

Experimentally, we assess the performance DiS across both unconditional and
class-conditional image generation tasks. In all settings, DiS demonstrate
comparative, if not superior, efficay when juxtaposed with CNN-based or
Transformer-based U-Nets of a similar size. Additionally, we provide evidence
showcasing that DiSs are scalable architectures for diffusion models, where a
discernible correlation is observed between network complexity and sample
quality. Moreover, experiments yield impressive results, with DiS achieving a
comparable FID scores in class-conditional image generation conducted on
ImageNet at a resolution of 256$\times$256 and 512$\times$512\. This study
endeavors to elucidate the criticality of architectural selections within the
domain of diffusion models, concurrently providing empirical benchmarks that
can inform forthcoming research in generative modeling. Our aspiration is that
the insights gleaned from DiS will serve as foundational knowledge for future
investigations concerning backbone architectures in diffusion models,
ultimately enriching the landscape of generative modeling, particularly in the
context of large-scale cross-modality datasets.

##  2 Methodology

###  2.1 Preliminaries

####  2.1.1 Diffusion models.

Prior to presenting our architecture, we provide a concise overview of
fundamental concepts relevant to our discussion. Diffusion models [26, 60]
gradually inject noise to data, and subsequently reverse this process to
generate data from noise. The noise injection process, also referred to as the
forward process, can be formalized as Markov chain as:

$\begin{split}q(\textbf{x}_{1:T}|\textbf{x}_{0})=\prod_{t=1}^{T}q(\textbf{x}_{t}|\textbf{x}_{t-1}),\\\
q(\textbf{x}_{t}|\textbf{x}_{t-1})=\mathcal{N}(\textbf{x}_{t}|\sqrt{\alpha_{t}}\textbf{x}_{t-1},\beta_{t}\textbf{I}),\end{split}$
(1)

Here $\textbf{x}_{0}$ is the data, $\alpha_{t}$ and $\beta_{t}$ represent the
noise schedule, ensuring $\alpha_{t}+\beta_{t}=1$. To reverse this process, a
Gaussian model
$p(\textbf{x}_{t-1}|x_{t})=\mathcal{N}(\textbf{x}_{t-1}|\mu_{t}(\textbf{x}_{t}),\sigma_{t}^{2}\textbf{I})$
is employed to approximate the ground truth reverse transition
$q(\textbf{x}_{t-1}|\textbf{x}_{t})$, where the learning is tantamount to a
noise prediction task. Formally, a noise prediction network
$\epsilon_{\theta}(\textbf{x}_{t},t)$ is incorporated by minimizing a noise
prediction objective, _i.e._ ,
$\min_{\theta}\mathbb{E}_{t,\textbf{x}_{0},\epsilon}||\epsilon-\epsilon_{\theta}(\textbf{x}_{t},t)||_{2}^{2}$,
where $t$ is uniformly distributed between 1 and $T$. To learn conditional
diffusion models, _e.g._ , class-conditional [9] or text-to-image [53, 2]
models, additional condition information is integrated into the noise
prediction objective as:

$\min_{\theta}\mathbb{E}_{t,\textbf{x}_{0},\textbf{c},\epsilon}||\eta-\eta_{\theta}(\textbf{x}_{t},t,\textbf{c})||^{2}_{2},$
(2)

where c is the condition index or its continuous embedding.

####  2.1.2 State space backbone.

State space models are conventionally defined as linear time-invariant systems
that builds a map stimulation $x(t)\in\mathbb{R}^{N}$ to response
$y(t)\in\mathbb{R}^{N}$ by a latent state $h(t)\in\mathbb{R}^{N}$. The process
can be formulated as:

$\displaystyle h^{\prime}(t)$ $\displaystyle=\textbf{A}h(t)+\textbf{B}x(t),$
(3) $\displaystyle y(t)$ $\displaystyle=\textbf{C}h(t),$

where $\textbf{A}\in\mathbb{R}^{N\times N}$ and
$\textbf{B},\textbf{C}\in\mathbb{R}^{N}$ denote the state matrix, input
matrix, and output matrix, respectively. In the quest to derive the output
sequence $y(t)$ at time $t$, the analytical solution for obtaining $h(t)$
proves to be a formidable challenge. Conversely, real-world data typically
manifests in discrete rather than continuous forms. As an alternative, we can
discretize the system in Equation 3 as follows:

$\begin{split}h_{t}&=\overline{\textbf{A}}h_{t-1}+\overline{\textbf{B}}x_{t},\\\
y_{t}&=\textbf{C}h_{t},\end{split}$ (4)

where $\overline{\textbf{A}}:=\exp(\Delta\cdot\textbf{A})$ and
$\overline{\textbf{B}}:=(\Delta\cdot\textbf{A})^{-1}(\exp(\Delta\cdot\textbf{A})-I)\cdot\Delta\textbf{B}$
are the discretized state parameters and $\Delta$ is the discretization step
size.

While SSMs boast rich theoretical properties, they are often afflicted by high
computational costs and numerical instability. These shortcomings have
prompted the development of structured state space sequence models (S4), which
seek to mitigate these challenges by imposing a structured format on the state
matrix A through the utilization of HIPPO matrices. This structural
enhancement has led to notable advancements in both performance and
efficiency. Notably, S4 has exhibited a substantial performance advantage over
Transformers, which necessitates effective modeling of long-range dependencies
[19]. More recently, Mamba [16] further advance its potential through an
input-dependent selection mechanism and a faster hardware-aware algorithm.

###  2.2 Model Structure Design

We introduce Diffusion State Space Models, denoted as DiS, a simple and
general architecture for diffusion models in image generation. Specifically,
DiS parameterizes the noise prediction network
$\epsilon_{\theta}(\textbf{x}_{t},t,\textbf{c})$, which takes the timestep
$t$, condition c and noised image $\textbf{x}_{t}$ as inputs and predicts the
noise injected into $\textbf{x}_{t}$. Our objective is to closely adhere to
the advanced state space architecture to preserve its scalability
characteristics, hence, DiS is based on the bidirectional Mamba [16, 72]
architecture which operates on sequences of tokens. Figure 1 provides an
overview of the complete DiS architecture. In this section, we expound upon
the forward pass of DiS, as well as the components of the design space of the
DiS class.

####  2.2.1 Image patchnify.

The first layer of DiS undertakes transformation of input image
$\textbf{I}\in\mathbb{R}^{H\times W\times C}$ into flatten 2-D patches
$\textbf{X}\in\mathbb{R}^{J\times(p^{2}\cdot C)}$. Subsequently, converts it
into a sequence of $J$ tokens, each of dimension $D$, by linearly embedding
each patch in the input. Consistent with [11], we apply learnable positional
embeddings to all input tokens. The number of tokens $J$ created by patchify
is determined by the hyperparameter patch size $p$ to $\frac{H\times
W}{p^{2}}$. The patchnify supports both raw pixel and latent space. We set
$p=2,4,8$ to the design space.

####  2.2.2 SSM block.

After embedding layer, the input tokens are processed by a sequence of SSM
blocks. In addition to noised image inputs, diffusion models sometimes process
additional conditional information such as noise timesteps $t$, condition c
such as class labels or natural language. Given that original Mamba block is
designed for 1-D sequence, we resort to [72], which incorporates bidirectional
sequence modeling tailed for vision tasks. The designs introduce subtle yet
pivotal modifications to the standard SSM block design. As shown in Figure 1
right part, forward of SSM blocks combine both forward and backward
directions.

####  2.2.3 Skip connection.

Given a series of $L$ SSM blocks, we categorize the stack SSM block into first
half $\lfloor\frac{L}{2}\rfloor$ shallow group, one middle layer, and second
half $\lfloor\frac{L}{2}\rfloor$ deep group. Let
$\textbf{h}_{shallow},\textbf{h}_{deep}\in\mathbb{R}^{J\times D}$ be the
hidden states from the main branch and long skip branch respectively. We
consider directly concatenating them and perform a linear projection, _i.e._ ,
Linear(Concat($\textbf{h}_{shallow},\textbf{h}_{deep}$)), before feeding them
to the next SSM block.

####  2.2.4 Linear decoder.

After the final SSM block, we need to decode our sequence of hidden states
into an output noise prediction and diagonal covariance prediction. Both of
these outputs retain a shape identical to the original spatial input. We use a
standard linear decoder, _i.e._ , we apply the final layer norm and linearly
decode each token into a $p^{2}\cdot C$ tensor. Finally, we rearrange the
decoded tokens into their original spatial layout to get the predicted noise
and covariance.

####  2.2.5 Condition incorporation.

To effectively integrate additional conditions, we adopt a straightforward
strategy of appending the vector embeddings of timestep $t$ and condition c as
two supplementary tokens in the input sequence. These tokens are treated
equivalently to the image tokens, akin to the approach with cls tokens in
Vision Transformers [11]. This approach enables the utilization of SSM blocks
without necessitating any modifications. Following the final block, the
conditioning tokens are removed from the sequence. We also explore to adaptive
normalization layer, where replace standard norm layer with adaptive norm
layer. That is, rather than directly learn dimension-wise scale and shift
parameters, we regress them from the sum of the embedding vectors of $t$ and
c, which will be discussed in the Experiments Section.

Table 1: Scaling law model size. The model sizes and detailed hyperparameters settings for scaling experiments.  | #Params | #Blocks $L$ | Hidden dim. $D$ | #Expand $E$ | Gflops
---|---|---|---|---|---
Small | 28.4M | 25 | 384 | 2 | 0.43
Base | 119.1M | 25 | 768 | 2 | 1.86
Medium | 229.4M | 49 | 768 | 2 | 3.70
Large | 404.0M | 49 | 1024 | 2 | 6.57
Huge | 900.6M | 49 | 1536 | 2 | 14.79

###  2.3 Computation Analysis

In summary, the hyper-parameters of our architecture encompass the following:
the number of blocks $L$, hidden state dimension $D$, expanded state dimension
$E$, and SSM dimension $N$. Various configurations of DiS are delineated in
Table 1. They cover a wide range of model sizes and flop allocations, from 28M
to 900M, thus affording comprehensive insights into the scalability
performance. Aligned with [51], Gflop metric is evaluated in 32$\times$32
unconditional image generation with patch size $p=4$ with thop python package.
Following [16] we also set the SSM dimension of all model variants to 16.

Both SSM block in DiS and self-attention in Transformer play a key role in
modeling long context adaptively. We further provide a theoretical analysis
pertaining to computation efficiency. Given a sequence
$\textbf{X}\in\mathbb{R}^{1\times J\times D}$ and the default setting $E=2$,
the computation complexity of a self-attention and SSM operation are
delineated as:

$\displaystyle\mathcal{O}(\text{SA})=4JD^{2}+2J^{2}D,$ (5)
$\displaystyle\mathcal{O}(\text{SSM})=3J(2D)N+J(2D)N^{2},$ (6)

where we can see that self-attention is quadratic to sequence length $J$, and
SSM is linear with respect to sequence length $J$. It is noteworthy that $N$
is a fixed parameter, typically set to 16 by default. This computational
efficiency renders DiS amenable to scalability in scenarios necessitating
generation with large sequence lengths, such as gigapixel applications.

##  3 Experiments

We delve into the design space and scrutinize the scaling properties of our
DiS model class. Our model is named according to their configs and patch size
$p$; for instance, DiS-L/2 refers to the Large version config and $p=2$.

###  3.1 Experimental Settings

####  3.1.1 Datasets.

For unconditional image learning, we consider CIFAR10 [37], which comprises
50K training images, and CelebA 64x64 [40], which contains 162,770 training
images of human faces. For class-conditional image learning, we consider
ImageNet [8] at 256$\times$256 and 512$\times$512 resolutions, which contains
1,281,167 training images across 1,000 different classes. The only data
augmentation is horizontal flips. We train 500K iterations on CIFAR10 and
CelebA 64$\times$64 with a batch size of 128. We train 500K and 1M iterations
on ImageNet 256$\times$256 and 512$\times$512, with a batch size of 1024.

####  3.1.2 Implementation details.

We use the AdamW optimizer [33] without weight decay across all datasets. We
maintain a learning rate of 1e-4 with a cosine schedule. In our early
experiments, we tried learning rates ranging from 1e-4 to 5e-4, and found that
loss may become NAN if the learning rate set is larger than 3e-4. Consistent
with practice in literature [51], we utilize an exponential moving average of
DiS weights over training with a decay of 0.9999. All results were reported
using the EMA model. Our models are trained on Nvidia A100 GPU. When trained
on ImageNet dataset at resolution of 256$\times$256 and 512$\times$512, we
adopt classifier-free guidance [28] following [54] and use a off-the-shelf
pre-trained variational autoencoder (VAE) model [35] from Stable Diffusion
[54] provided in huggingface111https://huggingface.co/stabilityai/sd-vae-ft-
ema. The VAE encoder has a downsampling factor of 8. We retrain diffusion
hyperparameters from [51], using a $t_{max}=1000$ linear variance schedule
ranging from $1\times 10^{-4}$ to $2\times 10^{-2}$ and parameterization of
the covariance.

Figure 2: Ablation experiments with DiS-S/2 on CIFAR10 dataset. We report FID
metrics on 10K generated samples. (a) Patch size. A smaller patch size can
improve the generation performance. (b) Long skip. Combining the long skip
branch can accelerate the training as well as optimize generated results.

####  3.1.3 Evaluation metrics.

We measure image generation performance with Frechet Inception Distance (FID)
[24], a widely adopted metric for assessing the quality of generated images.
We follow convention when comparing against prior works and report FID-50K
using 250 DDPM sampling steps [50] following the process of [9]. Notably, the
FID values presented in this section are calculated without classifier-free
guidance, unless specified otherwise. We additionally report Inception Score
[56], sFID [47] and Precision/Recall [38] as secondary metrics.

###  3.2 Model Analysis

In this section, we perform a systematical empirical investigation into the
fundamental components of DiS models. Specifically, we ablate on the CIFAR10
dataset, evaluate the FID score every 50K training iterations on 10K generated
samples, instead of 50K samples for efficiency identical to [1], and determine
the optimal default implementation details.

####  3.2.1 Effect of patch size.

We train patch size range over (8, 4, 2) with DiS-S configuration on the
CIFAR10 dataset. As depicted in Figure 2(a) shows FID metric fluctuates in
response to patch size decrease when maintaining model size consistent. We
observe that discernible FID enhancements are observed throughout the training
process by augmenting the number of tokens processed by DiS, holding
parameters approximately fixed. Consequently, optimal performance necessitates
a smaller patch size, such as 2. We postulate that this requirement stems from
the inherently low-level nature of the noise prediction task in diffusion
models, which favors smaller patches, in contrast to higher-level tasks like
classification. Moreover, given the computational costs associated with using
small patch sizes for high-resolution images, an alternative approach involves
transforming these images into low-dimensional latent representations,
subsequently also modeled using the DiS series.

Figure 3: Model analysis for different designs. (a) Variants of condition
incorporation. Feed time as token into the networks is effective. (b) Model
parameters scaling. As we expected, holding patch size constant, increasing
the model size can consistently improve the generation performance.

####  3.2.2 Effect of long skip.

For evaluating the efficiency of skipping operation, we examine three variants
including: (i) Concatenation, _i.e._ , Linear(Concat($\textbf{h}_{shallow}$,
$\textbf{h}_{deep}$)); (ii) Addition, _i.e._ ,
$\textbf{h}_{shallow}+\textbf{h}_{deep}$; (iii) No skip connection. As show in
Figure 2(b), directly adding the hidden states from the shallow and deep
layers does not yield any discernible benefits. Given that the SSM block
already incorporates residual skip connections internally, the hidden states
of deep layer inherently encapsulates some information from the shallow layer
in a linear fashion. Consequently, the addition operation in the long skip
merely amplifies the contribution of the shallow layer’s hidden state, without
fundamentally altering the network’s behavior. In contrast, employing
concatenation involves a learnable linear projection on the shallow hidden
states and effectively enhances performance compared to the absence of a long
skip connection.

Table 2: Benchmarking unconditional image generation on CIFAR10. DiS-S/2 model
obtains comparable results with fewer parameters.  Unconditional CIFAR10
---
Model | #Params |  FID$\downarrow$
_Diff. based on U-Net_
DDPM [26] | 36M | 3.17
EDM [30] | 56M | 1.97
_Diff. based on Transformer_ |
GenViT [67] | 11M | 20.20
U-ViT-S/2 [1] | 44M | 3.11
_Diff. based on SSM_ |
DiS-S/2 | 28M | 3.25

Table 3: Benchmarking unconditional image generation on CelebA 64$\times$64.
DiS-S/2 maintains a superior generation performance in small model settings.
Unconditional CelebA 64$\times$64
---
Model | #Params |  FID$\downarrow$
_Diff. based on U-Net_ |
DDIM [61] | 79M | 3.26
Soft Trunc. [31] | 62M | 1.90
_Diff. based on Transformer_ |
U-ViT-S/4 [1] | 44M | 2.87
_Diff. based on SSM_ |
DiS-S/2 | 28M | 2.05

####  3.2.3 Effect of condition combination.

We explore two ways to integrate the conditional timestep $t$ into the
network: (i) treat timestep $t$ as a token and directly concatenate it with
image patches in prefix. (ii) incorporate the time embedding after the layer
normalization in the SSM block, which is similar to the adaptive group
normalization [9] used in U-Net. The second way is referred to as adaptive
layer normalization (AdaLN). Formally,
AdaLN($h,s$)=$y_{s}\text{LayerNorm}(h)+y_{b}$, where $h$ is a hidden state
inside an SSM block, and $y_{s}$ and $y_{b}$ are obtained from a linear
projection of the time embedding. As illustrated in Figure 3 (a), while simple
and straightforward, the first way that treats time as a token performs better
than AdaLN.

####  3.2.4 Scaling model size.

We investigate scaling properties of DiS by studying the effect of depth,
_i.e._ , number of SSM layers, width, e.g. the hidden size. Specifically, we
train 5 DiT models on the ImageNet dataset with a resolution of
256$\times$256, spanning model configurations from small to huge as detailed
in Table 1, denoted as (S, B, M, L, H) for simple. As shown in Figure 3 (b),
the performance improves as the depth increase from 25 to 49. Similarly,
increasing the width from 387 to 768 yields performance gains. Notably, the
model’s performance has not yet reached convergence. Overall, across all five
configurations, significant enhancements in FID are observed across all
training stages by augmenting the depth and width of the SSM structure.

###  3.3 Main Results

Table 4: Benchmarking class-conditional image generation on ImageNet
256$\times$256. DiS-H/2 achieves state-of-the-art FID metrics towards best
competitors. Class-Conditional ImageNet 256$\times$256
---
Model | FID$\downarrow$ | sFID$\downarrow$ | IS$\uparrow$ | Precision$\uparrow$ | Recall$\uparrow$
_GAN_
BigGAN-deep [3] | 6.95 | 7.36 | 171.4 | 0.87 | 0.28
StyleGAN-XL [57] | 2.30 | 4.02 | 265.12 | 0.78 | 0.53
_Diff. based on U-Net_
ADM [9] | 10.94 | 6.02 | 100.98 | 0.69 | 0.63
ADM-U | 7.49 | 5.13 | 127.49 | 0.72 | 0.63
ADM-G | 4.59 | 5.25 | 186.70 | 0.82 | 0.52
ADM-G, ADM-U | 3.94 | 6.14 | 215.84 | 0.83 | 0.53
CDM [27] | 4.88 | - | 158.71 | - | -
LDM-8 [54] | 15.51 | - | 79.03 | 0.65 | 0.63
LDM-8-G | 7.76 | - | 209.52 | 0.84 | 0.35
LDM-4 | 10.56 | - | 103.49 | 0.71 | 0.62
LDM-4-G | 3.60 | - | 247.67 | 0.87 | 0.48
VDM++ [34] | 2.12 | - | 267.70 | - | -
_Diff. based on Transformer_
U-ViT-H/2 [1] | 2.29 | 5.68 | 263.88 | 0.82 | 0.57
DiT-XL/2 [51] | 2.27 | 4.60 | 278.24 | 0.83 | 0.57
_Diff. based on SSM_
DiS-H/2 | 2.10 | 4.55 | 271.32 | 0.82 | 0.58

####  3.3.1 Unconditional image generation.

We conduct a comparative analysis between the DiT-S/2 model and prior
diffusion models based on U-Net as well as Transformer-based diffusion models.
Consistent with previous literature [1], we evaluate image quality using the
FID score, based on 50K generated samples. As summarized in Tables 2 and 3,
DiS exhibits comparable performance to U-Net and U-ViT-S/2 on unconditional
CIFAR10 and CelebA 64$\times$64 datasets, meanwhile significantly
outperforming GenViT. Moreover, notably, our model achieves approximately a
50% reduction in model parameters.

####  3.3.2 Class-conditional image generation.

To demonstrate performance in latent space, where images are initially
converted to their latent representations before applying diffusion models, we
conducted additional experiments. The evaluation results are listed in Tables
4 and 5. On the class-conditional ImageNet 256$\times$256 dataset, our DiS-H/2
achieves an FID of 2.10, surpassing all prior models. Notably, DiS-H/2 also
outperforms DiT-XL/2, a competitive diffusion model employing a transformer as
its backbone. Furthermore, on the class-conditional ImageNet 512$\times$512
dataset, our DiS-H/2 surpasses ADM-G, a model directly modeling image pixels.
Finally, we also observe that DiS-H/2 achieves higher recall values across all
tested classifier-free guidance scales compared to LDM-4 and LDM-8 in various
settings.

###  3.4 Case Study

Figure 4 showcases a selection of samples from ImageNet datasets at
resolutions of 256$\times$256 and 512$\times$512, along with random samples
from other datasets, demonstrating clear semantics and high-quality
generation. For further exploration, additional generated samples, including
both class-conditional and random ones, are available on the project page.

Table 5: Benchmarking class-conditional image generation on ImageNet
512$\times$512. DiT-H/2 demonstrates a promising performance compared with
both CNN-based and Transformer-based UNet for diffusion. Class-Conditional
ImageNet 512$\times$512
---
Model | FID$\downarrow$ | sFID$\downarrow$ | IS$\uparrow$ | Precision$\uparrow$ | Recall$\uparrow$
_GAN_
BigGAN-deep [3] | 8.43 | 8.13 | 177.90 | 0.88 | 0.29
StyleGAN-XL [57] | 2.41 | 4.06 | 267.75 | 0.77 | 0.52
_Diff. based on U-Net_
ADM [9] | 23.24 | 10.19 | 58.06 | 0.73 | 0.60
ADM-U | 9.96 | 5.62 | 121.78 | 0.75 | 0.64
ADM-G | 7.72 | 6.57 | 172.71 | 0.87 | 0.42
ADM-G, ADM-U | 3.85 | 5.86 | 221.72 | 0.84 | 0.53
VDM++ [34] | 2.65 | - | 278.10 | - | -
_Diff. based on Transformer_
U-ViT-H/4 [1] | 4.05 | 6.44 | 263.79 | 0.84 | 0.48
DiT-XL/2 [51] | 3.04 | 5.02 | 240.82 | 0.84 | 0.54
_Diff. based on SSM_
DiS-H/2 | 2.88 | 4.74 | 272.33 | 0.84 | 0.56
Figure 4: Image results generated from DiS model.  Selected samples on
ImageNet 256$\times$256, and random samples on CIFAR10, CelebA 64$\times$64\.
We can see that DiS can generate high-quality images while keeping integrated
condition alignment.

##  4 Related Works

###  4.1 State Space Backbone for Sequence Modeling

State space models, a recent addition to the realm of deep learning, have
garnered attention for their capacity to transform state spaces [20, 19, 31].
Drawing inspiration from continuous state space models in control systems and
leveraging the HiPPO initialization [17], LSSL [20] has demonstrated promise
in addressing long-range dependency issues. However, the computational and
memory demands associated with the state representation render LSSL
impractical for real-world applications. To mitigate this challenge, S4 [19]
proposes normalizing parameters into a diagonal structure. Subsequently,
various iterations of structured state-space models have emerged, featuring
diverse structures such as complex-diagonal configurations [22, 18], multiple-
input multiple-output supporting [59], decomposition of diagonal plus low-rank
operations [23], selection mechanism [16]. These models are then integrated
into large representation models [44, 43, 15]. These models have been
integrated into large representation models [44, 43, 15], primarily focusing
on applications involving long-range and causal data such as language and
speech, including tasks like language understanding [43] and content-based
reasoning [16]. In this paper, we study the scaling properties of SSM when
used as the backbone of diffusion models of images.

###  4.2 Image Generation in Diffusions

Diffusion models [60, 7, 9, 26, 48] represent a category of generative
probabilistic models designed to approximate data distributions and facilitate
straightforward sampling procedures. Typically, these models operate by taking
a Gaussian noise input and iteratively denoising it through a series of
gradual steps until it conforms to the target distribution. The specifics of
the denoising process, including the number of steps and the transformation
parameterization, vary across different studies [60, 61, 42, 41, 39].
Recently, diffusion models have emerged as leading generators renowned for
their ability to learn intricate distributions and produce diverse, high-
quality samples. These models have found successful application across various
domains, including images [9, 48, 14, 55, 53, 54, 13, 12], video [58], 3D
scenes [46], and motion sequences [68, 65]. In this work, we investigate novel
structural enhancements for diffusion implementation. Here we explore a new
structure class for diffusion implementation.

##  5 Conclusion

This paper presents Diffusion State Space models (DiS), a simple and general
state space-based framework for image generation using diffusion models. DiS
adopts a unified approach to handle all inputs, including time, conditions,
and noisy image patches, treating them as concatenated tokens. Experimental
results indicate that DiS compares favorably, if not surpasses, prior CNN-
based or Transformer-based U-Net models while inheriting the remarkable
scalability characteristic of the state space model class. We posit that DiS
can offer valuable insights for future investigations into backbones within
diffusion models and contribute to advancements in generative modeling across
large-scale multimodal datasets. Given the encouraging scalability outcomes
presented in this study, future endeavors should focus on further scaling DiS
to larger models and token counts.

## References

  * [1] Bao, F., Nie, S., Xue, K., Cao, Y., Li, C., Su, H., Zhu, J.: All are worth words: A vit backbone for diffusion models. In: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition. pp. 22669–22679 (2023)
  * [2] Betker, J., Goh, G., Jing, L., Brooks, T., Wang, J., Li, L., Ouyang, L., Zhuang, J., Lee, J., Guo, Y., et al.: Improving image generation with better captions. Computer Science. https://cdn. openai. com/papers/dall-e-3. pdf 2(3), 8 (2023)
  * [3] Brock, A., Donahue, J., Simonyan, K.: Large scale gan training for high fidelity natural image synthesis. arXiv preprint arXiv:1809.11096 (2018)
  * [4] Cao, H., Tan, C., Gao, Z., Xu, Y., Chen, G., Heng, P.A., Li, S.Z.: A survey on generative diffusion models. IEEE Transactions on Knowledge and Data Engineering (2024)
  * [5] Chen, N., Zhang, Y., Zen, H., Weiss, R.J., Norouzi, M., Chan, W.: Wavegrad: Estimating gradients for waveform generation. arXiv preprint arXiv:2009.00713 (2020)
  * [6] Choi, J., Kim, S., Jeong, Y., Gwon, Y., Yoon, S.: Ilvr: Conditioning method for denoising diffusion probabilistic models. arXiv preprint arXiv:2108.02938 (2021)
  * [7] Croitoru, F.A., Hondru, V., Ionescu, R.T., Shah, M.: Diffusion models in vision: A survey. IEEE Transactions on Pattern Analysis and Machine Intelligence (2023)
  * [8] Deng, J., Dong, W., Socher, R., Li, L.J., Li, K., Fei-Fei, L.: Imagenet: A large-scale hierarchical image database. In: 2009 IEEE conference on computer vision and pattern recognition. pp. 248–255. Ieee (2009)
  * [9] Dhariwal, P., Nichol, A.: Diffusion models beat gans on image synthesis. Advances in neural information processing systems 34, 8780–8794 (2021)
  * [10] Dockhorn, T., Vahdat, A., Kreis, K.: Score-based generative modeling with critically-damped langevin diffusion. arXiv preprint arXiv:2112.07068 (2021)
  * [11] Dosovitskiy, A., Beyer, L., Kolesnikov, A., Weissenborn, D., Zhai, X., Unterthiner, T., Dehghani, M., Minderer, M., Heigold, G., Gelly, S., et al.: An image is worth 16x16 words: Transformers for image recognition at scale. arXiv preprint arXiv:2010.11929 (2020)
  * [12] Duan, X., Cui, S., Kang, G., Zhang, B., Fei, Z., Fan, M., Huang, J.: Tuning-free inversion-enhanced control for consistent image editing. arXiv preprint arXiv:2312.14611 (2023)
  * [13] Fei, Z., Fan, M., Huang, J.: Gradient-free textual inversion. arXiv preprint arXiv:2304.05818 (2023)
  * [14] Fei, Z., Fan, M., Huang, J., Wei, X., Wei, X.: Progressive denoising model for fine-grained text-to-image generation. arXiv preprint arXiv:2210.02291 (2022)
  * [15] Fu, D.Y., Dao, T., Saab, K.K., Thomas, A.W., Rudra, A., Ré, C.: Hungry hungry hippos: Towards language modeling with state space models. arXiv preprint arXiv:2212.14052 (2022)
  * [16] Gu, A., Dao, T.: Mamba: Linear-time sequence modeling with selective state spaces. arXiv preprint arXiv:2312.00752 (2023)
  * [17] Gu, A., Dao, T., Ermon, S., Rudra, A., Ré, C.: Hippo: Recurrent memory with optimal polynomial projections. Advances in neural information processing systems 33, 1474–1487 (2020)
  * [18] Gu, A., Goel, K., Gupta, A., Ré, C.: On the parameterization and initialization of diagonal state space models. Advances in Neural Information Processing Systems 35, 35971–35983 (2022)
  * [19] Gu, A., Goel, K., Ré, C.: Efficiently modeling long sequences with structured state spaces. arXiv preprint arXiv:2111.00396 (2021)
  * [20] Gu, A., Johnson, I., Goel, K., Saab, K., Dao, T., Rudra, A., Ré, C.: Combining recurrent, convolutional, and continuous-time models with linear state space layers. Advances in neural information processing systems 34, 572–585 (2021)
  * [21] Gu, S., Chen, D., Bao, J., Wen, F., Zhang, B., Chen, D., Yuan, L., Guo, B.: Vector quantized diffusion model for text-to-image synthesis. In: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition. pp. 10696–10706 (2022)
  * [22] Gupta, A., Gu, A., Berant, J.: Diagonal state spaces are as effective as structured state spaces. Advances in Neural Information Processing Systems 35, 22982–22994 (2022)
  * [23] Hasani, R., Lechner, M., Wang, T.H., Chahine, M., Amini, A., Rus, D.: Liquid structural state-space models. arXiv preprint arXiv:2209.12951 (2022)
  * [24] Heusel, M., Ramsauer, H., Unterthiner, T., Nessler, B., Hochreiter, S.: Gans trained by a two time-scale update rule converge to a local nash equilibrium. Advances in neural information processing systems 30 (2017)
  * [25] Ho, J., Chan, W., Saharia, C., Whang, J., Gao, R., Gritsenko, A., Kingma, D.P., Poole, B., Norouzi, M., Fleet, D.J., et al.: Imagen video: High definition video generation with diffusion models. arXiv preprint arXiv:2210.02303 (2022)
  * [26] Ho, J., Jain, A., Abbeel, P.: Denoising diffusion probabilistic models. Advances in neural information processing systems 33, 6840–6851 (2020)
  * [27] Ho, J., Saharia, C., Chan, W., Fleet, D.J., Norouzi, M., Salimans, T.: Cascaded diffusion models for high fidelity image generation. The Journal of Machine Learning Research 23(1), 2249–2281 (2022)
  * [28] Ho, J., Salimans, T.: Classifier-free diffusion guidance. arXiv preprint arXiv:2207.12598 (2022)
  * [29] Kalman, R.E.: A new approach to linear filtering and prediction problems (1960)
  * [30] Karras, T., Aittala, M., Aila, T., Laine, S.: Elucidating the design space of diffusion-based generative models. Advances in Neural Information Processing Systems 35, 26565–26577 (2022)
  * [31] Kim, D., Shin, S., Song, K., Kang, W., Moon, I.C.: Soft truncation: A universal training technique of score-based diffusion model for high precision score estimation. arXiv preprint arXiv:2106.05527 (2021)
  * [32] Kingma, D., Salimans, T., Poole, B., Ho, J.: Variational diffusion models. Advances in neural information processing systems 34, 21696–21707 (2021)
  * [33] Kingma, D.P., Ba, J.: Adam: A method for stochastic optimization. arXiv preprint arXiv:1412.6980 (2014)
  * [34] Kingma, D.P., Gao, R.: Understanding diffusion objectives as the elbo with simple data augmentation. In: Thirty-seventh Conference on Neural Information Processing Systems (2023)
  * [35] Kingma, D.P., Welling, M.: Auto-encoding variational bayes. arXiv preprint arXiv:1312.6114 (2013)
  * [36] Kong, Z., Ping, W., Huang, J., Zhao, K., Catanzaro, B.: Diffwave: A versatile diffusion model for audio synthesis. arXiv preprint arXiv:2009.09761 (2020)
  * [37] Krizhevsky, A., Hinton, G., et al.: Learning multiple layers of features from tiny images (2009)
  * [38] Kynkäänniemi, T., Karras, T., Laine, S., Lehtinen, J., Aila, T.: Improved precision and recall metric for assessing generative models. Advances in Neural Information Processing Systems 32 (2019)
  * [39] Liu, L., Ren, Y., Lin, Z., Zhao, Z.: Pseudo numerical methods for diffusion models on manifolds. arXiv preprint arXiv:2202.09778 (2022)
  * [40] Liu, Z., Luo, P., Wang, X., Tang, X.: Deep learning face attributes in the wild. In: Proceedings of the IEEE international conference on computer vision. pp. 3730–3738 (2015)
  * [41] Lu, C., Zheng, K., Bao, F., Chen, J., Li, C., Zhu, J.: Maximum likelihood training for score-based diffusion odes by high order denoising score matching. In: International Conference on Machine Learning. pp. 14429–14460. PMLR (2022)
  * [42] Lu, C., Zhou, Y., Bao, F., Chen, J., Li, C., Zhu, J.: Dpm-solver: A fast ode solver for diffusion probabilistic model sampling in around 10 steps. Advances in Neural Information Processing Systems 35, 5775–5787 (2022)
  * [43] Ma, X., Zhou, C., Kong, X., He, J., Gui, L., Neubig, G., May, J., Zettlemoyer, L.: Mega: moving average equipped gated attention. arXiv preprint arXiv:2209.10655 (2022)
  * [44] Mehta, H., Gupta, A., Cutkosky, A., Neyshabur, B.: Long range language modeling via gated state spaces. arXiv preprint arXiv:2206.13947 (2022)
  * [45] Mei, K., Patel, V.: Vidm: Video implicit diffusion models. In: Proceedings of the AAAI Conference on Artificial Intelligence. vol. 37, pp. 9117–9125 (2023)
  * [46] Müller, N., Siddiqui, Y., Porzi, L., Bulo, S.R., Kontschieder, P., Nießner, M.: Diffrf: Rendering-guided 3d radiance field diffusion. In: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition. pp. 4328–4338 (2023)
  * [47] Nash, C., Menick, J., Dieleman, S., Battaglia, P.W.: Generating images with sparse representations. arXiv preprint arXiv:2103.03841 (2021)
  * [48] Nichol, A.Q., Dhariwal, P.: Improved denoising diffusion probabilistic models. In: International Conference on Machine Learning. pp. 8162–8171. PMLR (2021)
  * [49] Orvieto, A., Smith, S.L., Gu, A., Fernando, A., Gulcehre, C., Pascanu, R., De, S.: Resurrecting recurrent neural networks for long sequences. arXiv preprint arXiv:2303.06349 (2023)
  * [50] Parmar, G., Zhang, R., Zhu, J.Y.: On aliased resizing and surprising subtleties in gan evaluation. In: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition. pp. 11410–11420 (2022)
  * [51] Peebles, W., Xie, S.: Scalable diffusion models with transformers. In: Proceedings of the IEEE/CVF International Conference on Computer Vision. pp. 4195–4205 (2023)
  * [52] Poole, B., Jain, A., Barron, J.T., Mildenhall, B.: Dreamfusion: Text-to-3d using 2d diffusion. arXiv preprint arXiv:2209.14988 (2022)
  * [53] Ramesh, A., Dhariwal, P., Nichol, A., Chu, C., Chen, M.: Hierarchical text-conditional image generation with clip latents. arXiv preprint arXiv:2204.06125 1(2), 3 (2022)
  * [54] Rombach, R., Blattmann, A., Lorenz, D., Esser, P., Ommer, B.: High-resolution image synthesis with latent diffusion models. In: Proceedings of the IEEE/CVF conference on computer vision and pattern recognition. pp. 10684–10695 (2022)
  * [55] Saharia, C., Chan, W., Saxena, S., Li, L., Whang, J., Denton, E.L., Ghasemipour, K., Gontijo Lopes, R., Karagol Ayan, B., Salimans, T., et al.: Photorealistic text-to-image diffusion models with deep language understanding. Advances in Neural Information Processing Systems 35, 36479–36494 (2022)
  * [56] Salimans, T., Goodfellow, I., Zaremba, W., Cheung, V., Radford, A., Chen, X.: Improved techniques for training gans. Advances in neural information processing systems 29 (2016)
  * [57] Sauer, A., Schwarz, K., Geiger, A.: Stylegan-xl: Scaling stylegan to large diverse datasets. In: ACM SIGGRAPH 2022 conference proceedings. pp. 1–10 (2022)
  * [58] Singer, U., Polyak, A., Hayes, T., Yin, X., An, J., Zhang, S., Hu, Q., Yang, H., Ashual, O., Gafni, O., et al.: Make-a-video: Text-to-video generation without text-video data. arXiv preprint arXiv:2209.14792 (2022)
  * [59] Smith, J.T., Warrington, A., Linderman, S.W.: Simplified state space layers for sequence modeling. arXiv preprint arXiv:2208.04933 (2022)
  * [60] Sohl-Dickstein, J., Weiss, E., Maheswaranathan, N., Ganguli, S.: Deep unsupervised learning using nonequilibrium thermodynamics. In: International conference on machine learning. pp. 2256–2265. PMLR (2015)
  * [61] Song, J., Meng, C., Ermon, S.: Denoising diffusion implicit models. arXiv preprint arXiv:2010.02502 (2020)
  * [62] Song, Y., Durkan, C., Murray, I., Ermon, S.: Maximum likelihood training of score-based diffusion models. Advances in Neural Information Processing Systems 34, 1415–1428 (2021)
  * [63] Song, Y., Ermon, S.: Generative modeling by estimating gradients of the data distribution. Advances in neural information processing systems 32 (2019)
  * [64] Song, Y., Sohl-Dickstein, J., Kingma, D.P., Kumar, A., Ermon, S., Poole, B.: Score-based generative modeling through stochastic differential equations. arXiv preprint arXiv:2011.13456 (2020)
  * [65] Tevet, G., Raab, S., Gordon, B., Shafir, Y., Cohen-Or, D., Bermano, A.H.: Human motion diffusion model. arXiv preprint arXiv:2209.14916 (2022)
  * [66] Vahdat, A., Kreis, K., Kautz, J.: Score-based generative modeling in latent space. Advances in Neural Information Processing Systems 34, 11287–11302 (2021)
  * [67] Yang, X., Shih, S.M., Fu, Y., Zhao, X., Ji, S.: Your vit is secretly a hybrid discriminative-generative diffusion model. arXiv preprint arXiv:2208.07791 (2022)
  * [68] Yuan, Y., Song, J., Iqbal, U., Vahdat, A., Kautz, J.: Physdiff: Physics-guided human motion diffusion model. In: Proceedings of the IEEE/CVF International Conference on Computer Vision. pp. 16010–16021 (2023)
  * [69] Zhang, L., Rao, A., Agrawala, M.: Adding conditional control to text-to-image diffusion models. In: Proceedings of the IEEE/CVF International Conference on Computer Vision. pp. 3836–3847 (2023)
  * [70] Zhao, M., Bao, F., Li, C., Zhu, J.: Egsde: Unpaired image-to-image translation via energy-guided stochastic differential equations. Advances in Neural Information Processing Systems 35, 3609–3623 (2022)
  * [71] Zhou, Z., Tulsiani, S.: Sparsefusion: Distilling view-conditioned diffusion for 3d reconstruction. In: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition. pp. 12588–12597 (2023)
  * [72] Zhu, L., Liao, B., Zhang, Q., Wang, X., Liu, W., Wang, X.: Vision mamba: Efficient visual representation learning with bidirectional state space model. arXiv preprint arXiv:2401.09417 (2024)


1The Swedish Defence Research Agency, Stockholm, Sweden
(E-mail<EMAIL_ADDRESS><EMAIL_ADDRESS>
# tinySLAM-based exploration with a swarm of nano-UAVs

Mattias Vikgren†,1 and Johan Markdahl1

###### Abstract

This paper concerns SLAM and exploration for a swarm of nano-UAVs. The laser
range finder-based tinySLAM algorithm is used to build maps of the
environment. The maps are synchronized using an iterative closest point
algorithm. The UAVs then explore the map by steering to points selected by a
modified dynamic coverage algorithm, for which we prove a stability result.
Both algorithms inform each other, allowing the UAVs to map out new areas of
the environment and move into them for exploration. Experimental findings
using the nano-UAV Crazyflie 2.1 platform are presented. A key challenge is to
implement all algorithms on the hardware limited experimental platform.

###### keywords:

Swarm robotics, multi-agent systems, SLAM, exploration, nano-UAV, Crazyflie.

##  1 INTRODUCTION

Autonomous swarms of flying robots is an emerging technology with promising
potential. There is an advancing trend towards smaller robots due to
developments in hardware. In particular, nano-robots are well-suited to travel
in narrow and cluttered indoor environments. They are also safe to operate
around humans. Applications that leverage these strengths include tasks like
surveillance [1], search-and-rescue operations [2], and gas leak localization
[3]. To realize such tasks in autonomous swarms, a number of control
algorithms must be employed, including GPS-denied localization, mapping, and
exploration. In this paper we demonstrate such capabilities in experiments
with a swarm of Crazyflie 2.1 which are unmanned aerial vehicles (UAVs) of
size 27 g, 1 dm2.

We couple a low memory and low computational complexity laser range finder
(LRF) based simultaneous localization and mapping (SLAM) algorithm [4] with a
modified swarm exploration algorithm [5] and implement them on-board the
hardware limited Crazyflie 2.1 nano-UAV platform. The two algorithms
synergize; the SLAM maps the environment and the exploration moves the swarm
to previously occluded areas, which the SLAM algorithm can survey. This
behavior is realized in an experiment in two adjacent rooms. These
capabilities have not previously been demonstrated using on-board calculations
in a nano-UAV swarm without an external positioning system [6].

###  1.1 Literature review

The survey paper [6] concerns swarming with micro air vehicles. A key
distinction is made between centralized control _vs._ distributed control.
Likewise there is a difference between position sensing that rely on an
external system such as GPS or a motion caption system _vs._ on-board sensing.
The case of distributed control using only on-board sensing is the most
difficult, with state-of-the-art number of swarm members being in the single
digits [6].

Low complexity distributed control/on-board sensing algorithms used for maze
solving with tiny robots are referred to as bug algorithms [7]. The work [2]
implements a wall following bug algorithm with a swarm of Crazyflie 2.0 UAVs.
The UAVs travel through a complex environment along their preferential
direction and later return to their departure point by following a beacon. The
work [3] uses a smell sensor on a swarm of Crazyflie 2.1 UAVs to detect and
localize gas leaks. Compared to [2, 3] we implement tweaked, light versions of
more complex algorithms on the same Crazyflie hardware.

The area coverage task requires robots to move their sensors to cover an area
of interest. Individual coverage algorithms have a long history (see
_e.g.,_[8] and the survey [9]), and there also exist swarm versions ([10] is
the multi-robot version of [8]). A classic exploration algorithm is the mutual
information based work [11]. That work does not consider flying robots but is
still similar in scope to our setup. An algorithm for persistent exploration
is given in [5]. We prefer a modified version of the algorithm [5] to more
systematic algorithms like [10]. This is due to the algorithm [5] being a
feedback law whereas [10] is open loop. Moreover, in certain scenarios it can
be preferable that the algorithm in [5] is unpredictable, as that makes it
more difficult for mobile adversarial agents to avoid discovery.

Figure 1: The Crazyflie 2.1. It weighs $27$ g and measures $92\textrm{
mm}\times 92\textrm{ mm}\times 29\textrm{ mm}$.

###  1.2 Contribution

The main contribution of this paper is to couple the tiny-SLAM algorithm [4]
with an exploration algorithm [5] such that the two algorithms inform each
other about the environment that a swarm of nano-UAVs operates in. Distributed
algorithms of this complexity that only use on-board calculations and sensing
have not previously been implemented on nano-UAV platforms. There are examples
of complex swarm algorithms being executed on the Crazyflie 2.1 nano-UAV, see
_e.g.,_[12], but they rely on a motion capture system and model predictive
control calculations done on a remote computer. We also modify the exploration
algorithm [5] and prove a stability result for our version.

##  2 Model

All UAVs are quadrotors with planar position and orientation states
$({\mathbf{x}}_{i},\theta_{i})\in{{\mathbb{R}}}^{2}\times(-\pi,\pi]$, where
$i\in\\{1,\ldots,N\\}$ refers to the $i$th UAV. A Crazyflie 2.1 can be
actuated by position and velocity commands. We mainly use velocity commands,
_i.e.,_ we assign $\dot{\mathbf{x}}_{i}={\mathbf{u}}_{i,x}$ and
$\dot{\theta}_{i}=u_{i,\theta}$. The input signals are saturated,
$\|{\mathbf{u}}_{i,x}\|\leq S_{x},\,|u_{i,\theta}|\leq S_{\theta}$, as a
safety measure and to prevent resetting of the Crazyflie platform (in the
experiments, we observe that the extended Kalman filter in a Crazyflie is
reset to zero if the states change quickly).

The camera’s field of view is identified with the forward direction in the
UAV’s body fixed frame, see Fig. 2. The coordinates of the set of points
corresponding to the camera image is a rectangular area centered around a
point
${\mathbf{x}}_{i,c}={\mathbf{x}}_{i}+r_{c}[\cos\theta_{i}\,\sin\theta_{i}]^{T}$,
where $r_{c}$ is a camera parameter. The width of the focus area is $w$ and
the length of the focus area is $l$. The field of view is
${\mathcal{V}}_{i}=\\{{\mathbf{y}}\in{\mathbb{R}}^{2}:\,|\langle{\mathbf{y}}-{\mathbf{x}}_{i,c},{\mathbf{v}}_{i,1}\rangle|\leq
w/2,|\langle{\mathbf{y}}-{\mathbf{x}}_{i,c},{\mathbf{v}}_{i,2}\rangle|\leq
l/2\\},$ where ${\mathbf{v}}_{i,1}=[\cos\theta_{i}\,\sin\theta_{i}]^{T}$,
${\mathbf{v}}_{i,2}=[-\sin\theta_{i}\,\cos\theta_{i}]^{T}$. Note that the
cameras are used for exploration, but not for SLAM.

##  3 Collision avoidance

Each Crazyflie 2.1 is equipped with six LRFs that are used to detect
obstacles, walls, the ceiling, and the floor. Four are directed to the front,
back, left, and right of the UAV. The UAV is programmed to move away from any
obstacle it detects when the distance to the obstacle is lower than a preset
threshold. Because the UAVs are controlled to move with the camera sensor
directed forward, the front LRF detects any obstacle in its path. This type of
naive collision avoidance algorithm may lead to UAVs getting stuck, which can
be avoided by path planning or by only moving towards points that are in the
visible line of sight.

The multi-agent collision avoidance we use in this paper is similar to the
potential-field based method in [5], and indeed in many of other work, see
_e.g.,_[13]. The key difference is that rather than assuming that all agents
can sense the relative localization ${\mathbf{x}}_{i}-{\mathbf{x}}_{j}$, each
agent broadcasts its position ${\mathbf{x}}_{i}$ in a map they have built
using SLAM. This requires the maps to be synchronized, which we achieve by
matching point clouds using an iterative closest point algorithm. The
knowledge of ${\mathbf{x}}_{i}$ and ${\mathbf{x}}_{j}$ allows an agent to
calculate ${\mathbf{x}}_{i}-{\mathbf{x}}_{j}$, which is used for collision
avoidance and swarm algorithms.

##  4 SLAM

The Crazyflie 2.1 platform restricts the choice of SLAM implementation due to
its severe limitations in terms of computational power, sensory information,
and memory. Therefore, the tinySLAM algorithm [4] is chosen for its
lightweight implementation and compatibility with LRF sensors. Due to the
restrictions on the available computational resources, the algorithm is used
without the recommended particle filter as a single hypothesis SLAM method.

Figure 2: Camera model. The vectors ${\mathbf{v}}_{i,1}$ and
${\mathbf{v}}_{i,2}$ are the basis vectors of the $i$th body frame expressed
in the world frame which is aligned with the horizontal and vertical axes.

The tinySLAM map representation has been adopted to the Crazyflie 2.1 system.
The SLAM map ${\mathbf{n}}_{i}$ of UAV $i$ is of occupancy grid type and
discretizes the environment into equally sized square cells, $10\textrm{
cm}\times 10\textrm{ cm}$ in size. Each cell is associated with an 8-bit
integer $({\mathbf{n}}_{i})_{ab}$, _i.e.,_ a number in the range $0$ to $255$,
which represents a degree of confidence in whether a cell is occupied by an
obstacle or not. The number $0$ represents occupied and $255$ free. The cells
are initialized at $127$, representing uncertainty about their state. The map
is a $100\times 100$ matrix, occupying $10^{4}$ out of a total of $5\cdot
10^{4}$ bytes in the Crazyflie 2.1 memory.

The map update is based on Bresenham’s line algorithm but has been modified by
the creator of the tinySLAM algorithm to facilitate faster convergence of the
scan-matching method used for localization [4]. The main modification is what
the authors describe as the ‘hole function’, that blurs measured obstacles to
the neighboring cells in the direction of the scan. This gives the walls a
thickness larger than one cell. The update of the cell values can be tuned
using a confidence parameter $\alpha\in\\{0,\ldots,255\\}$ and a parameter
that sets the hole width. The confidence parameter $\alpha$ defines how the
update of each cell weighs the state of the existing map compared to the new
measurement $y\in\\{0,255\\}$ as

$\displaystyle({\mathbf{n}}_{i})_{ab}(k+1)=((255-\alpha)({\mathbf{n}}_{i})_{ab}(k)+\alpha
y)/255.$ (1)

The tinySLAM algorithm utilizes a Monte Carlo scan matching localization
method. The scan matching compares the latest laser scan to the current map in
different positions decided by a random walk starting in the current estimated
position. The position that the scan matcher finds to best fit to the map is
fed into the extended Kalman filter used for sensor fusion and state
estimation of the UAV.

##  5 Relative localization

Multi-agent systems require relative localization for data fusion and
collision avoidance between the agents. Our method for relative localization
is based on an iterative closest point (ICP) algorithm commonly used for map
merging or scan matching in multi-agent SLAM systems [14]. The ICP algorithm
is an iterative point registry method which tries to find the rigid
transformation that minimizes the Euclidean distance between two point sets
with unknown correspondence. The two point sets correspond to obstacle
coordinates in two maps ${\mathbf{n}}_{i}$ and ${\mathbf{n}}_{j}$ created by
two different UAVs. The ICP algorithm provides an estimate of the
transformation between the two agent’s local maps. It is used to transform the
locally estimated position being communicated between the agents for relative
localization.

Due to the limited communication capabilities and memory available on the
Crazyflie 2.1 we suggest an approach for selecting a few dozen points
representing the obstacles in map ${\mathbf{n}}_{i}$ created by the SLAM
algorithm, rather than using all points marked as obstacles in the map. These
points are extracted by iterating over the map to find obstacle cells in the
map. The coordinates of the obstacle cell is added to the point set if there
are no points already added to the set within a certain distance threshold.
The threshold has to be tuned to find equally spaced points covering the
entire area that has been mapped by the SLAM algorithm without consuming too
much memory.

The transform is found through a matching and transformation stage, which is
repeated until convergence. The matching stage is started with finding the
correspondence between the two point sets. Each point ${\mathbf{p}}_{i,k}$ in
the set ${\mathcal{S}}_{i}$ of UAV $i$ is associated with the point
${\mathbf{q}}_{j,l}$ in the set ${\mathcal{S}}_{j}$ of UAV $j$ having the
minimum Euclidean distance $\|{\mathbf{p}}_{i,k}-{\mathbf{q}}_{j,l}\|$. If the
distance is larger than a certain tolerance, then the point is ignored and it
is seen as an outlier in this iteration.

The transformation stage is based on a least square solver to find the optimal
transformation between the two point sets. First calculate the centroids
$\widebar{{\mathbf{p}}}_{i}$ and $\widebar{{\mathbf{q}}}_{j}$

$\displaystyle\widebar{{\mathbf{p}}}_{i}=\frac{1}{K}\sum_{k=1}^{K}{\mathbf{p}}_{i,k},\quad\widebar{{\mathbf{q}}}_{j}=\frac{1}{K}\sum_{l=1}^{K}{\mathbf{q}}_{j,l},$

and the distance of each point to the centroid, ${\mathbf{p}}_{i,k}^{\prime}$
and ${\mathbf{q}}_{j,l}^{\prime}$,

$\displaystyle{\mathbf{p}}_{i,k}^{\prime}={\mathbf{p}}_{i,k}-\widebar{{\mathbf{p}}}_{i},\quad{\mathbf{q}}_{j,l}^{\prime}={\mathbf{q}}_{j,l}-\widebar{{\mathbf{q}}}_{j}.$

Use the distances to the centroids to form matrices

$\displaystyle{\mathbf{P}}_{i}^{\prime}=\begin{bmatrix}{\mathbf{p}}_{i,1}^{\prime}\,{\mathbf{p}}_{i,2}^{\prime}\,\ldots{\mathbf{p}}_{i,K}^{\prime}\end{bmatrix},\quad{\mathbf{Q}}_{j}^{\prime}=\begin{bmatrix}{\mathbf{q}}_{j,1}^{\prime}\,{\mathbf{q}}_{j,2}^{\prime}\,\ldots{\mathbf{q}}_{j,K}^{\prime}\end{bmatrix}.$

Then calculate the cross-variance matrix ${\mathbf{H}}_{ij}$,

$\displaystyle{\mathbf{H}}_{ij}={\mathbf{\textasciiacute
P}}_{i}^{\prime}{\mathbf{Q}}_{j}^{\prime\mathsf{T}}\in{\mathbb{R}}^{2\times
2}.$

The matrix ${\mathbf{H}}_{ij}$ yields the rotation angle $\vartheta_{ij}$ as

$\displaystyle\vartheta_{ij}=\mathrm{atan}2(({\mathbf{H}}_{ij})_{12}-({\mathbf{H}}_{ij})_{21},({\mathbf{H}}_{ij})_{11}+({\mathbf{H}}_{ij})_{22}).$

Now the points can be transformed using the rotation matrix
${\mathbf{R}}(\vartheta_{ij})$ and the translational shift
${\mathbf{t}}_{ij}$, where

$\displaystyle{\mathbf{R}}(\vartheta_{ij})=\begin{bmatrix}\cos(\vartheta_{ij})&-\sin(\vartheta_{ij})\\\
\sin(\vartheta_{ij})&\cos(\vartheta_{ij})\end{bmatrix},{\mathbf{t}}_{ij}=\widebar{{\mathbf{q}}}_{j}-{\mathbf{R}}(\vartheta_{ij})\widebar{{\mathbf{p}}}_{i}.$

The rotation matrix ${\mathbf{R}}(\vartheta_{ij})$ and translation
${\mathbf{t}}_{ij}$ is applied to points in set ${\mathcal{S}}_{j}$ and the
next iteration starts from there.

Each UAV in the swarm has a predetermined neighbor with which it is supposed
to synchronize its map. The UAVs that synchronize their maps are chosen such
that a chain of transformations are created and communicated. This enables all
UAVs to transform coordinates from all other UAVs to its own map coordinates.
Note that this allows for a distributed solution, without all-to-all
communication. Moreover, the pairs of communicating neighbors could be
determined dynamically based on current relative distances rather be
predetermined.

##  6 Exploration control design

Part of the exploration control design is taken from the dynamic coverage
algorithm [5]. However, we modify it for use on a low memory, low
computational resources hardware. Moreover, we add an orientation state and
design the control law to couple the position and orientation errors.

###  6.1 The camera coverage map

Each robot has two maps of its own. The map ${\mathbf{n}}_{i}$ created by the
SLAM algorithm keeps track of obstacles and walls in the environment. A second
map ${\mathbf{m}}_{i}$ is used to keep track of when points in the environment
were last seen by the on-board cameras. The camera coverage map
${\mathbf{m}}_{i}$ is an $A\times B$ matrix where each element
$({\mathbf{m}}_{i})_{ab}$ is an 8-bit integer, _i.e.,_
${\mathbf{m}}_{i}\in\\{0,1,\ldots,255\\}^{A\times B}$.

A discrete time-step $k\in\\{1,\ldots,255\\}$ is used to time-stamp the map.
Each time step has a length $\Delta t$. Let $t$ denote the current time. For
$t\in[\Delta tk,\Delta t(k+1))$, the map is updated as $({\mathbf{m}}_{i})=k$
for all ${\mathbf{y}}$ in the field of view ${\mathcal{V}}_{i}$. The map
${\mathbf{m}}_{i}$ can be interpreted as a landscape of time-stamps that keeps
track of the age of information. This limits the total time of exploration to
$255\Delta t$, but with $\Delta t\approx 1$ that is still larger than the
battery allows.

Our approach to exploration is based on [5]. However, their algorithm is not
designed for nano-UAVs. In particular, they require any pair of agents $(i,j)$
that are sufficiently close to one-another to exchange their maps
${\mathbf{m}}_{i}$ and ${\mathbf{m}}_{j}$. The Crazyflie 2.1 communication
bandwidth is too limited for this. Instead, we assume that agents exchange
only their current position and orientation in a local map. This suffices for
collision avoidance and agent $i$ can timestamp its map ${\mathbf{m}}_{i}$
based on receiving $({\mathbf{x}}_{j},\theta_{j})$ for all $j$ (or just a
subset of neighbors for a distributed solution).

###  6.2 Position control design

The exploration algorithm strives to make the UAV see its desired camera point
${\mathbf{x}}_{i,c,d}\in{{\mathbb{R}}}^{2}$ inside its field of view
${\mathcal{V}}_{i}$. This is done by solving an optimization problem using the
brute force method. The calculation requires ${\mathcal{O}}(AB)$ operations,
where $A\times B$ is the dimension of the map ${\mathbf{m}}_{i}$. This is the
most expensive calculation of our algorithm, but it is only done infrequently
when updating ${\mathbf{x}}_{i,c,d}$ as it enters the field of view
${\mathcal{V}}_{i}$.

As detailed in [5], the next desired camera point ${\mathbf{x}}_{i,c,d}$ is
calculated as

$\displaystyle{\mathbf{x}}_{i,c,d}$
$\displaystyle=\operatorname*{arg\,min}_{{\mathbf{y}}}(t-({\mathbf{m}}_{i})_{ab})(\lambda+(1-\lambda)f_{i}({\mathbf{y}})),$
(2) $\displaystyle f_{i}({\mathbf{y}})$
$\displaystyle=\exp(-\sigma_{1}\|{\mathbf{y}}-{\mathbf{x}}_{i,c}\|-\sigma_{2}\|{\mathbf{y}}-{\mathbf{x}}_{i,c,d}^{\prime}\|).$

The function $f_{i}$ is used to penalize points ${\mathbf{y}}$ for being
either too far from the UAV’s current camera position ${\mathbf{x}}_{i,c}$ or
too far from the previous desired camera point
${\mathbf{x}}_{i,c,d}^{\prime}$. The parameter $\lambda\in[0,1]$ determines if
the UAVs travel as a group or spread out. In particular, for $\lambda=1$, all
the UAVs will calculate the same next desired camera point. The current
desired UAV orientation $\theta_{i,d}$ is calculated based on
${\mathbf{x}}_{i,c}$ and ${\mathbf{x}}_{i,c,d}$.

The camera position error, heading error, exploration control law, and closed
loop system are defined by

$\displaystyle{\mathbf{e}}_{i,c}:\\!$
$\displaystyle={\mathbf{x}}_{i,c,d}-{\mathbf{x}}_{i,c}={\mathbf{x}}_{i,c,d}-{\mathbf{x}}_{i}-r_{c}{\mathbf{v}}_{i,1},$
$\displaystyle e_{i,\theta}:\\!$
$\displaystyle=\theta_{i,d}-\theta_{i}\in(-\pi,\pi],$
$\displaystyle{\mathbf{u}}_{i}:\\!$
$\displaystyle=k_{c}\,{\mathrm{sech}}^{2}(k_{s}e_{i,\theta})\tanh(k_{t}\|{\mathbf{e}}_{i,c}\|)\,{\mathbf{e}}_{i,c}/\|{\mathbf{e}}_{i,c}\|,$
(3) $\displaystyle\dot{\mathbf{x}}_{i}$ $\displaystyle={\mathbf{u}}_{i},$
$\displaystyle\dot{\mathbf{e}}_{i,c}$
$\displaystyle=-{\mathbf{u}}_{i}-\dot{\theta}_{i}r_{c}{\mathbf{v}}_{i,2},$

where $k_{t},k_{s},k_{c}\in(0,\infty)$. The control ${\mathbf{u}}_{i,c}$ is
designed to steer the center of its camera field of view,
${\mathbf{x}}_{i,c}$, to the desired point ${\mathbf{x}}_{i,c,d}$. The unit
vector and $\tanh$ function make $k_{c}$ easier to tune. The UAVs should
rotate before they translate. This allows the camera and LRF to align with the
direction of motion. Note that this also excludes the use of a feedforward
term in ${\mathbf{u}}_{i}$ to cancel the ${\mathbf{v}}_{i,2}$ term in
$\dot{\mathbf{e}}_{i,c}$. To achieve fast convergence to $\theta_{d}$, we
could use a high-gain control law. However, we observe that high angular
velocities lead to resetting of the Kalman filter of the Crazyflie 2.1, which
interferes with the SLAM algorithm. We hence choose a control design that
couples the position and heading errors, where the $\mathrm{sech}$ function is
shaped like a Gaussian and makes translational velocities small when
$k_{s}|e_{i,\theta}|\gg 0$.

The control (3) differs from that in [5] in that it couples the position and
orientation error (note that there is no orientation state in [5]). To gain
intuition for our approach, we could consider a control action that consists
of two separate maneuvers. First, a rotation to align the UAV orientation
${\mathbf{v}}_{i,1}$ with the line that connects ${\mathbf{x}}_{i,c}$ and
${\mathbf{x}}_{i,c,d}$. Second, a translation that takes ${\mathbf{x}}_{i,c}$
to ${\mathbf{x}}_{i,c,d}$. However, that sequence of maneuvers would be an
open loop control which would perform poorly due to errors and disturbances.
Our feedback control (3) fuses these two separate maneuvers into a single
smooth motion. The notion of two separate motions is still helpful; it is also
the intuition for the proof of our main theoretical result, Theorem 1 in the
next section.

###  6.3 Heading control design

The heading $\theta_{i}\in(-\pi,\pi]$ or yaw angle of the $i$th UAV is
controlled to align with the desired direction of travel
$\theta_{i,d}={\mathrm{atan}}2(({\mathbf{e}}_{i,c})_{2},({\mathbf{e}}_{i,c})_{1})$.
The Crazyflie 2.1 platform allows control by angular velocity commands,
$\dot{\theta}_{i}=u_{i,\theta}$, where $u_{i,\theta}\in{{\mathbb{R}}}$ is the
input. The heading error is

$\displaystyle e_{i,\theta}=\theta_{i,d}-\theta_{i}\in(-\pi,\pi].$

The heading control and closed-loop system is given by

$\displaystyle u_{i,\theta}:=k_{\theta}\frac{\sin e_{i,\theta}}{\sqrt{1+\cos
e_{i,\theta}}},\quad\dot{\theta}_{i}=u_{i,\theta}.$ (4)

The trigonometric functions $\sin$ and $\cos$ remove the discontinuity due to
$e_{i,\theta}\in(-\pi,\pi]$. The control (4) is the ${\mathsf{SO}}(2)$ version
of the large angle control law on ${\mathsf{SO}}(3)$ proposed in [15]. It
achieves its highest gain at $e_{i,\theta}\approx\pi$, leading to fast
convergence in the case of large errors. This is desirable since our position
control is slow for large rotation errors $e_{i,\theta}$. Note that the
control (4) is discontinuous at $e_{i,\theta}=\pi$.

##  7 Main results

Figure 3: Local views of the environment given by the SLAM maps
${\mathbf{n}}_{i}$ (heatmaps) and the three UAV position and orientation
states (triangles and colored paths). The local perspectives are those of the
blue UAV (left), the red UAV (middle) and the gray UAV (right) as can be
inferred from the yellow areas of obstacle free grid points surrounding their
trajectories. The initial positions of all UAVs are given by black circles and
the final position by black squares.

The results of this paper are: a theoretical proof of stability of the
modified exploration algorithm based on [5] and an experiment with a nano-UAV
swarm with the Crazyflie 2.1 platform. Note that we have run many more
experiments than those displayed here. We strive to provide results that guide
intuition rather than quantitative characterization.

###  7.1 Stability analysis

A number of theoretical results about the control algorithm on which our
approach is based are given in [5]. We extend the algorithm of [5] to include
control of the UAV orientation. The dynamics of the error terms are given by

$\displaystyle\dot{\mathbf{e}}_{i,c}={}$ $\displaystyle-
k_{c}\,{\mathrm{sech}}^{2}(k_{s}e_{i,\theta})\tanh(k_{t}\|{\mathbf{e}}_{i,c}\|)\,{\mathbf{e}}_{i,c}/\|{\mathbf{e}}_{i,c}\|+$
$\displaystyle-k_{\theta}r_{c}\frac{\sin e_{i,\theta}}{\sqrt{1+\cos
e_{i,\theta}}}\begin{bmatrix}-\sin(\theta_{i,d}-e_{i,\theta})\\\
\cos(\theta_{i,d}-e_{i,\theta})\end{bmatrix},$ (5)
$\displaystyle\dot{e}_{i,\theta}={}$ $\displaystyle-k_{\theta}\sin
e_{i,\theta}/\sqrt{1+\cos e_{i,\theta}}.$ (6)

The system has a triangular structure where $\dot{\mathbf{e}}_{i,c}$ depends
on ${\mathbf{e}}_{i,c}$ and $e_{i,\theta}$, whereas $\dot{e}_{i,\theta}$
depends on $e_{i,\theta}$ but not on ${\mathbf{e}}_{i,c}$. We leverage this
structure to prove a stability result for our version of the algorithm [5].
First we need some theory.

######  Lemma 1 (Vidyasagar [16]).

The origin $({\mathbf{0}},{\mathbf{0}})$ of a system

$\displaystyle\dot{\mathbf{x}}={\mathbf{f}}({\mathbf{x}},{\mathbf{y}}),\quad\dot{\mathbf{y}}={\mathbf{g}}({\mathbf{y}}),$
(7)

is asymptotically stable if and only if ${\mathbf{0}}$ and ${\mathbf{0}}$ are
asymptotically stable equilibria of
$\dot{\mathbf{x}}={\mathbf{f}}({\mathbf{x}},{\mathbf{0}})$ and
$\dot{\mathbf{y}}={\mathbf{g}}({\mathbf{y}})$.

The notion of almost global stability of an equilibrium corresponds to the
idea of stability and convergence from all initial conditions except for a set
that is negligible in size.

######  Definition 1.

An equilibrium ${\mathbf{y}}\in\mathcal{X}$ of a dynamical system
$\dot{\mathbf{x}}={\mathbf{f}}({\mathbf{x}})$ on a space $\mathcal{X}$ is
_almost globally asymptotically stable_ if it is stable and
$\lim_{t\rightarrow\infty}{\mathbf{x}}(t;{\mathbf{x}}_{0})={\mathbf{y}}$ for
all ${\mathbf{x}}_{0}\in\mathcal{X}\backslash\mathcal{N}$, where $\mathcal{N}$
is a set of Lebesgue measure zero.

What follows is our main theoretical result:

######  Theorem 1.

The origin $({\mathbf{0}},0)$ is an almost globally asymptotically stable
equilibrium of the error dynamics (5)–(6).

###### Proof.

The error dynamics (5)–(6) are on the triangular form (7). Moreover, the two
decoupled subsystems,

$\displaystyle\dot{\mathbf{e}}_{i,c}|_{e_{\theta}=0}$
$\displaystyle=-k_{c}\tanh(k_{t}\|{\mathbf{e}}_{i,c}\|)\,{\mathbf{e}}_{i,c}/\|{\mathbf{e}}_{i,c}\|,$
$\displaystyle\dot{e}_{i,\theta}$ $\displaystyle=-k_{\theta}\sin
e_{i,\theta}/\sqrt{1+\cos e_{i,\theta}},$

both have an asymptotically stable equilibrium in ${\mathbf{0}}$ and $0$
respectively. By Lemma 1, $({\mathbf{0}},0)$ is an asymptotically stable
equilibrium of (5)–(6). In particular, $({\mathbf{0}},0)$ is stable.

It remains to prove almost global attractiveness of $({\mathbf{0}},0)$. We
will use the comparison principle [17] for a Lyapunov function based on
${\mathbf{e}}_{i,c}$. It is straightforward to show that $e_{i,\theta}$
converges to $0$ for all initial conditions except $\mathcal{N}:=\\{\pi\\}$,
_i.e.,_ from almost all initial conditions. Note that Equation (6) can be
solved for the state $e_{i,\theta}$ as a function of time. This allows us to
substitute $e_{i,\theta}(t)$ into (5) and consider the evolution of
${\mathbf{e}}_{i,c}$ as a non-autonomous system on its own.

Introduce a candidate Lyapunov function
$U:=\frac{1}{2}\|{\mathbf{e}}_{i,c}\|^{2}>0$ which satisfies

$\displaystyle\dot{U}=$
$\displaystyle{}-k_{c}\,{\mathrm{sech}}^{2}(k_{s}e_{i,\theta})\tanh(k_{t}\|{\mathbf{e}}_{i,c}\|)\|{\mathbf{e}}_{i,c}\|+$
$\displaystyle-\frac{k_{\theta}r_{c}}{2}\frac{\sin e_{i,\theta}}{\sqrt{1+\cos
e_{i,\theta}}}\left\langle{\mathbf{e}}_{i,c},\begin{bmatrix}-\sin(\theta_{i,d}-e_{i,\theta})\\\
\cos(\theta_{i,d}-e_{i,\theta})\end{bmatrix}\right\rangle.$

Let $V=\sqrt{U}=\|{\mathbf{e}}_{i,c}\|/\sqrt{2}$, yielding
$\dot{V}=\dot{U}/(2V)$, _i.e.,_

$\displaystyle\dot{V}=$
$\displaystyle{}-\frac{k_{c}}{2}\,{\mathrm{sech}}^{2}(k_{s}e_{i,\theta})\tanh(k_{t}V)$
$\displaystyle-\frac{k_{\theta}r_{c}}{4}\frac{\sin e_{i,\theta}}{\sqrt{1+\cos
e_{i,\theta}}}\left\langle\frac{{\mathbf{e}}_{i,c}}{V},\begin{bmatrix}-\sin(\theta_{i,d}-e_{i,\theta})\\\
\cos(\theta_{i,d}-e_{i,\theta})\end{bmatrix}\right\rangle$ $\displaystyle\leq$
$\displaystyle{}-\frac{k_{c}}{2}\,{\mathrm{sech}}^{2}(k_{s}e_{i,\theta})\tanh(k_{t}V)+\frac{k_{\theta}r_{c}}{2\sqrt{2}}\frac{\sin|e_{i,\theta}|}{\sqrt{1+\cos
e_{i,\theta}}},$

where we used the Cauchy-Schwarz inequality.

By definition of $e_{i,\theta}\rightarrow 0$ as $t\rightarrow\infty$, for any
$\delta>0$ there is a time $T$ such that $|e_{i,\theta}(t)|\leq\delta$ for all
$t\geq T$. In particular, for any small $\delta\in(0,1)$ there is a
$T(\delta)$ such that

$\displaystyle\dot{V}\leq-\frac{k_{c}}{2}\,\tanh(k_{t}V)+\frac{k_{\theta}r_{c}}{2\sqrt{2}}\delta,$

for all $t\geq T(\delta)$. Introduce a variable $W$ given by

$\displaystyle\dot{W}:=-\frac{k_{c}}{2}\,\tanh(k_{t}W)+\frac{k_{\theta}r_{c}}{2\sqrt{2}}\delta,$

Let $V(0)\leq W(0;\delta)$, then $V(t)\leq W(t;\delta)$ by the comparison
principle [17]. For small $\delta$, it can be shown that

$\displaystyle\lim_{t\rightarrow\infty}W(t;\delta)=\frac{1}{k_{t}}{\mathrm{atanh}}\,\frac{k_{\theta}r_{c}\delta}{\sqrt{2}k_{c}}.$

Note that $\lim_{\delta\rightarrow 0}W(t;\delta)=0$. It follows that for
sufficiently large $t$, $V(t)$ is bounded above by an arbitrarily small number
and bounded below by $0$. Hence $\lim_{t\rightarrow\infty}V(t)=0$.∎

######  Remark 1.

As desired, $({\mathbf{x}}_{i,c,d},\theta_{d})$ is an almost globally
asymptotically stable equilibrium of the UAV system

$\displaystyle\dot{\mathbf{x}}_{i,c}={\mathbf{u}}_{i,c},\quad\dot{\theta}_{i}=u_{i,\theta},$

in the case of a constant ${\mathbf{x}}_{i,c,d}$. This means that the dynamics
(5)–(6) and our preceding analysis concerns only one maneuver from
$({\mathbf{x}}_{i,c},\theta_{i})$ to $({\mathbf{x}}_{i,c,d},\theta_{i,d})$.
Note that in practice, we have non-constant
$({\mathbf{x}}_{i,c,d},\theta_{i,d})$ and require ${\mathbf{x}}_{i,c}$ to
converge to a set, the camera field of view ${\mathcal{V}}_{i}$, rather than a
point ${\mathbf{x}}_{i,c,d}$. The field of view ${\mathcal{V}}_{i}$ will
contain the desired equilibrium ${\mathbf{x}}_{i,c,d}$ in finite time. Then,
as ${\mathbf{x}}_{i,c,d}$ is updated, ${\mathcal{V}}_{i}$ will move towards a
new set containing the next equilibrium.

###  7.2 Nano-UAV swarm experiment

The combined SLAM, ICP, and exploration algorithms are implemented on a swarm
of three Crazyflie 2.1 nano-UAVs. The environment that the swarm operates
consists of two adjacent rooms. The UAVs start out close to each other,
aligned with a nearby wall in the environment. The nearby wall provides a
source of LRF sensor readings that the SLAM and ICP algorithms can exploit for
initial synchronization of the three SLAM maps. The UAVs are stabilized to
different heights as to not interfere with each other’s LRF.

The SLAM map and trajectories taken by the UAVs are displayed in Fig. 3. The
UAVs start out close to each other on a line and are by default located in the
center of their maps. After the initial SLAM and ICP phases, the UAVs start to
explore the environment. The local view that each UAV has of its environment
is displayed in the map. The UAV detects the walls of the environment and
explores both its inside and boundaries. Two of the UAVs stay in the first
room. The third UAV discovers the adjacent second room and moves into it for
exploration.

The camera coverage map ${\mathbf{m}}_{i}$ with the timestamped field of view
${\mathcal{V}}_{i}$ of the three UAVs is displayed in Fig. 4. The three
timestamp maps are different since the belief that each UAV has about the
positions of the other members of the swarm differs somewhat, see Fig. 3. Fig.
4 also displays the total coverage obtained as the percentage of all strictly
positive map elements $({\mathbf{m}}_{i})_{ab}>0$. The increase in coverage is
approximately linear. Note that the camera coverage maps does not account for
walls in the environment.

Figure 4: The camera coverage maps ${\mathbf{m}}_{i}$ for the three UAVs. The
total coverage of ${\mathbf{m}}_{i}$ as a function of time (bottom right). The
coverage is similar for each UAV, the colors of the three curves are the same
as those in Fig. 3.

##  8 Conclusions

This paper is a proof of concept study that demonstrates SLAM and exploration
in a swarm of three nano-UAVs. In particular, SLAM has not been implemented
onboard the Crazyflie 2.1 platform previously, let alone in a swarm. We couple
the SLAM algorithm with an exploration algorithm that has been modified for
implementation on the limited hardware of the experimental platform. The two
algorithms synergize; the SLAM algorithm maps out new areas and the
exploration algorithm moves the UAV swarm into these areas from where they can
map out further areas. The agents interact and solve the exploration task in a
cooperative fashion. We also prove an almost global stability result for the
our modified UAV position and orientation control algorithm.

## References

  * [1] M. Saska, V. Vonásek, J. Chudoba, J. Thomas, G. Loianno, and V. Kumar, “Swarm distribution and deployment for cooperative surveillance by micro-aerial vehicles,” Journal of Intelligent & Robotic Systems, vol. 84, no. 1, pp. 469–492, 2016.
  * [2] K. McGuire, C. De Wagter, K. Tuyls, H. Kappen, and G. de Croon, “Minimal navigation solution for a swarm of tiny flying robots to explore an unknown environment,” Science Robotics, vol. 4, no. 35, 2019.
  * [3] B. Duisterhof, S. Li, J. Burgués, and G. Reddi, V.J.and de Croon, “Sniffy bug: A fully autonomous swarm of gas-seeking nano quadcopters in cluttered environments,” in Proceedings of the 34th IEEE/RSJ International Conference on Intelligent Robots and Systems, pp. 9099–9106, 2021.
  * [4] B. Steux and O. Hamzaoui, “tinySLAM: A SLAM algorithm in less than 200 lines C-language program,” in Proceedings of the 11th International Conference on Control Automation Robotics Vision, pp. 1975–1979, 2010.
  * [5] N. Ganganath, C.-T. Cheng, and K. Chi, “Distributed antiflocking algorithms for dynamic coverage of mobile sensor networks,” IEEE Transactions on Industrial Informatics, vol. 12, no. 5, pp. 1795–1805, 2016.
  * [6] M. Coppola, K. McGuire, C. De Wagter, and G. De Croon, “A survey on swarming with micro air vehicles: Fundamental challenges and constraints,” Frontiers in Robotics and AI, vol. 7, p. 18, 2020.
  * [7] K. McGuire, G. de Croon, and K. Tuyls, “A comparative study of bug algorithms for robot navigation,” Robotics and Autonomous Systems, vol. 121, p. 103261, 2019.
  * [8] Y. Gabriely and E. Rimon, “Spanning-tree based coverage of continuous areas by a mobile robot,” Annals of Mathematics and Artificial Intelligence, vol. 31, no. 1, pp. 77–98, 2001.
  * [9] E. Galceran and M. Carreras, “A survey on coverage path planning for robotics,” Robotics and Autonomous systems, vol. 61, no. 12, pp. 1258–1276, 2013.
  * [10] N. Hazon and G. Kaminka, “Redundancy, efficiency and robustness in multi-robot coverage,” in Proceedings of the 22th IEEE International Conference on Robotics and Automation, pp. 735–741, 2005.
  * [11] W. Burgard, M. Moors, C. Stachniss, and F. E. Schneider, “Coordinated multi-robot exploration,” IEEE Transactions on Robotics, vol. 21, no. 3, pp. 376–386, 2005.
  * [12] E. Soria, F. Schiano, and D. Floreano, “Predictive control of aerial swarms in cluttered environments,” Nature Machine Intelligence, vol. 3, no. 6, pp. 545–554, 2021.
  * [13] R. Olfati-Saber, “Flocking for multi-agent dynamic systems: Algorithms and theory,” IEEE Transactions on Automatic Control, vol. 51, no. 3, pp. 401–420, 2006.
  * [14] H. Sobreira, C. Costa, I. Sousa, L. Rocha, J. Lima, P. Farias, P. Costa, and A. Moreira, “Map-matching algorithms for robot self-localization: A comparison between perfect match, iterative closest point and normal distributions transform,” Journal of Intelligent & Robotic Systems, vol. 93, 03 2019.
  * [15] T. Lee, “Exponential stability of an attitude tracking control system on so (3) for large-angle rotational maneuvers,” Systems & Control Letters, vol. 61, no. 1, pp. 231–237, 2012.
  * [16] M. Vidyasagar, “Decomposition techniques for large-scale systems with nonadditive interactions: Stability and stabilizability,” IEEE Transactions on Automatic Control, vol. 25, no. 4, pp. 773–779, 1980.
  * [17] H. Khalil, Nonlinear Systems.  Prentice Hall, 2002.