# Multi-agent Systems for Quadcopters

<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract.

Unmanned Aerial Vehicles (UAVs) have been increasingly used in the context of
remote sensing missions such as target search and tracking, mapping, or
surveillance monitoring. In the first part of our paper we consider agent
dynamics, network topologies, and collective behaviors. The objective is to
enable multiple UAVs to collaborate toward a common goal, as one would find in
a remote sensing setting. An agreement protocol is carried out by the multi-
agents using local information, and without external user input. The second
part of the paper focuses on the equations of motion for a specific type of
UAV, the quadcopter, and expresses them as an affine nonlinear control system.
Finally, we illustrate our work with a simulation of an agreement protocol for
dynamically sound quadcopters augmenting the particle graph theoretic approach
with orientation and a proper dynamics for quadcopters.

###### Key words and phrases:

Quadcopter, rigid body dynamics, graph Laplacian, multi-agent coordination,
network consensus, weighted network, decentralized control

###### 1991 Mathematics Subject Classification:

Primary: 70E60; 93A14; Secondary: 70E55.

∗ Corresponding author: M. Chyba

Richard Carney, Monique Chyba∗, Chris Gray, and George Wilkens

Department of Mathematics, University of Hawaii at Manoa

2565 McCarthy Mall

Honolulu, Hawaii 96822, USA

Corey Shanbrom

Department of Mathematics and Statistics

California State University, Sacramento

6000 J St., Sacramento, CA 95819, USA

##  1\. Introduction

Unmanned Aerial Vehicles (UAVs) present unique opportunities to investigate
unsafe or otherwise unreachable locations. UAVs are already being used in
numerous agricultural applications ([20, 17]), search and rescue ([18, 33]),
and conservation biology ([34]). While there are still many open issues in the
use of UAVs, they also have much to offer when it comes to disaster
preparation and response ([19]).

In Spring 2018, motivated by the 2018 Kilauea volcano eruption, we conducted
work in collaboration with NASA Jet Propulsion Laboratory (JPL) to develop a
real-time data relay system for UAVs. A relay system augmented with an
interactive visual aid for data viewing was developed and implemented to
exploit the live stream, including image superposition between sequential
flights to rapidly provide information about changes to the situation under
study. In September 2018, we demonstrated the concept at the University of
Hawai’i at Mānoa with live streaming of images from Bellows Air Force Station
and video from Fissure 8 on Hawai’i Island. We then developed on-board
capability for image analysis for real time situation assessment ([12]). In
parallel, we also started to study motion planning for a swarm of UAVs ([15]),
which is the problem we focus on in this paper. See [2] for a video
demonstrating some of the work.

Implementing a coordinated group of UAVs, rather than a single vehicle, is
naturally more efficient for many mission types. These include tracking,
surveillance, and mapping ([14]). In this regard, there are numerous reasons,
driven by the variety of applications, why collective control of a group of
UAVs is important. Common objectives include synchronization, formation,
flocking, consensus, and rendezvous. Rendezvous is the primary focus here
([21]). Moreover, autonomy of individual UAVs provides a flexibility and
adaptability that can be beneficial in many kinds of applications.
Optimization techniques, as an approach to collective control, have the
general drawback of scaling poorly with the number of vehicles. We shall focus
on a common alternative: the so-called graph theoretic methods. The relative
simplicity of these methods makes them an attractive choice in terms of both
scalability and the desire for real-time control ([25, 27]). Moreover, they
apply quite naturally to situations where there is an existing network of
communications between the UAVs, although the existence of such is also an
additional assumption being made. The benefits include control that is
decentralized in that it occurs over the network without any advanced
planning, and certain objectives, such as rendezvous, can be carried out with
only relative information being shared across the network.

In this paper we focus on path-planning of the most common type of UAV,
quadcopters, and we focus on the dynamics of this specific type of UAV in
Sections 3 and 4. Section 2, however, applies to any multi-agent UAV path
configuration since we treat them as particle points. We also occasionally
refer to any UAV by the common, though imprecise, more general term “drone”.
These small versatile vehicles have proved to be a major asset in the
increasing demand for complex autonomous robotics missions, particularly those
involving survey and object detection. Our end-goal is developing, for a swarm
of quadcopters, the ability to autonomously make decisions based on shared
information.

The outline of the paper is as follows. In Section 2 we discuss the multi-
agent autonomous agreement. Agents are considered as particles and the
methodology is based on graph theory. Section 3 introduces the dynamics of
quadcopters and their navigation. In Section 4 we approximate our particle
multi-agent trajectories by implementable ones satisfying the quadcopter’s
dynamic equations of motion.

It is a genuine pleasure to have our article appear in this special issue that
honors Professor Bloch. The depth and breadth of his prodigious body of work
inspires and enriches each of us. On those occasions when the authors feel
like explorers venturing into uncharted territory, we often find a flag
marking Tony’s prior visitation. Indeed, Tony’s contributions to the subject
of geometric control of mechanical systems are keenly felt in this paper. His
recent work on multi-agent formations with obstacle and inter-agent avoidance
in [13] aims toward a similar goal as ours. Although the mathematical approach
taken differs from the one here, some of his work in [29] does apply a graph
theoretic methodology. In addition, the rigid body dynamics that we highlight
in the second half of the paper is core to much of his research ([3, 5, 4,
30]). With sincere appreciation and gratitude, the authors offer Professor
Bloch their very best wishes.

##  2\. Multi-Agent Autonomous Agreement

Suppose a mission with multiple agents requires that at some point the
vehicles come together at a single location, for instance at the end of a
mission. With batteries possibly running low, a main objective is to conserve
power and rely on short-range communication. Autonomous agreement is one of
the fundamental problems in multi-agent coordination, where a collection of
agents are to agree on a joint state value.

###  2.1. Agreement Protocol

Our multi-agent systems will be viewed as networks, as described in this
section. Our primary goal is to highlight the existence of an intricate
relationship between the convergence properties of the so-called agreement
protocol on one hand, and the structure of the underlying agent network
connections on the other.

######  Definition 2.1.

We define a network of agents to be a finite collection of agents possessing
some inter-connectivity via relative information exchange links. The network
is modeled by a graph, denoted $G(V,E)$, where $V$ is the set of vertices and
$E$ is the set of edges. The agents correspond to the vertices of $G(V,E)$ and
the relative information exchange links correspond to the edges.

In the sequel we denote by $V=\\{v_{1},v_{2},\dots,v_{n}\\}$ the vertices
corresponding to the $n$ agents, labeled $1,2,\dots,n$. The set of edges is a
subset $E\subset V\times V$; each edge has the form
$(v_{i},v_{j})=v_{i}v_{j}$. We assume $v_{i}v_{j}\in E$ implies $v_{j}v_{i}\in
E$ and that $v_{i}v_{i}\notin E$.

We shall enhance this initial model to accommodate time-varying weighted
networks, motivated by the desire to account for the strength of
communications between agents, or possibly to take into account how the
strength of communications evolves in the rendezvous problem as distances
between agents change. Interestingly, suitable assumptions on the weights
imply the agreement value is an invariant, and only the trajectories leading
to it differ.

The agreement protocol assumes the rate of change of each agent’s state is
governed by a sum of relative states with respect to a subset of neighboring
agents. Assume the scalar function $q^{i}(t)$ denotes the state of agent $i$
for $i=1,\dots,n$, and define $q(t)=(q^{1}(t),\cdots,q^{n}(t))^{T}$. It
follows that

$\dot{q}^{i}(t)=-\sum_{j\in N(i)}\big{(}q^{i}(t)-q^{j}(t)\big{)},\qquad 1\leq
i\leq n,$ (1)

where $N(i)$ is the neighborhood of a vertex $v_{i}$, i.e., the collection of
adjacent vertices $N(i)=\\{v_{j}\in V\mid v_{i}v_{j}\in E\\}$. If we adopt the
notion that adjacency is symmetric, i.e., $v_{i}v_{j}\in E\iff v_{j}v_{i}\in
E$, then we say the network is undirected and we can encapsulate the above
system in a single matrix equation:

$\dot{q}(t)=-L(G)q(t),$ (2)

where $L(G)\in\mathbb{R}^{n\times n}$ is called the graph Laplacian. The graph
Laplacian is uniquely determined by the structure of $G(V,E)$.

######  Definition 2.2.

We refer to Equation (2) as the _agreement dynamics_ , the _agreement
protocol_ , or the _consensus protocol_.

###  2.2. Consensus Joint Value

The graph Laplacian has desirable properties when $G(V,E)$ is a connected,
undirected network. In the rest of this section we will always assume $G(V,E)$
is a connected, undirected graph unless specified otherwise.

######  Proposition 1.

Assume $G(V,E)$ is a connected, undirected graph. The graph Laplacian $L(G)$
is a symmetric positive semidefinite matrix with eigenvalues
$0=\lambda_{1}<\lambda_{2}\leq\dots\leq\lambda_{n}$. Moreover, the
$n$-dimensional vector of all ones, $\vec{1}$, is an eigenvector corresponding
to the zero eigenvalue $\lambda_{1}$, i.e., $L(G)\vec{1}=\vec{0}$ .

###### Proof.

This is a standard result and can be found in [22, 27, 28]. ∎

Let $\Lambda={\rm diag}(0,\lambda_{2},\dots,\lambda_{n})$ and let
$U=[\hat{u}_{1}\ \hat{u}_{2}\ \dots\ \hat{u}_{n}]$ (the hat notation indicates
a unit vector) be an $n\times n$ matrix consisting of orthonormal eigenvectors
of $L(G)$ corresponding to the eigenvalue ordering above, i.e., $U^{T}\,U=I$
and $L(G)\,U=U\,\Lambda$. Since we can factor the Laplacian as $L(G)=U\Lambda
U^{T}$, the solution to Equation (2) with initial condition $q(0)=q_{0}$ is
given by

$q(t)=e^{-L(G)t}q_{0}=e^{-U\Lambda U^{T}t}q_{0}=Ue^{-\Lambda t}U^{T}q_{0},$

which can be written as:

$q(t)=\sum_{i=1}^{n}e^{-\lambda_{i}t}(\hat{u}^{T}_{i}q_{0})\hat{u}_{i}.$ (3)

######  Definition 2.3.

Given a set of $n$ multi-agents, its _agreement set_
$\mathcal{A}\subset\mathbb{R}^{n}$ is defined as the subspace
$span\\{\vec{1}\\}$. A _consensus joint state value_ is an element of
$\mathcal{A}$, and the common value of its (identical) components is called
its _agreed state value_.

The next proposition shows that when the state of a network of multi-agents
satisfies the agreement protocol, the state converges to the consensus joint
state value whose agreed state value is obtained by averaging the state’s
initial components.

######  Proposition 2.

Let $G(V,E)$ be a network of $n$ agents whose state $q(t)$ satisfies the
agreement protocol $\dot{q}=-L(G)q$. Then
$\lim_{t\rightarrow\infty}q(t)=\alpha\vec{1}$, where
$\alpha=\frac{\vec{1}^{T}q(0)}{n}\in\mathbb{R}$ is the average of the initial
states. In addition, since $\lambda_{2}$ is $L(G)$’s smallest nonzero
eigenvalue, its magnitude will determine the rate of convergence to
$\alpha\vec{1}$ ([28, 22]).

###### Proof.

Since $\lambda_{1}=0$ and $\lambda_{i}>0$ for $i\geq 2$ we have, from Equation
(3), that $q(t)$ converges to the consensus joint state value
$(\hat{u}^{T}_{1}q_{0})\hat{u}_{1}=\frac{\vec{1}^{T}q_{0}}{n}\vec{1}=\alpha\vec{1}\in\mathcal{A}$
as $t\rightarrow\infty$, where
$\frac{\vec{1}^{T}q_{0}}{n}=\alpha\in\mathbb{R}$ is its agreed state value.
Moreover, Equation (3) implies that $e^{\lambda_{2}t}(q^{i}(t)-\alpha)$
converges, as $t\rightarrow\infty$, for all $1\leq i\leq n$. ∎

The point
$\displaystyle\lim_{t\rightarrow\infty}q(t)=\frac{\vec{1}^{T}q_{0}}{n}\vec{1}=\alpha\vec{1}\in\mathbb{R}^{n}$
is precisely the orthogonal projection of $q_{0}$ onto the agreement subspace
and therefore minimizes the quantity $\lVert q-q_{0}\rVert$ over all possible
$q\in\mathcal{A}$, where $\lVert\cdot\rVert$ is the standard Euclidean
distance. Furthermore, $\frac{\vec{1}^{T}q(t)}{n}$ is a constant of motion
since

$\displaystyle\frac{d}{dt}\left(\frac{\vec{1}^{T}q(t)}{n}\right)=$
$\displaystyle\frac{\vec{1}^{T}\dot{q}(t)}{n}=\frac{\vec{1}^{T}(-L(G)q(t))}{n}=-\frac{q^{T}(t)L^{T}(G)\vec{1}}{n}$
(4) $\displaystyle=$
$\displaystyle-\frac{q^{T}(t)L(G)\vec{1}}{n}=-\frac{q^{T}(t)\vec{0}}{n}=0$ (5)

for arbitrary values of $t$. Hence the agreed state value
$\alpha=\frac{\vec{1}^{T}q_{0}}{n}$ is completely determined at $t=0$.

If the state of each agent is a vector quantity rather than a scalar value,
the agreement protocols for the components can still be written in a single
matrix equation. Indeed, assume we desire agreement on $r$ scalar values. This
means that to each agent $i$ we associate a vector value
$q^{i}(t)=(q_{1}^{i}(t),\cdots,q_{r}^{i}(t))$ and we let $Q$ be the $n\times
r$ matrix whose rows are given by the $q^{i}$,

$Q=\left[\begin{array}[]{ccc}q_{1}^{1}&\cdots&q_{r}^{1}\\\
q_{1}^{2}&\cdots&q_{r}^{2}\\\ \vdots&\vdots&\vdots\\\
q_{1}^{n}&\cdots&q_{r}^{n}\\\ \end{array}\right]$ (6)

and we can write the agreement dynamics as:

$\dot{Q}(t)=-L(G)Q(t).$ (7)

We introduce $q^{i}(0)=(q^{i}_{1}(0),\cdots,q^{i}_{r}(0))$ the initial row
vector for the state of agent $i$ and, more importantly,
$q_{k}(0)=(q_{k}^{1}(0),\cdots,q_{k}^{n}(0))^{T}$ the initial column vector
for a given state value over all agents. Since $L(G)=U\Lambda U^{T}$, we have
that the solution of (7) is given by: $Q(t)=e^{-L(G)t}Q(0)=e^{-U\Lambda
U^{T}t}Q(0)=Ue^{-\Lambda t}U^{T}Q(0)$, which can be written as:

$Q(t)=\Bigg{[}\sum_{i=1}^{n}e^{-\lambda_{i}t}\hat{u}_{i}^{T}q_{1}(0)\hat{u}_{i}\,\,\cdots\,\,\sum_{i=1}^{n}e^{-\lambda_{i}t}\hat{u}_{i}^{T}q_{r}(0)\hat{u}_{i}\Bigg{]}$
(8)
$=[e^{-\lambda_{1}t}\hat{u}_{1}^{T}q_{1}(0)\hat{u}_{1}\,\,\cdots\,\,e^{-\lambda_{1}t}\hat{u}_{1}^{T}q_{r}(0)\hat{u}_{1}]+\sum_{i=2}^{n}[e^{-\lambda_{i}t}\hat{u}_{i}^{T}q_{1}(0)\hat{u}_{i}\,\,\cdots\,\,e^{-\lambda_{i}t}\hat{u}_{i}^{T}q_{r}(0)\hat{u}_{i}]$
$=\Bigg{[}\frac{\vec{1}^{T}q_{1}(0)}{n}\vec{1}\,\,\cdots\,\,\frac{\vec{1}^{T}q_{r}(0)}{n}\vec{1}\Bigg{]}+\sum_{i=2}^{n}[e^{-\lambda_{i}t}\hat{u}_{i}^{T}q_{1}(0)\hat{u}_{i}\,\,\cdots\,\,e^{-\lambda_{i}t}\hat{u}_{i}^{T}q_{r}(0)\hat{u}_{i}]$
$=\Big{[}\alpha_{1}\vec{1}\,\,\cdots\,\,\alpha_{r}\vec{1}\Big{]}+\sum_{i=2}^{n}[e^{-\lambda_{i}t}\hat{u}_{i}^{T}q_{1}(0)\hat{u}_{i}\,\,\cdots\,\,e^{-\lambda_{i}t}\hat{u}_{i}^{T}q_{r}(0)\hat{u}_{i}]$
(9)

where $\frac{\vec{1}^{T}q_{k}(0)}{n}=\alpha_{k}\in\mathbb{R}$ is the agreed
state value for each scalar state, $1\leq k\leq r$.

Note that here we assume the states $q_{k}$, measured by the agents, are
shared over the same network $G(V,E)$. Hence the states share a common graph
Laplacian, $L(G)$. In more general situations each state variable could have a
distinct network, leading to an agreement protocol of the form:
$\dot{x}(t)=-L(G_{x})x(t),\dots,\dot{z}(t)=-L(G_{z})z(t)$.

###  2.3. Example:

Consider a collection of $n$ agents required to meet at a single location, not
given in advance, and the agents do not have access to their global positions.
Rather, all they can measure is their relative distance with respect to their
neighbors and they have to agree on their spatial coordinates $x,y$ and $z$.
In this case, $r=3$ and
$q^{i}(t)=(q^{i}_{1}(t),q^{i}_{2}(t),q^{i}_{3}(t))^{T}=(x^{i}(t),y^{i}(t),z^{i}(t))^{T}\in\mathbb{R}^{3}$
with $Q\in\mathbb{R}^{n\times 3}$ as an example of the agreement protocol in
action. By executing the agreement protocol, the convergence to a unique
rendezvous point $(\alpha_{x},\alpha_{y},\alpha_{z})$ is guaranteed for any
arbitrary initial conditions:
$(x^{1}(0),y^{1}(0),z^{1}(0)),\dots,(x^{n}(0),y^{n}(0),z^{n}(0))$ (i.e., any
initial spatial configuration for the $n$ agents) so long as $G(V,E)$ is a
connected graph.

###  2.4. Multi-agent Trajectories

In this section, $\hat{e}_{j}$ denotes the standard unit vector of all $0$’s
with a single $1$ in the $j$th component. We use $u_{1}$ to denote the vector
$\vec{1}$, which is an eigenvector of $L(G)$ with eigenvalue $\lambda_{1}=0$.
With this notation, $\hat{u}_{1}=\dfrac{\vec{1}}{\sqrt{n}}$.

######  Lemma 2.4.

Assume agent $j$ is connected to every other agent, and set $u_{1}=\vec{1}$.
Then, $u_{1}-n\hat{e}_{j}$ is an eigenvector for $L(G)$ with corresponding
eigenvalue $\lambda_{n}=n$:

$L(G)(u_{1}-n\hat{e}_{j})=n(u_{1}-n\hat{e}_{j}).$ (10)

###### Proof.

We have
$L(G)(u_{1}-n\hat{e}_{j})=L(G)u_{1}-nL(G)\hat{e}_{j}=0-nL(G)_{j}=-nL(G)_{j}$,
where $L(G)_{j}$ is the $j$th column of the Laplacian $L(G)$. But since
$G(V,E)$ has the $j$th node with $n-1$ connections we know
$L(G)_{j}=(-1,-1,\dots,n-1,\dots,-1,-1)^{T}=n\hat{e}_{j}-u_{1}$, i.e., all
$-1$’s except an $n-1$ in the $j$th component. Hence
$-nL(G)_{j}=n(u_{1}-n\hat{e}_{j})$ as claimed.

∎

######  Definition 2.5.

We say that vertex $v_{i}$ is _fully connected_ if there exists an edge
between $v_{i}$ and every other vertex $v_{j}\in V,j\neq i$. Therefore a fully
connected vertex has $n-1$ edges.

Recall that $\hat{u}_{1}=\frac{\vec{1}}{\sqrt{n}}$ and $\lambda_{1}=0$.
Furthermore, it is proven in [1] that the eigenvalues of $L(G)$ are less than
or equal to $n$ and the multiplicity of the eigenvalue $\lambda=n$ is equal to
the number of fully connected vertices.

######  Theorem 2.6.

Suppose vertex $v_{j}$ is fully connected and that $q(t)$ satisfies the
agreement protocol $\dot{q}=-L(G)q$. Then $q^{j}(t)$ is a straight line.

###### Proof.

We know that $\displaystyle\sum_{i=1}^{n}q^{i}=n\alpha$ is a constant. We also
know that:

$\dot{q}^{j}=-n\,q^{j}+\sum q^{i}=-n(q^{j}-\alpha).$ (11)

This simple ODE (Newton’s Law of Cooling) implies that

$q^{j}-\alpha=(q^{j}(0)-\alpha)e^{-nt}$ (12)

which is a straight line. ∎

Note that the proof is the same regardless of whether $q^{j}$ is a scalar or
row $j$ of an $n\times r$ state matrix $Q(t)$, provided it satisfies the
agreement protocol $\dot{Q}=-L(G)Q$.

####  2.4.1. Rendezvous Mission, Unweighted Network

In this section we illustrate some of our prior results with a rendezvous
mission between four agents. We let $q^{i}(t)=(x^{i}(t),y^{i}(t),z^{i}(t))$
denote the position of agent $i$ in $\mathbb{R}^{3}$. We assume the agents to
be initially distributed as follows:

$q^{1}(0)=(4,17,24),\;q^{2}(0)=(18,10,32),\;q^{3}(0)=(15,10,26),\;q^{4}(0)=(4,2,35).$

In this simulation we assume the communication network between the agents is
fixed and unweighted. From our algorithm, the vector of agreed state values
(the agreed rendezvous position) is given by:

$q^{*}=(10.25,9.75,29.25).$ (13)

In Figure 1 we illustrate how two different communication networks produce
different trajectories, even though the agreed rendezvous position is the same
for each network since it depends only on the initial spatial positions of the
agents. In the first scenario (solid curve), we assume the following: agent 1
is connected to agent 2; agent 2 is connected to agents 1, 3 and 4; agent 3 is
connected to agents 2, 4; agent 4 is connected to agents 2, 3. This is
represented by the following graph Laplacian:

$L_{1}=\left(\begin{array}[]{cccc}1&-1&0&0\\\ -1&3&-1&-1\\\ 0&-1&2&-1\\\
0&-1&-1&2\end{array}\right).$ (14)

In the second scenario (dashed curve), we assume: agent 1 is connected to
agents 3, 4; agent 2 is connected to agent 3; agent 3 is connected to agents
1, 2; agent 4 is connected to agent 1:

$L_{2}=\left(\begin{array}[]{cccc}2&0&-1&-1\\\ 0&1&-1&0\\\ -1&-1&2&0\\\
-1&0&0&1\end{array}\right).$ (15)

In the first scenario agent 2 moves along a straight line since it is
connected to all other agents.

Figure 1. Rendezvous Missions with Unweighted Network. Displays two different
communication network scenarios for a 4-agent rendezvous mission. Agreement
positions coincide but trajectories differ.

In Figure 2 we compare the trajectories for agent 1 for each spatial
coordinate. The eigenvalues for scenario 1 are given by $\\{0,1,3,4\\}$ and
the eigenvalues for scenario 2 are given by $\\{0,0.586,2,3.414\\}$.

Figure 2. Comparison of the $x,y,z$-motions for agent 1 for the scenarios of
Fig. 1 corresponding to rendezvous missions with Unweighted Network.

###  2.5. Weighted Network

We can imagine a scenario where it is desirable to account for the
communications strength between vehicles, or for how that strength evolves
with time. Interestingly, the agreed state value is an invariant, and only the
trajectories leading to it are affected [28]. Note that not all extensions are
immediately straightforward, however, since the explicit solution may be
difficult to write down. Not all properties can be expected to carry over; for
example, the linear trajectory proof relied on cancellation due to the integer
eigenvalues.

We now assume that the edges are weighted, i.e., edges carry a numerical
weight. If the graph $G_{w}$ is a weighted graph, then the off-diagonal
entries of the Laplacian $L(G_{w})$ will not be limited to only $0$ and $-1$.
Instead, the $-1$’s corresponding to edges in the graph would be replaced with
$-w_{ij}=-w_{ji}$ to account for the weights. The diagonal entries will once
again be the absolute value of the row or column sum of the non-diagonal
entries. Equations (1) and (2) therefore become:

$\displaystyle\dot{q}^{i}(t)$ $\displaystyle=-\sum_{j\in
N(i)}w_{ij}\big{(}q^{i}(t)-q^{j}(t)\big{)}$ (16) $\displaystyle\dot{q}(t)$
$\displaystyle=-L(G_{w})q(t).$ (17)

The proof that the consensus joint state value is given by
$\alpha=\frac{\vec{1}q_{0}}{n}\vec{1}$ comes from the same logic as before
since $L(G_{w})$ is a symmetric positive semidefinite matrix with eigenvalues
$0=\lambda_{1}<\lambda_{2}\leq\dots\leq\lambda_{n}$ and
$L(G_{w})\vec{1}=\vec{0}$. This is a consequence of the fact that weighting
preserves the properties that $G_{w}(V,E)$ is connected and undirected.
Moreover, each row and column of $L(G_{w})$ still sums to zero regardless of
whether or not the edges of $G_{w}(V,E)$ are weighted. Hence, the values of
$\lambda_{2},\dots\lambda_{n}$ and $\hat{u}_{2}\dots\hat{u}_{n}$, will change,
but we still have $\lambda_{1}=0$ and $\hat{u}_{1}=\frac{\vec{1}}{\sqrt{n}}$.
Finally, in full generality:

######  Theorem 2.7.

Consider the dynamical system:

$\dot{q}^{i}(t)=-\sum_{j=1}^{n}w_{ij}(t)[q^{i}(t)-q^{j}(t)]\quad(1\leq i\leq
n)$ (18)

with initial conditions $q(0)=(q^{1}(0),\cdots q^{n}(0))^{T}$.

For each $i,j$, assume $w_{ij}:[0,\infty)\to[0,\infty)$ and
$w_{ij}(t)=w_{ji}(t)$.

Define the Laplacian matrix $L(G_{w})(t)$ by

$L(G_{w})(t)=\left[\begin{array}[]{cccc}\sum_{j=1}^{n}w_{1j}(t)&-w_{12}(t)&\dots&-w_{1n}(t)\\\
-w_{12}(t)&\sum_{j=1}^{n}w_{2j}(t)&\dots&-w_{2n}(t)\\\
\vdots&\vdots&\vdots&\vdots\\\
-w_{1n}(t)&-w_{2n}(t)&\dots&\sum_{j=1}^{n}w_{nj}(t)\\\ \end{array}\right].$
(19)

Then the system of first order equations can be rewritten as:
$\dot{q}(t)=-L(G_{w})(t)q(t)$. If the rank of $L(G_{w})(t)$ is $n-1$ for all
$t$ then we have:

$\lim_{t\rightarrow\infty}q(t)=\bigg{[}\frac{1}{n}\sum_{i=1}^{n}q^{i}(0)\bigg{]}\vec{1}.$

###### Proof.

By assumption, the matrix $L(G_{w})$ is real and symmetric, therefore
Hermitian. Furthermore, it is diagonally dominant with non-negative diagonal
entries. Therefore it is positive semi-definite. As for the unweighted case,
$\lambda_{1}=0$ with corresponding eigenvector $\vec{1}$. Since we assume the
graph to be connected throughout the motion, the rank of $L$ is always $n-1$
and the multiplicity of the zero eigenvalue is one. We have for any time $t$:
$0=\lambda_{1}<\lambda_{2}(t)\leq\dots\leq\lambda_{n}(t)$ with a corresponding
set of orthonormal eigenvectors
$\hat{u}_{1}=\frac{1}{\sqrt{n}}\vec{1},\hat{u}_{2}(t),\dots,\hat{u}_{n}(t)$.

Let $q^{*}=\bigg{[}\frac{1}{n}\sum_{i=1}^{n}q^{i}(0)\bigg{]}\vec{1}$. Even if
the $w_{ij}:[0,\infty)\to[0,\infty)$ are time-varying, we can show with
initial conditions $q(0)$ that the trajectories $q(t)\to q^{*}$ as
$t\to\infty$. A useful observation is that

$\dot{q}(t)=-L(G_{w})(t)q(t)\quad\Rightarrow\quad
q(t)-q(0)=-\int_{0}^{t}L(G_{w})(s)q(s)ds$
$\Rightarrow\quad\vec{1}^{T}\big{(}q(t)-q(0)\big{)}=-\vec{1}^{T}\bigg{(}\int_{0}^{t}L(G_{w})(s)q(s)ds\bigg{)}$
$=-\int_{0}^{t}\vec{1}^{T}L(s)q(s)ds=-\int_{0}^{t}\vec{0}^{T}q(s)ds=-\int_{0}^{t}0ds=0,$

so that we have $\vec{1}^{T}q(t)=\vec{1}^{T}q(0)$ for all $t\in[0,\infty)$.

Let $V:\mathbb{R}^{n}\to\mathbb{R}$ be defined by:

$q\mapsto V(q)=\frac{1}{2}\|q-q^{*}\|^{2}=\frac{1}{2}(q-q^{*})^{T}(q-q^{*}).$

Clearly, we have:

  1. (1)

$V(q)\geq 0$,

  2. (2)

$V(q)=0\iff q=q^{*}$,

  3. (3)

$\dot{V}(q)=-(q-q^{*})^{T}L(G_{w})(q-q^{*})\leq 0$

since $-v^{T}L(G_{w})v\leq 0$ for any $v\in\mathbb{R}^{n}$ since $L(G_{w})$ is
positive semi-definite. Clearly $q=q^{*}$ implies $\dot{V}(q)=0$, but also
$\dot{V}(q)=0$ implies $q=q^{*}$. Let
$\mathcal{A}=span\\{\vec{1}\\}=\mathcal{N}(L(G_{w}))$ represent the null space
of $L(G_{w})$ and suppose $(q-q^{*})^{T}L(G_{w})(q-q^{*})=0$. This implies
$q-q^{*}\in\mathcal{A}$. But we know that $\vec{1}^{T}q(t)=\vec{1}^{T}q(0)$,
and therefore $q-q^{*}\perp\mathcal{A}$ since

$\vec{1}^{T}(q-q^{*})=\vec{1}^{T}\big{(}q-\frac{\vec{1}^{T}q(0)}{n}\vec{1}\big{)}=\vec{1}^{T}q-\frac{\vec{1}^{T}q(0)}{n}\vec{1}^{T}\vec{1}=\vec{1}^{T}q-\vec{1}^{T}q(0)=0.$
(20)

Finally, $q-q^{*}\in\mathcal{A}$ and $q-q^{*}\perp\mathcal{A}$ imply
$q-q^{*}=0\Rightarrow q=q^{*}$. We conclude:

$\dot{V}(q)=-(q-q^{*})^{T}L(G_{w})(q-q^{*})=0\iff q=q^{*}.$ (21)

Therefore 1, 2, 3 and Equation (21) imply that $V(q)$ is a non-negative
function which is strictly decreasing for all trajectories $q(t)\neq q^{*}$
hence $q(t)\to q^{*}$ as $t\to\infty$, i.e.,

$\lim_{t\rightarrow\infty}q(t)=\bigg{[}\frac{1}{n}\sum_{i=1}^{n}q_{i}(0)\bigg{]}\vec{1}$

as claimed.

∎

####  2.5.1. Rendezvous Mission, Weighted Network

Here we imagine a scenario taking into account the strength of communications
between agents. As previously stated, the agreement value is an invariant, and
only the trajectories leading to it will differ compared to not taking the
strength of the signal into account. We assume four agents to be initially
distributed as follows:

$q^{1}(0)=(2,16,29),\;q^{2}(0)=(14,20,31),\;q^{3}(0)=(10,1,25),\;q^{4}(0)=(10,14,36).$

The agreement rendezvous position is given by:

$q^{*}=(9,12.75,30.25).$ (22)

In Figure 3 we illustrate how a weighted network impacts the trajectories. For
both scenarios we assume: agent 1 is connected to agent 2; agent 2 is
connected to agents 1, 3 and 4; agent 3 is connected to agents 2, 4; agent 4
is connected to agents 2, 3 (as in Section 2.4.1):

$L_{1}=\left(\begin{array}[]{cccc}1&-1&0&0\\\ -1&3&-1&-1\\\ 0&-1&2&-1\\\
0&-1&-1&2\end{array}\right).$ (23)

In the second scenario (dashed curves), we assume a weighted network with
weights due to initial distance:

$L_{2}=\left(\begin{array}[]{cccc}12.8062&-12.8062&0&0\\\
-12.8062&41.9036&-20.3224&-8.7750\\\ 0&-20.3224&37.3518&-17.0294\\\
0&-8.7750&-17.0294&25.8044\end{array}\right).$ (24)

It can be observed that for the weighted network, even though agent 2 is
connected to all other agents, it does not move along a straight line anymore.

Figure 3. Comparison between trajectories on rendezvous missions with
unweighted and weighted networks. The solid curves represents the trajectories
for the unweighted network and the dashed ones for the weighted network. The
two scenarios converge to the same agreement value.

In Figure 4 we compare the trajectories for agent 1 for each of the
coordinates of motion. The eigenvalues for the unweighted network are given by
$\\{0,1,3,4\\}$ and the eigenvalues for the weighted network are given by
$\\{0,13.060,42.248,62.558\\}$.

Figure 4. Comparison of the $x,y,z$-motions for agent 1 for the rendezvous
missions displayed in Figure 3. Observe that each component converges much
more rapidly for the weighted network, corresponding to the difference in
their $\lambda_{2}$ values: $13.060>1$.

####  2.5.2. Rendezvous Mission, Time-varying Weighted Network

The goal of this rendezvous mission is to compare a fixed network determined
by initial proximity versus a time-varying weighted network, where the weights
are the inter-agent distances that vary in time, and also edges are added as
the preset proximity threshold is reached.

We assume four agents to be initially distributed as follows:

$q^{1}(0)=(16,5,36),\;q^{2}(0)=(19,19,29),\;q^{3}(0)=(12,16,33),\;q^{4}(0)=(14,1,26).$

The agreement rendezvous position is given by:

$q^{*}=(15.25,10.25,31).$ (25)

In Figure 5 we illustrate how a weighted network impacts the trajectories. For
the fixed network (solid curves) we have: agent 1 is connected to agents 3, 4;
agent 2 is connected to agent 3; agent 3 is connected to agents 1, 3; agent 4
is connected to agent 1:

$L_{1}=\left(\begin{array}[]{cccc}2&0&-1&-1\\\ 0&1&-1&0\\\ -1&-1&2&0\\\
-1&0&0&1\end{array}\right).$ (26)

In the time-varying scenario (dashed lines), the network changes throughout
the trajectories. Here we provide the initial Laplacian and the final one:

$L_{2_{i}}=\left(\begin{array}[]{cccc}23.0375&0&-12.0830&-10.9545\\\
0&8.6023&-8.6023&0\\\ -12.0830&-8.6023&20.6854&0\\\
-10.9545&0&0&10.9545\end{array}\right),$ (27)
$L_{2_{f}}=\left(\begin{array}[]{cccc}0.1465&-0.0523&-0.0482&-0.0460\\\
-0.0523&0.1501&-0.0450&-0.0528\\\ -0.0482&-0.0450&0.1457&-0.0525\\\
-0.0460&-0.0528&-0.0525&0.1513\end{array}\right).$ (28)

Initially the connectivity network is the same as for the fixed network but
with weights, however as the agents get closer it changes, and all agents are
connected in the end. The trajectories corresponding to the time-varying
network here would be easier to approximate for a real quadcopter since they
are “straighter” (the sudden jump in information leads to straighter paths
toward the rendezvous, however in this case, the weights prevent the perfect
cancellation that would lead to true linearity.)

Figure 5. Rendezvous missions comparing and unweighted network (solid curves)
to a time-varying weighted one (dashed curves). They both agree on the
consensus joint value. The points where the dashed curves diverge from the
solid ones correspond to the addition of edges.

In Figure 6 we compare the trajectories for agent 1 for each coordinate of
motion. The eigenvalues for the unweighted network are given by
$\\{0,1,3,4\\}$ and the eigenvalues for the weighted network are initially
given by $\\{0,6.234,19.385,37.655\\}$, and near the rendezvous agreement are
given by $\\{0,0.190,0.196,0.207\\}$.

Figure 6. Comparison of the $x,y,z$-motions for agent 1 for the rendezvous
mission of Figure 5. The components of the trajectory corresponding to the
time-varying network clearly converge more rapidly.

##  3\. Dynamics of quadcopters

In order to successfully implement multi-agent motion planning, we need
dynamically sound equations of motion governing the movement of each drone. To
this end, in this section we will derive such equations under certain
realistic symmetry assumptions regarding the shape of the drone body.
Utilizing the geometric formalism of [10], we present these equations in both
first-order and second-order forms, and provide a differential geometric
version in terms of the Levi-Civita connection on our configuration space
equipped with a natural Riemannian metric. Finally, we give the equations
explicitly in coordinates. Note that the notation in this section is
independent from that in the previous section.

###  3.1. Setup

In this subsection we provide the basic physical notation and standard
derivation of the equations of motion for a rigid body in the body-fixed
frame. A rigid body can model a number of vehicles in space or in a fluid,
such as satellites or underwater vehicles, which can be stabilized by spinning
rotors ([5, 6]). Here, we assume our drone is nothing more than a rigid body
in a fluid; later we will make certain additional assumptions explicit. We
closely follow the derivations in [16, 31].

The configuration space for a rigid body moving in a fluid ([3]) is

$Q=SE(3)\cong\mathbb{R}^{3}\times SO(3).$

Here we equip $Q$ with coordinates $(b,R)$ where $b\in\mathbb{R}^{3}$
represents the position of the center of mass of the body in space and
$R\in\mathbb{R}^{3\times 3}$ is a rotation matrix representing the orientation
of the body aligned with the principal axes of inertia. This coordinate system
corresponds to some inertial frame of reference.

We regularly use the notation of the hat map, the Lie algebra isomorphism

$\hat{}:(\mathbb{R}^{3},\times)\to(\mathfrak{so}(3),[,])$

given by $\hat{y}z=y\times z$, equivalently
$\hat{y}=\begin{pmatrix}0&-y_{3}&y_{2}\\\ y_{3}&0&-y_{1}\\\
-y_{2}&y_{1}&0\end{pmatrix}$.

First, choose a non-inertial body-fixed frame and set the following notation:

$\displaystyle v$ $\displaystyle=\text{translational velocity in body frame}$
$\displaystyle\Omega$ $\displaystyle=\text{angular velocity in body frame}$
$\displaystyle p$ $\displaystyle=\text{translational momentum in inertial
frame}$ $\displaystyle\pi$ $\displaystyle=\text{angular momentum in inertial
frame}$ $\displaystyle P$ $\displaystyle=\text{translational momentum in body
frame}$ $\displaystyle\Pi$ $\displaystyle=\text{angular momentum in body
frame}.$

We seek to derive the dynamic equations of motion for the system in the body
frame. Here this will take the form of a system of ODEs in the variables
$(b,R,v,\Omega)\in TQ$. Our kinematic equations are:

$\displaystyle\dot{b}$ $\displaystyle=Rv$ (29) $\displaystyle\dot{R}$
$\displaystyle=R\hat{\Omega}.$ (30)

We transform momenta between the inertial and body frames according to:

$\displaystyle p$ $\displaystyle=RP$ (31) $\displaystyle\pi$
$\displaystyle=R\Pi+\hat{b}p$ (32)

Begin by differentiating equations (31) and (32):

$\displaystyle\dot{p}$ $\displaystyle=\dot{R}P+R\dot{P}$ (33)
$\displaystyle\dot{\pi}$ $\displaystyle=\dot{R}\Pi+R\dot{\Pi}+\dot{b}\times
p+b\times\dot{p}.$ (34)

To rewrite these equations in a single frame, we must introduce additional
quantities. In the inertial frame, let

$\displaystyle f_{i}$ $\displaystyle=\text{external forces},\quad i=1,\dots,k$
$\displaystyle\tau_{i}$ $\displaystyle=\text{external torques},\quad
i=1,\dots,l,$

so that

$\displaystyle E_{F}$ $\displaystyle=\sum_{i=1}^{k}R^{t}f_{i}=\text{total
external force in body frame}$ (35) $\displaystyle E_{T}$
$\displaystyle=\sum_{i=1}^{l}R^{t}\tau_{i}=\text{total external torque in body
frame}.$ (36)

Now in the inertial frame, the dynamics are solely due to external forces and
torques:

$\displaystyle\dot{p}$ $\displaystyle=\sum_{i=1}^{k}f_{i}$ (37)
$\displaystyle\dot{\pi}$
$\displaystyle=\sum_{i=1}^{k}\hat{x}_{i}f_{i}+\sum_{i=1}^{l}\tau_{i},$ (38)

where $x_{i}$ is the vector from inertial frame origin to line of action of
$f_{i}$.

We then solve (33), (34) for $\dot{P}$ and $\dot{\Pi}$ and simplify using
equations (29), (30), (31), (35), (36), (37), (38):

$\displaystyle\dot{P}$ $\displaystyle=\hat{P}\Omega+E_{F}$ (39)
$\displaystyle\dot{\Pi}$
$\displaystyle=\hat{\Pi}\Omega+\hat{P}v+E_{T}+R^{t}\sum_{i}(x_{i}-b)\times
f_{i}.$ (40)

This gives the evolution of the momenta in the body frame, but the equations
mix momenta and velocities. In order to arrive at equations in terms of solely
momenta or velocities, we must find explicit formulas relating the two. This
can be acheived by the Legendre transform

$P=\frac{\partial T}{\partial v},\qquad\Pi=\frac{\partial T}{\partial\Omega}.$

To make this explicit, we introduce the kinetic energies

$\displaystyle T_{body}$ $\displaystyle=\frac{1}{2}\begin{pmatrix}v\\\
\Omega\end{pmatrix}^{t}\begin{pmatrix}mI_{3}&-m\hat{r}_{C_{G}}\\\
m\hat{r}_{C_{G}}&J_{b}\end{pmatrix}\begin{pmatrix}v\\\ \Omega\end{pmatrix}$
(41) $\displaystyle T_{fluid}$ $\displaystyle=\frac{1}{2}\begin{pmatrix}v\\\
\Omega\end{pmatrix}^{t}\begin{pmatrix}M_{f}&C_{f}^{t}\\\
C_{f}&J_{f}\end{pmatrix}\begin{pmatrix}v\\\ \Omega\end{pmatrix}$ (42)
$\displaystyle T$ $\displaystyle=T_{body}+T_{fluid}$
$\displaystyle=\frac{1}{2}\begin{pmatrix}v\\\
\Omega\end{pmatrix}^{t}\begin{pmatrix}mI_{3}+M_{f}&-m\hat{r}_{C_{G}}+C_{f}^{t}\\\
m\hat{r}_{C_{G}}+C_{f}&J_{b}+J_{f}\end{pmatrix}\begin{pmatrix}v\\\
\Omega\end{pmatrix}$ (43)

where $m$ is the mass of the body, $r_{C_{G}}$ is vector from center of
gravity to body frame origin, $J_{b}$ is the body inertia tensor, $J_{f}$ is
the added mass inertia tensor, $M_{f}$ is the added mass, $C_{f}$ is the added
cross terms, and $I_{3}$ is the $3\times 3$ identity matrix.

Then by differentiating the total kinetic energy (43), we find

$\displaystyle P=\frac{\partial T}{\partial v}$
$\displaystyle=(mI_{3}+M_{f})v+(-m\hat{r}_{C_{G}}+C_{f}^{t})\Omega$ (44)
$\displaystyle\Pi=\frac{\partial T}{\partial\Omega}$
$\displaystyle=(m\hat{r}_{C_{G}}+C_{f})v+(J_{b}+J_{f})\Omega.$ (45)

These equations transform velocities into momenta.

We now assume the body has an especially nice shape, which is realistic for
most drones.

Assumption 1. Assume the body has three planes of symmetry and the principle
axes of inertial coincide with the body frame axes. This implies
$J_{b},J_{f},$ and $M_{f}$ are all diagonal and $C_{f}=0$. When we consider
gravity in Section 3.2, we also assume that this symmetry extends to the mass
density of the body.

Assumption 2. Assume the center of gravity of the body coincides with the
origin of the body frame. This implies $\hat{r}_{C_{G}}=0$.

With these two assumptions, equations (44) and (45) reduce to

$\displaystyle P$ $\displaystyle=Mv$ (46) $\displaystyle\Pi$
$\displaystyle=J\Omega,$ (47)

where $M=mI_{3}+M_{f}$ is the total mass of the system and $J=J_{b}+J_{f}$ is
the moment of inertia of the system.

Now differentiate (46) and (47) and simplify using (39), (40), (46) and (47)
to obtain the derivatives of the velocities expressed solely in terms of
positions and velocities, without momenta:

$\displaystyle m\dot{v}$ $\displaystyle=Mv\times\Omega+E_{F}$ (48)
$\displaystyle J\dot{\Omega}$ $\displaystyle=J\Omega\times\Omega+Mv\times
v+E_{T}+R^{t}\sum_{i}(x_{i}-b)\times f_{i}.$ (49)

Finally, since for drones our fluid is air, we make the following simplifying
assumption.

Assumption 3. Assume the added fluid is negligible. This implies
$M_{f}=J_{f}=0$. Then $J=J_{b}$, and since $M=mI_{3}$ is a scalar matrix, the
fictitious force $Mv\times v=0$.

Combining this with (29), (30), (48), (49), we obtain the complete equations
of motion for a rigid drone in air as a first order system on $TQ$
coordinatized by $(b,R,v,\Omega)$:

$\displaystyle\dot{b}$ $\displaystyle=Rv$ (50) $\displaystyle\dot{R}$
$\displaystyle=R\hat{\Omega}$ (51) $\displaystyle m\dot{v}$
$\displaystyle=mv\times\Omega+E_{F}$ (52) $\displaystyle J\dot{\Omega}$
$\displaystyle=J\Omega\times\Omega+E_{T}+R^{t}\sum_{i=1}^{k}(x_{i}-b)\times
f_{i}.$ (53)

###  3.2. External forces and torques

In this section we explicitly work out the external terms $E_{F}$ and $E_{T}$
for a drone in the body frame. We have four rotors producing thrust in a
quadcopter, with the $i$th rotor possessing angular velocity $\omega_{i}$
producing thrust $t_{i}=K_{r}\omega_{i}^{2}$, where $K_{r}$ is the thrust
coefficient ([26]). Here we are implicitly treating the $\omega_{i}$ as our
controls.

For the external forces, we follow [32]. The force due to drag is

$f_{1}=-\text{diag}(v_{1}|v_{1}|,v_{2}|v_{2}|,v_{3}|v_{3}|)C_{D},$

where $v_{j}$ is the $j$th component of the linear velocity $v$ and $C_{D}$ is
the vector of translational drag coefficients in the body frame. The force due
to gravity is

$f_{2}=-mgR^{t}e_{3}^{I},$

where $g$ is the gravitational constant and the third inertial frame basis
vector $e_{3}^{I}$ represents the opposite direction of gravitational
attraction. Lastly, the force produced by thrust is

$f_{3}=\sum_{i=1}^{4}e_{3}^{B}t_{i},$

where the third body frame basis vector $e_{3}^{B}$ represents the direction
in which each of the motors produce thrust. Thus $E_{F}=f_{1}+f_{2}+f_{3}$.

Ostensibly, each of these forces induces a torque, represented by the last
term $R^{t}\sum_{i=1}^{k}(x_{i}-b)\times f_{i}$ in (53). However, for $f_{1}$
and $f_{2}$, no torque is induced as the line of force passes through the
center of mass of the body due to our symmetry assumptions. The force $f_{3}$,
however, represents the sum of four individual forces located at each of the
four rotors. Each of these forces induces a torque, and their sum gives the
net torque generated by the rotors:

$\tau_{f}=\begin{pmatrix}K_{r}d(\omega_{3}^{2}-\omega_{1}^{2})\\\
K_{r}d(\omega_{4}^{2}-\omega_{2}^{2})\\\
K_{d}\sum_{i=1}^{4}(-1)^{i+1}\omega_{i}^{2}\end{pmatrix},$

where $K_{d}$ is the propeller drag coefficient and $d$ is the distance from
the drone center of mass to the rotation axis of each rotor.

The external torques are computed as in [8, 32]. The torque in the body frame
generated by drag is

$\tau_{1}=-\text{diag}(\Omega_{1}|\Omega_{1}|,\Omega_{2}|\Omega_{2}|,\Omega_{3}|\Omega_{3}|)C_{\tau},$

where $\Omega_{j}$ is the $j$th component of the angular velocity $\Omega$ and
$C_{\tau}$ is the vector of rotational drag coefficients in the body frame.
The torque due to gyroscopic effects is

$\tau_{2}=\sum_{i=1}^{4}\Omega\times(-1)^{i+1}J_{r}(0,0,\omega_{i})^{t},$

where $J_{r}$ is the moment of inertia for a rotor. Assuming the rotor has
nontrivial inertia about only the vertical axis, $J_{r}$ has only one nonzero
entry which we denote by $\bar{J_{r}}$, so that
$J_{r}=\text{diag}(0,0,\bar{J_{r}})$. We then have $E_{T}=\tau_{1}+\tau_{2}$.

We can now give our updated equations of motion for a quadcopter:

$\displaystyle\dot{b}$ $\displaystyle=Rv$ (54) $\displaystyle\dot{R}$
$\displaystyle=R\hat{\Omega}$ (55) $\displaystyle m\dot{v}$
$\displaystyle=mv\times\Omega+f_{1}+f_{2}+f_{3}$ (56) $\displaystyle
J\dot{\Omega}$ $\displaystyle=J\Omega\times\Omega+\tau_{f}+\tau_{1}+\tau_{2}.$
(57)

###  3.3. Coordinate expression

Here we work in coordinates $\eta=(b_{1},b_{2},b_{3},\phi,\theta,\psi)$ on
$Q=SE(3)$ representing the position and orientation of the drone in the
inertial frame. The positions $b_{1},b_{2},b_{3}$ are the standard coordinates
for $\mathbb{R}^{3}$. The angles $\phi,\theta,\psi$ are Tait-Bryan angles,
known in aeronautics as roll, pitch, and yaw respectively, and sometimes
referred to as Euler angles (although not the proper or classical kind). The
rotation angles are given intrinsically by the sequence
$z-y^{\prime}-x^{\prime\prime}$ or extrinsically by the sequence $x-y-z$. Here
we take $(\phi,\theta,\psi)\in(-\pi,\pi)\times(-\pi/2,\pi/2)\times(-\pi,\pi).$
These coordinates on $Q$ induce natural coordinates for the velocities
$(v,\Omega)=(v_{1},v_{2},v_{3},\Omega_{1},\Omega_{2},\Omega_{3})$.

Now let

$R(\eta)=\begin{pmatrix}C_{\psi}C_{\theta}&C_{\psi}S_{\theta}S_{\phi}-S_{\psi}C_{\phi}&C_{\psi}S_{\theta}C_{\phi}+S_{\psi}S_{\phi}\\\
S_{\psi}C_{\theta}&S_{\psi}S_{\theta}S_{\phi}+C_{\psi}C_{\phi}&S_{\psi}S_{\theta}C_{\phi}-C_{\psi}S_{\phi}\\\
-S_{\theta}&C_{\theta}S_{\phi}&C_{\theta}C_{\phi}\end{pmatrix},$

and

$\Theta(\eta)=\begin{pmatrix}1&S_{\phi}T_{\theta}&C_{\phi}T_{\theta}\\\
0&C_{\phi}&-S_{\phi}\\\
0&\frac{S_{\phi}}{C_{\theta}}&\frac{C_{\phi}}{C_{\theta}}\end{pmatrix},$

where $C_{x}=\cos(x),$ $S_{x}=\sin(x)$ and $T_{x}=\tan(x)$. We can write $R$
as the composition of three pure rotations as follows. Let $R_{i,\alpha}$
denote rotation by $\alpha$ radians about the $i$th principal axis of inertia.
Explicitly,

$\displaystyle R_{1,\phi}=\begin{pmatrix}1&0&0\\\ 0&C_{\phi}&-S_{\phi}\\\
0&S_{\phi}&C_{\phi}\end{pmatrix},$ $\displaystyle
R_{2,\theta}=\begin{pmatrix}C_{\theta}&0&S_{\theta}\\\ 0&1&0\\\
-S_{\theta}&0&C_{\theta}\end{pmatrix},$ $\displaystyle
R_{3,\psi}=\begin{pmatrix}C_{\psi}&-S_{\psi}&0\\\ S_{\psi}&C_{\psi}&0\\\
0&0&1\end{pmatrix}.$

Then $R=R_{3,\psi}R_{2,\theta}R_{1,\phi}$.

Now $R$ transforms linear velocity in the body frame to linear velocity in the
inertial frame, and $\Theta$ transforms the body-fixed angular velocity
$(\Omega_{1},\Omega_{2},\Omega_{3})$ into the Euler rate vector
$(\dot{\phi},\dot{\theta},\dot{\psi})$. Thus our kinematic equations (54) and
(55) take the form

$\dot{\eta}=\begin{pmatrix}R(\eta)&0\\\
0&\Theta(\eta)\end{pmatrix}\begin{pmatrix}v\\\ \Omega\end{pmatrix}.$

Substituting $R$ into equations (56) and (57) and expanding yields the
following.

######  Lemma 3.1.

In coordinates developed above, the equations of motion for a quadcopter take
the form

$\displaystyle\dot{b}_{1}$
$\displaystyle=v_{1}C_{\psi}C_{\theta}+v_{2}(C_{\psi}S_{\theta}S_{\phi}-S_{\psi}C_{\phi})+v_{3}(C_{\psi}S_{\theta}C_{\phi}+S_{\psi}S_{\phi})$
(58) $\displaystyle\dot{b}_{2}$
$\displaystyle=v_{1}S_{\psi}C_{\theta}+v_{2}(S_{\psi}S_{\theta}S_{\phi}+C_{\psi}C_{\phi})+v_{3}(S_{\psi}S_{\theta}C_{\phi}-C_{\psi}S_{\phi})$
(59) $\displaystyle\dot{b}_{3}$
$\displaystyle=-v_{1}S_{\theta}+v_{2}C_{\theta}S_{\phi}+v_{3}C_{\theta}C_{\phi}$
(60) $\displaystyle\dot{\phi}$
$\displaystyle=\Omega_{1}+\Omega_{2}S_{\phi}T_{\theta}+\Omega_{3}C_{\phi}T_{\theta}$
(61) $\displaystyle\dot{\theta}$
$\displaystyle=\Omega_{2}C_{\phi}-\Omega_{3}S_{\phi}$ (62)
$\displaystyle\dot{\psi}$
$\displaystyle=\Omega_{2}\frac{S_{\phi}}{C_{\theta}}+\Omega_{3}\frac{C_{\phi}}{C_{\theta}}$
(63) $\displaystyle\dot{v}_{1}$
$\displaystyle=v_{2}\Omega_{3}-v_{3}\Omega_{2}-\frac{1}{m}v_{1}|v_{1}|C_{D_{1}}+gS_{\theta}$
(64) $\displaystyle\dot{v}_{2}$
$\displaystyle=v_{3}\Omega_{1}-v_{1}\Omega_{3}-\frac{1}{m}v_{2}|v_{2}|C_{D_{2}}-gC_{\theta}S_{\phi}$
(65) $\displaystyle\dot{v}_{3}$
$\displaystyle=v_{1}\Omega_{2}-v_{2}\Omega_{1}+\frac{1}{m}\Bigg{(}\sum_{i=1}^{4}t_{i}-v_{3}|v_{3}|C_{D_{3}}\Bigg{)}-gC_{\theta}C_{\phi}$
(66) $\displaystyle\dot{\Omega}_{1}$
$\displaystyle=\frac{1}{J_{1}}\left[(J_{2}-J_{3})\Omega_{2}\Omega_{3}+\bar{J}_{r}\sum_{i=1}^{4}(-1)^{i+1}\Omega_{2}\omega_{i}+K_{r}d(\omega_{3}^{2}-\omega_{1}^{2})-\Omega_{1}|\Omega_{1}|C_{\tau_{1}}\right]$
(67) $\displaystyle\dot{\Omega}_{2}$
$\displaystyle=\frac{1}{J_{2}}\left[(J_{3}-J_{1})\Omega_{1}\Omega_{3}-\bar{J}_{r}\sum_{i=1}^{4}(-1)^{i+1}\Omega_{1}\omega_{i}+K_{r}d(\omega_{4}^{2}-\omega_{2}^{2})-\Omega_{2}|\Omega_{2}|C_{\tau_{2}}\right]$
(68) $\displaystyle\dot{\Omega}_{3}$
$\displaystyle=\frac{1}{J_{3}}\left[(J_{1}-J_{2})\Omega_{1}\Omega_{2}+K_{d}\sum_{i=1}^{4}(-1)^{i+1}\omega_{i}^{2}-\Omega_{3}|\Omega_{3}|C_{\tau_{3}}\right].$
(69)

###  3.4. Geometric formulation

Here we will equip $Q$ with a natural metric allowing us to characterize the
equations of motion as a second order system on $Q$. We then lift the dynamics
to a first order system on $TQ$.

Under our assumptions, the kinetic energy (43) can be encoded as the matrix

$\mathbb{G}=\begin{pmatrix}mI_{3}&0\\\ 0&J\end{pmatrix}$

where $m$ and $J$ denote the mass and inertia tensor of the drone,
respectively. This is a (diagonal) positive definite symmetric matrix and
induces a Riemannian metric on $Q$, which in turn yields the Levi-Civita
connection $\nabla$. This allows us to define the acceleration of a curve in
$Q$ as $\nabla_{\dot{\gamma}}\dot{\gamma}$, where a curve
$\gamma(t)=(b(t),R(t))$ has velocity
$\dot{\gamma}=(b(t),R(t),v(t),\Omega(t))$. Explicitly, we compute

$\nabla_{\dot{\gamma}}\dot{\gamma}=\begin{pmatrix}\dot{v}+\Omega\times v\\\
\dot{\Omega}+J^{-1}(\Omega\times J\Omega)\end{pmatrix}.$

This connection evidently contains the fictitious forces, and is independent
of choice of coordinates. The geodesics for $\nabla$ (curves with zero
acceleration) represent motions of a drone subject to no external forces.

Now in general, a Riemannian metric $\mathbb{G}$ on a manifold $Q$ induces the
“musical isomorphisms” $\mathbb{G}^{\sharp}:T^{*}Q\to TQ$ and
$\mathbb{G}^{\flat}:TQ\to T^{*}Q$. Newton’s equation can be expressed
geometrically as
$F(\gamma)=\mathbb{G}^{\flat}\nabla_{\dot{\gamma}}\dot{\gamma}$ where $F$ is
some force. In our setting, this produces the geometric version of our
equations of motion:

$\nabla_{\dot{\gamma}}\dot{\gamma}=\mathbb{G}^{\sharp}F(\gamma)$ (70)

where

$\mathbb{G}^{\sharp}=\begin{pmatrix}m^{-1}I_{3}&0\\\ 0&J^{-1}\end{pmatrix}$

and

$F(b,R)=\begin{pmatrix}f_{1}+f_{2}+f_{3}\\\
\tau_{f}+\tau_{1}+\tau_{2}\end{pmatrix}.$

Equation (70) represents a second order system on $Q$. We now express this as
a first order system on $TQ$. To that end, we invoke the vertical lift
$\text{vlft}_{w}:T_{q}Q\to T_{w}TQ$ where $w\in T_{q}Q$, and the geodesic
spray $S:TQ\to TTQ$ ([23, 16]). For our system, we can explicitly compute

$S(b,R,v,\Omega)=\begin{pmatrix}v\\\ \Omega\\\ v\times\Omega\\\
J^{-1}(J\Omega\times\Omega)\end{pmatrix}.$

Here the components are expressed relative to the standard left-invariant
basis of vector fields $\\{Y_{k}\\}$ on $TQ=TSE(3)$, where at a point
$(b,R)\in SE(3)$ we have $Y_{i}=(Re_{i},0)$ for $i=1,2,3$ and
$Y_{j+3}=(0,R\hat{e}_{j})$ for $j=1,2,3$.

Then our equations of motion can be expressed by the first order system

$\dot{\Gamma}=X(\Gamma)$ (71)

where $\Gamma$ is a curve in $TQ$ and $X$ is the vector field on $TQ$ given by

$X=S+\text{vlft}(\mathbb{G}^{\sharp}F).$ (72)

Next, we express this system as an affine nonlinear control system as in [3].
We can rewrite equation (71) as

$\dot{\Gamma}=\begin{pmatrix}v\\\ \Omega\\\
v\times\Omega+m^{-1}f_{1}+m^{-1}f_{2}\\\
J^{-1}(J\Omega\times\Omega)+J^{-1}\tau_{1}\end{pmatrix}+\begin{pmatrix}0\\\
0\\\ m^{-1}f_{3}\\\ J^{-1}\tau_{f}+J^{-1}\tau_{2}\end{pmatrix},$ (73)

where the first summand is the drift, and the second summand contains all the
controls $\omega_{i}$. Explicitly, we have

$\displaystyle m^{-1}f_{3}$
$\displaystyle=m^{-1}\sum_{i=1}^{4}K_{r}\omega^{2}_{i}e_{3}^{B}$
$\displaystyle J^{-1}\tau_{f}+J^{-1}\tau_{2}$
$\displaystyle=J^{-1}\begin{pmatrix}K_{r}d(\omega_{3}^{2}-\omega_{1}^{2})\\\
K_{r}d(\omega_{4}^{2}-\omega_{2}^{2})\\\
K_{d}\sum_{i=1}^{4}(-1)^{i+1}\omega_{i}^{2}\end{pmatrix}+J^{-1}\sum_{i=1}^{4}\Omega\times(-1)^{i+1}J_{r}\begin{pmatrix}0\\\
0\\\ \omega_{i}\end{pmatrix}.$

Note that the torque due to gyroscopic effects is linear in $\omega_{i}$,
while the force and torque due to thrust are linear in $\omega_{i}^{2}$.
Consequently we cannot write this system as an affine control system in
independent controls. As is common in the literature (e.g. [26, 32]) we now
simplify by ignoring the gyroscopic effects in $\tau_{2}$. This simplification
does not affect the simple motion planning in Section 4.1, but the gyroscopic
effects are present in the simulations of Section 4.2. Consider the controls

$u_{i}=\omega_{i}^{2}\quad\text{for}\quad i=1,2,3,4.$

######  Proposition 3.

Ignoring gyroscopic effects, the quadcopter can be expressed as an affine
nonlinear control system as:

$\dot{\Gamma}=f(\Gamma)+\sum_{i=1}^{4}g_{i}(\Gamma)u_{i}$ (74)

where the drift vector field is

$f(\Gamma)=\begin{pmatrix}v\\\ \Omega\\\
v\times\Omega+m^{-1}f_{1}+m^{-1}f_{2}\\\
J^{-1}(J\Omega\times\Omega)+J^{-1}\tau_{1}\end{pmatrix}$

and the control vector fields are

$\displaystyle g_{1}=\begin{pmatrix}0\\\ 0\\\ (0,0,m^{-1}K_{r})^{t}\\\
(-K_{r}d,0,K_{d})^{t}\end{pmatrix},\quad g_{2}=\begin{pmatrix}0\\\ 0\\\
(0,0,m^{-1}K_{r})^{t}\\\ (0,-K_{r}d,-K_{d})^{t}\end{pmatrix},$ $\displaystyle
g_{3}=\begin{pmatrix}0\\\ 0\\\ (0,0,m^{-1}K_{r})^{t}\\\
(K_{r}d,0,K_{d})^{t}\end{pmatrix},\quad g_{4}=\begin{pmatrix}0\\\ 0\\\
(0,0,m^{-1}K_{r})^{t}\\\ (0,K_{r}d,-K_{d})^{t}\end{pmatrix}.$

##  4\. Quadcopter multi-agent dynamics

###  4.1. Single-agent motion planning

We first explore some simple motions for a single quadcopter. We use the
coordinates and notation from Section 3.

######  Lemma 4.1.

Suppose a quadcopter undergoes yaw only, with no translational or other
rotational motion. Then

$\omega_{1}=\omega_{3}=\sqrt{\frac{mg}{2K_{r}}-\omega_{2}^{2}}\qquad\text{and}\qquad\omega_{4}=\omega_{2},$

where $\omega_{2}$ is a free parameter.

###### Proof.

The desired motion keeps all variables constant except $\Omega_{3}$ and
$\psi$. This imposes $\omega_{1}=\omega_{3}$ and $\omega_{2}=\omega_{4}$ as
well as

$\dot{b}=v_{1}=v_{2}=\theta=\phi=\Omega_{1}=\Omega_{2}=0$

for the duration of the motion.

The equations of motion in Lemma 3.1 simplify considerably; in particular, we
have $\dot{b}_{3}=v_{3}$ and $\dot{v}_{3}=\frac{1}{m}\sum_{i=1}^{4}t_{i}-g$.
Assuming $\dot{b}_{3}=0$ gives the expected result for hovering:
$\sum_{i=1}^{4}t_{i}=mg$. Using $\omega_{1}=\omega_{3}$ and
$\omega_{2}=\omega_{4}$ with $t_{i}=K_{r}\omega_{i}^{2}$ allows us to solve
for $\omega_{1}=\sqrt{\frac{mg}{2K_{r}}-\omega_{2}^{2}}.$

∎

######  Lemma 4.2.

Suppose a quadcopter moves only in the direction of the first body axis. Then

$\omega_{2}=\frac{1}{3}\left(\sqrt{\frac{6a-8K_{r}\omega_{4}^{2}}{K_{r}}}-\omega_{4}\right)\qquad\text{and}\qquad\omega_{1}=\omega_{3}=\frac{1}{2}(\omega_{2}+\omega_{4}),$

where

$a=m\left(\frac{1}{C_{\theta}}(\dot{v}_{1}S_{\theta}+v_{1}C_{\theta}\Omega_{2}+v_{3}S_{\theta}\Omega_{2})-v_{1}\Omega_{2}+gC_{\theta}\right)+v_{3}|v_{3}|C_{D_{3}}$

and $\omega_{4}$ is a free parameter.

###### Proof.

We desire to move only along the body $x$-axis with no motion in the body $y$
or $z$ direction, requiring the body to pitch without any yaw or roll. This
imposes $\omega_{1}=\omega_{3}$ as well as

$\dot{b}_{3}=\phi=V_{2}=\Omega_{3}=\Omega_{1}=0$

for the duration of the motion.

The equations of motion in Lemma 3.1 simplify considerably; in particular, we
have $\dot{\theta}=\Omega_{2}$ and
$\dot{b}_{3}=-V_{1}S_{\theta}+V_{3}C_{\theta}.$ Using $\omega_{1}=\omega_{3}$
and setting $\dot{\Omega}_{1}=0$ forces
$\omega_{1}=\omega_{3}=\frac{1}{2}(\omega_{2}+\omega_{4})$. Setting
$\ddot{b}_{3}=0$ and using $\dot{\theta}=\Omega_{2}$ gives

$\dot{v}_{3}=\frac{1}{C_{\theta}}(\dot{v}_{1}S_{\theta}+v_{1}C_{\theta}\Omega_{2}+v_{3}S_{\theta}\Omega_{2}).$

Setting this equal to the equation of motion for $\dot{v}_{3}$, we find

$\sum_{i=1}^{4}t_{i}=m\left(\frac{1}{C_{\theta}}(\dot{v}_{1}S_{\theta}+v_{1}C_{\theta}\Omega_{2}+v_{3}S_{\theta}\Omega_{2})-v_{1}\Omega_{2}+gC_{\theta}\right)+v_{3}|v_{3}|C_{D_{3}}.$

Denote the right side of this equation by $a$. Finally, using
$t_{i}=K_{r}\omega_{i}^{2}$, we can solve for

$\omega_{2}=\frac{1}{3}\left(\sqrt{\frac{6a-8K_{r}\omega_{4}^{2}}{K_{r}}}-\omega_{4}\right).$

∎

Note that in Lemma 4.2 we could alternatively treat $\omega_{4}$ as the free
parameter and solve for an identical expression for $\omega_{4}$ in terms of
$\omega_{2}$. We have similar choices for the yaw motion in Lemma 4.1.

Moreover, a nearly identical calculation gives path planning along the second
body axis. Treating $\omega_{3}$ as the free parameter, we can move only in
the body $y$ direction using

$\omega_{1}=\frac{1}{3}\left(\sqrt{\frac{6a-8K_{r}\omega_{3}^{2}}{K_{r}}}-\omega_{3}\right)\qquad\text{and}\qquad\omega_{2}=\frac{1}{2}(\omega_{1}+\omega_{3}).$

In order to move in another direction in the body $xy$-plane (with no motion
in the $z$ direction), we can simply combine the motions in Lemmas 4.1 and
4.2.

###  4.2. Simulations

####  4.2.1. Single agent

We demonstrate the simple motions from Section 4.1 in Figures 7 and 8. In the
coordinates of Section 3.3, the drone begins at $(0,0,0,0,0,0)\in SE(3)$ with
motors producing the necessary amount of thrust to remain at a hovering
equilibrium. Over the first 12 seconds, the drone maintains orientation while
rising vertically in the body (and inertial) $z$-direction, ending at another
hover. Over the next four seconds, the drone maintains position while rotating
$\pi/2$ radians counterclockwise about the body $z$-axis as described in Lemma
4.1. Finally, for the last four seconds the drone executes the motion in Lemma
4.2, flying directly along the body $x$-axis. Figure 7 shows the path of the
drone in space as well as the orientation angles over time. Figure 8 shows the
controls used to produce this motion.

Figure 7. A quadcopter rising straight up, then yawing while hovering, then
flying straight along the body $x$-axis. Top: position in space. Bottom:
orientation angles over time.

Figure 8. The controls used to produce the motion in Figure 7. Top: angular
velocities of each of the four motors over time. Bottom: total thrust over
time.

####  4.2.2. Multi-agent

Here we combine most of the work completed above in order to simulate an
agreement protocol for dynamically sound quadcopters. Consider three drones on
a rendezvous mission with given initial conditions. We assume each drone
begins hovering with all translational and angular velocities zero. This
forces each initial pitch and roll to also be zero, but we choose different
initial yaws. For simplicity, we assume all three drones start in the $z=0$
plane (which need not correspond to the ground) but choose different initial
$x,y$ positions. This simplifying assumption does not sacrifice much
generality: if the three drones started at three different nonzero altitudes,
we would simply find the average initial altitude and execute simple motions
bringing each drone up or down to this height.

In the coordinates $(b_{1},b_{2},b_{3},\phi,\theta,\psi)\in SE(3)$ from
Section 3.3, we choose our three drones to have initial configurations
$(0,0,0,0,0,0)$, $(0,9,0,0,0,-\pi/4)$, and $(15,9,0,0,0,\pi/2)$.

For this example we will assume a complete network $G=K_{3}$, so by Theorem
2.6 we have that all trajectories are straight lines. In particular, each
drone will traverse the line segment connecting their initial position in
$\mathbb{R}^{3}$ to the rendezvous position, which is simply the component-
wise average of the initial conditions: $(5,6,0).$ Figure 9 shows these
trajectories.

Figure 9. Three drones start at position $(0,0,0)$, $(0,9,0)$ and $(15,9,0)$
with initial yaw angles $0$, $-\pi/4$, and $\pi/2$. Their trajectories to the
rendezvous position are shown.

There are, however, two different methods for parametrizing these line
segments and hence planning the drones’ flights. First, we could use the
methods from Section 2, obtaining trajectories which exponentially approach
the rendezvous position. In this setting, the drones’ orientations are ignored
completely, as are all dynamics, giving a simple yet unrealistic path.

A main objective of this project was to investigate a second method: implement
the path-planning from Section 4.1 to find controls bringing the drones along
the desired line segments in a manner that respects the realistic dynamics
from Section 3. This requires each drone beginning its flight by executing a
pure yaw rotation to orient its body $x$-axis toward the rendezvous position,
then fly straight in this direction. See Figure 10 to see the yaw for each
drone over time; note that the total flight times differ for each drone. Also
note that for each drone, the roll over time is identically zero, but the
pitch is non-zero in order to produce lateral displacement (compare with
Figure 7).

Figure 10. The yaw $\psi$ over time for each drone. Total flight times differ
for each drone.

In order to create the desired trajectories in accordance with our equations
of motion, we utilize Lemmas 4.1 and 4.2. The resulting controls – the four
motor speeds for each of the three drones – are shown in Figure 11, as is the
total thrust for each drone. The parameters used in this simulation were taken
from [7], and appear in Table 1.

Table 1. Parameters used in simulations Constant | Symbol | Value
---|---|---
drone mass | $m$ | 0.468
drone inertia | $J$ | diag($(3.8278,3.8288,7.6566)\cdot 10^{-3}$)
rotor inertia | $J_{r}$ | diag$(0,0,2.8385\cdot 10^{-5})$
distance to rotor | $d$ | 0.25
thrust coefficient | $K_{r}$ | $2.9842\cdot 10^{-5}$
translational drag | $C_{D}$ | $(5.5670,5.5670,6.3540)\cdot 10^{-4}$
rotational drag | $C_{\tau}$ | $(5.5670,5.5670,6.3540)\cdot 10^{-4}$
propeller drag | $K_{d}$ | $3.2320\cdot 10^{-7}$

Figure 11. The controls used to produce the motions in Figures 9 and 10. We
show the four motor speeds and total thrust for each drone as functions of
time.

Finally, in Figure 12 we compare the paths generated by the two methods above.
The solid curves represent the simple but unrealistic paths derived in Section
2. The dashed curves represent the dynamically sound paths in Figures 10 and
11. Both methods trace the same curves shown in Figure 9, but differ in their
parametrizations. The top plot shows the $x$-coordinate over time; the solid
curves for drones 1 and 2 coincide as both start with an initial $x$-value of
0. The bottom plot shows the $y$-coordinate over time; the solid curves for
drones 2 and 3 coincide as both start with an initial $y$-value of 9. Note
that the solid curves approach (but do not reach) the rendezvous position much
faster than the dashed curves as they are unconstrained by orientation or
realistic acceleration. To move from the simulation environment to a more
realistic experiment with actual drones would require a more careful
consideration of the observer design as in [9].

Figure 12. Comparison of two methods of motion planning. Dashed curves include
dynamics; solid curves do not. Left: $x$-coordinate over time. Right:
$y$-coordinate over time.

##  5\. Conclusion

The primary goals of the paper were twofold. First, we introduce the agreement
protocol from a graph theoretic viewpoint and apply it to a rendezvous mission
for a set of particle agents without orientation. Second, we derive the
equations of motion for a quadcopter to take into account its orientation and
dynamics. We then apply the quadcopter dynamics to a rendezvous mission and
compare the trajectories with the ones from the network multi-agents
algorithm. This is a first step towards developing more complex multi-agent
systems for quadcopters. A next step will be to analyze the decoupling vector
fields for the quadcopter from its affine connection control system
formulation and use these trajectories to approximate geometric trajectories
obtained from the multi-agents network approach. We focused here on the
agreement protocol, but ongoing work shows that it can be extended to
formation flying for multi-agents and more complex missions, such as search
and rescue with a large team of quadcopters including leaders and followers.
We are also planning to run real drone flights to demonstrate the
applicability of our techniques.

## References

  * [1] W. Anderson Jr. and T. Morley,  Eigenvalues of the Laplacian of a graph,  _Linear and Multilinear Algebra_ , 18 (1985), 141–145.
  * [2] _AUVs on-board capability for image analysis for real time situation assessment_ ,  https://www.youtube.com/watch?v=qo0VWWvTgdU
  * [3] A. Bloch,  _Nonholonomic Mechanics and Control_ , $2^{nd}$ edition,  Springer-Verlag, New York, 2015.
  * [4] A. Bloch, I. Hussein, M. Leok, and A. Sanyal,  Geometric structure-preserving optimal control of the rigid body,  _Journal of Dynamical and Control Systems_ , 15 (2009), 307–330.
  * [5] A. Bloch, P. Krishnaprasad, J. Marsden, and G. Sanchez De Alvarez,  Stabilization of rigid body dynamics by internal and external torques,  _Automatica_ , 28 (1992), 745–756.
  * [6] A. Bloch, N. Leonard, and J. Marsden,  Controlled Lagrangians and the stabilization of mechanical systems I: The first matching theorem,  _IEEE Trans. Automat. Control_ , 45 (2000), 2253–2270.
  * [7] H. Bouadi, M. Bouchoucha, and M. Tadjine,  Sliding mode control based on backstepping approach for an UAV type-quadrotor,  _International Journal of Mechanical and Mechatronics Engineering_ , 1 (2007), 39–44.
  * [8] H. Bouadi and M. Tadjine,  Nonlinear observer design and sliding mode control of four rotors helicopter,  _International Journal of Aerospace and Mechanical Engineering_ , 1 (2007), 354–359.
  * [9] F. Bullo, J. Cortes, and S. Martinez,  _Distributed Control of Robotic Networks_ ,  Princeton University Press, 2009.
  * [10] F. Bullo and A. Lewis,  _Geometric Control of Mechanical Systems_ ,  Springer-Verlag, New York, 2004.
  * [11] Y. Cao, W. Yu, W. Ren and G. Chen,  An overview of recent progress in the study of distributed multi-agent coordination,  _IEEE Transactions on Industrial Informatics_ , 9 (2013), 427–438.
  * [12] R. Carney, M. Chyba et al.,  Multi-agents path planning for a swarm of unmanned aerial vehicles,  _2020 IEEE International Geoscience and Remote Sensing Symposium_ , (2020).
  * [13] R. Chandrasekaran, L. Colombo, M. Camarinha, R. Banavar, and A. Bloch,  Variational collision and obstacle avoidance of multi-agent systems on Riemannian manifolds,  _2020 European Control Conference_ , (2020), 1689–1694.
  * [14] S. Chung, A. Paranjape, P. Dames, S. Shen, and V. Kumar,  A survey on aerial swarm robotics,  _IEEE Transactions on Robotics_ , 34 (2018), 837–855.
  * [15] M. Chyba, R. Carney, C. Gray, and Z. Trimble,  Increasing small unmanned aerial system real-time autonomy,  _2020 IEEE International Geoscience and Remote Sensing Symposium_ , (2020).
  * [16] M. Chyba, T. Haberkorn, R. Smith, and G. Wilkens,  A geometric analysis of trajectory design for underwater vehicles,  _Discrete Contin. Dyn. Syst. Ser. B_ , 11 (2009), 233–262.
  * [17] F. Costa, J. Ueyama, T. Braun, G. Pessin, F. Osório and P. Vargas,  The use of unmanned aerial vehicles and wireless sensor network in agricultural applications,  _IEEE International Symposium on Geoscience and Remote Sensing_ , (2012).
  * [18] P. Doherty and P. Rudol,  A UAV search and rescue scenario with human body detection and geolocalization,  _Australian Conference on Artificial Intelligence_ , (2007).
  * [19] M. Erdelj and E. Natalizio,  UAV-assisted disaster management: applications and open issues,  _2016 International Conference on Computing, Networking and Communications (ICNC)_ , (2016).
  * [20] J. Kim, S. Kim, C. Ju and H. Son,  Unmanned aerial vehicles in agriculture: a review of perspective of platform, control, and applications,  _IEEE Access_ , 7 (2019).
  * [21] S. Knorn, Z. Chen and R. Middleton,  Collective control of multiagent systems,  _IEEE Transactions on Control of Network Systems_ , 3 (2016), 334–347.
  * [22] G. Lafferriere, A. Williams, J. Caughman, and J. Veerman,  Decentralized control of vehicle formations,  _Systems & Control Letters_, 54 (2005), 899–910.
  * [23] A. Lewis and R. Murray,  Configuration controllability of simple mechanical control systems,  _SIAM Journal on Control and Optimization_ , 35 (1997), 766–790.
  * [24] T. Luukkonen,  _Modelling and control of quadcopter_ ,  Independent research project, Aalto University in Espoo, Finland, 2011.
  * [25] M. Mesbahi and M. Egerstedt,  _Graph Theoretic Methods in Multiagent Networks_ ,  Princeton University Press, 2010.
  * [26] M. Mueller and T. D’Andrea,  Stability and control of a quadrocopter despite the complete loss of one, two, or three propellers,  _2014 IEEE international conference on robotics and automation (ICRA)_ , (2014), 45–52.
  * [27] I. Okoloko,  Path planning for multiple spacecraft using consensus with LMI avoidance constraints,  _IEEE Aerospace Conference_ , (2012), 1–8.
  * [28] R. Olfati-Saber, J. Fax, and R. Murray,  Consensus and cooperation in networked multi-agent systems,  _Proceedings of the IEEE_ , 95 (2007), 215–233.
  * [29] P. Rooney, A. Bloch, and C. Rangan,  Trees, forests, and stationary states of quantum Lindblad systems,  _arXiv:1810.11144_ , (2018).
  * [30] A. Sanyal, A. Bloch, and N. McClamroch,  Dynamics of multibody systems in planar motion in a central gravitational field,  _Dynamical Systems_ ,19 (2004),303–343.
  * [31] R. Smith, M. Chyba, G. Wilkens, and C. Catone  A geometrical approach to the motion planning problem for a submerged rigid body,  _International Journal of Control_ , 82 (2009), 1641–1656.
  * [32] V. Stepanyan and K. Krishnakumar,  Estimation, navigation and control of multi-rotor drones in an urban wind field,  _AIAA Information Systems-AIAA Infotech @ Aerospace_ , (2017).
  * [33] S. Waharte and N. Trigoni,  Supporting search and rescue operations with UAVs,  _2010 International Conference on Emerging Security Technologies_ , (2010).
  * [34] S. Wich and L. Koh,  _Conservation Drones: Mapping and Monitoring Biodiversity_ ,  Oxford University Press, New York, 2018.

Received xxxx 20xx; revised xxxx 20xx.


inkscapelatex=false

# Wireless Sensor Networks for Optimisation of Search and Rescue Management in
Floods

Harshil Bhatt  Department of Electronics
and Communication Engineering
Manipal Institute of Technology
Manipal, India
<EMAIL_ADDRESS>Pranesh G  Department of Mechatronics
Engineering
Manipal Institute of Technology
Manipal, India
<EMAIL_ADDRESS>Samarth Shankar  Shriyash Haralikar  Department of
Electronics
and Communication Engineering
Manipal Institute of Technology
Manipal, India
<EMAIL_ADDRESS>Department of Electronics
and Instrumentation Engineering
Manipal Institute of Technology
Manipal, India
<EMAIL_ADDRESS>
###### Abstract

We propose a novel search-and-rescue management method that relies on the
aerial deployment of Wireless Sensor Network (WSN) for locating victims after
floods. The sensor nodes will collect vital information such as heat
signatures for detecting human presence and location, the flow of flood. The
sensor modules are packed in a portable floating buoy with a user interface to
convey emergency messages to the base station. Sensor nodes are designed based
on disaster conditions, cost-effectiveness and deployed in the affected region
by a centrifugal dispersion system from a helicopter. A mobile ad-hoc network
is set up by modifying the Low Energy Adaptive Cluster Hierarchy (LEACH)
protocol for greater efficiency and adoption of multi-hop of Cluster Heads for
long-distance communication to Base Station. The model metrics have been
defined considering previous rural floods in India. The efficiency and power
characteristics of the network are are compared to other protocols via
simulations. The sensor data from the network makes resource management,
rescue planning and emergency priority more efficient, thus saving more lives
from floods.

###### Index Terms:

Wireless Sensor Network, Disaster Management, Search and Rescue operation,
Energy efficient networks

††publicationid: pubid:  978-1-6654-2849-1/21/$31.00 © 2021 IEEE

##  I Introduction

Floods are the worst natural disaster that claims more lives and casualties
than any other disasters in the world. Every year, an average of 1700 people
are killed, and 39 lakh people are affected by the ravages of the flood in
India [1]. The frequency of massive floods has increased due to climate change
and global warming. Recurrent floods have led to high causalities and huge
economic loss to the country. Due to floods, the estimated loss for the
country is close to 1.5 per cent of GDP in 1915-2015 [2]. Rural floods are
more devastating since most of rural population live close to rivers. The poor
infrastructure and remote location have been a critical challenge in disaster
management in rural areas. With the development of WSN, their application in
disaster management has received significant attention since they can collect
data more accurately with minimum human intervention. The WSN can help in
collection of on-ground information and effectively coordinate the rescue
operation. WSNs are resilient and remain online even during infrastructural
collapse due to disaster unlike their wired counterparts. In addition to
resilience, they are also are easy to scale and configure according to the
requirement. The previous research on WSN in disaster management has generally
focused on the ease of deployment [3] and utilizing UAV to increase the
surveillance [4]. To effectively harness the full potential of WSNs, one needs
to consider the main challenges faced by WSN for disaster management.

The primary challenge with WSN during a disaster are

Battery source and management: WSNs consume much more power to transceive data
than wired counterparts. The power consumption increases with the amount of
data it sends. Also, a vast number of nodes increases the complexity of power
management, decreasing the network’s efficiency.

Network Architecture: It plays a dominant role in the complexity and effective
routing of the data to the base station. The simple architecture is
challenging to scale and complex architecture are computationally expensive.
Also, the regular protocol does not account for the efficient power
consumption of the node. There is a trade-off between computational complexity
and efficiency in data transmission in the selection of network architecture.

Sensors: Since there is a limitation on the quantity of data transmitted,
sensors need to be carefully selected by trading off between the nodes’
resolution and power efficiency. Besides the restriction to the amount of
data, the minimum number of sensors needs to be determined to get the required
data to keep the cost low and increase deployment scale.

Sensor Node deployment: The deployment needs to ensure uniform distribution,
cover the complete area of disaster, ensure safe landing and retrieval of
sensor nodes. The mode of deployment is also crucial since this will directly
affect the time taken for deployment and survival rate of each node.

Complexity: The time taken for the relief to arrive once WSNs are deployed
crucial. Management of a large number of sensor nodes and data is a difficult
task.

Considering these challenges, We propose a solution based on deployment of WSN
during the time of flood to acquire real-time data from the affected region
that could help the response team act better. The primary contribution of this
paper are

  * •

Designed WSN and sensor node for flood relief application in rural areas with
minimum infrastructure or damaged infrastructure. The design ensures
scalability, cost-effectiveness, ease of deployment, and retrieval of the
sensor node.

  * •

Proposed a novel Modified LEACH (MLEACH) protocol for effective routing to
ensure minimum loss of data and speedy delivery of data packets to a base
station.

  * •

Proposed a layered method with the abstraction of data between each layer to
decrease complexity and improve data management in the network.

  * •

Proposed a unique method to filter out irrelevant and stationary data within
the cluster and thus preventing the data traffic at the base station.

Thus, we propose power and cost-effective WSN for transmitting vital on-ground
data like the victim’s location, the flow of flood, and the priority of
emergency will better optimize the rescue operation with limited response time
and resources.

Figure 1: The WSN deployment for collecting data of the victims stuck in flood
for optimizing the relief operation

Rest of the paper is organised as follows: Section II summarises the relevant
works followed by section III, discusses the deployment and node
specifications. Subsequently, Section IV summarises the Modified LEACH
(MLEACH) protocol. Section V puts forth the simulation results along with the
comparison of proposed Communication model against traditional model and is
concluded with Section VI and Section VII.

##  II Literature Survey

Our work builds on crucial previous work in WSN and remote sensing. Interest
in Early Warning Systems and Remote Sensing has picked up in the last few
years. Balis et al. used WSNs for real-time remote sensing for Early Warning
Systems and presented a framework to facilitate early warning systems [5].
Mojaddadi et al. using remote sensing data and Geographic Information System
proposed an ensemble method to create flood probability indices [6]. These
models are resource-intensive as they use satellite data which generates large
amounts of data. Supercomputers are used to manage and work on the data to
predict floods. These resources are not feasible for under-developed or
developing countries.

A significant challenge is the lifetime of the WSN, as it should last for the
entire duration of the search and rescue operation. Since energy is a scarce
resource that limits the lifetime of WSNs, several research works have
proposed energy-efficient routing protocols. Heinzelman et al. proposed Low
Energy Adaptive Clustering Hierarchy (LEACH). This clustering-based protocol
utilises the randomised rotation of local cluster heads to distribute energy
load among sensors nodes in the network evenly. The significant drawbacks of
this protocol were the single-hop communication between cluster heads and the
Base Station, and remaining energy in sensor nodes was not considered during
the election of cluster heads [7]. Sabet et al. proposed an energy-efficient
decentralised hierarchical cluster-based multi-hop routing algorithm for WSNs.
In their approach, clustering and multi-hop routing algorithms perform
simultaneously to reduce the energy consumption due to control packets [8].
Singh et al. compared different variations of the LEACH routing protocol based
on parameters like clustering, overhead, scalability, load balancing, and
complexity for WSNs [9]. Kodali et al. presented a multi-level hierarchical
protocol based on LEACH protocol which improved the energy efficiency and
lifetime of the network [10]. Since the cluster head (CH) in LEACH is random,
the probability of a node becoming a CH is equal for a node with low energy
and high energy. This affects the robustness and lifetime of the network.
Since LEACH follows single-hop communication between the CH and Base Station
(BS), the energy required by CH increases after a certain distance to transmit
data to the BS. This leads to uneven dissipation, where nodes far from the BS
die quickly; thus, reducing the effective sensing area. These issues are
addressed by implementing multi-hop communication between the CH and BS along
with beneficial attributes of existing LEACH protocol. Also, the proposed
model improves the random selection of CH by preventing the same node from
becoming a CH for the next $n$ rounds.

Survivor detection during disaster management is a popular topic amongst
researchers, and several sensors and algorithms have been proposed. Portmann
et al. proposed a tracking algorithm based on particle filter combined with
background subtraction for detecting and tracking people in thermal-infrared
images. [11]. Byunghun et al. analyse the performance and applicability of
passive infrared motion sensors in security systems and propose a region-based
human tracking algorithm [12]. While these algorithms substantially increase
human detection accuracy and detection speed, it is computationally expensive
to deploy on each sensor node. The advantage of using passive infrared sensors
is the reduction in processing power and cost.

Aerial deployment is one of the easiest, effective and robust way to deploy
WSN during a disaster. The Deployment method needs to achieve maximum coverage
in a minimum time. Sharma et al. proposed a Centrifugal Cannon based Sprinkler
random deployment algorithm to scatter sensor nodes from a helicopter moving
over the candidate region. This model is time efficient and scatters over vast
and unreachable areas while minimising the number of scans over the area [13].
In contrast, random dispersal mechanism may cause coverage holes right below
the vehicle, thus making targeted dropping nearly impossible.

Figure 2: CAD model of the Sensor Node

##  III Sensor Node Specification

The design of the proposed senor node is presented with utmost considerations
of both the safety and robustness aspect of the node. The overview of the
proposed solution is given by fig 1 The whole section consists of three parts,
Design and Deployment, Hardware Architecture as mention in 3 and Software
Architecture.

###  III-A Design and Deployment

After a Base Station is setup at a suitable location, the deployment of the
WSN is carried out by the helicopter with a centrifugal canon sprinkler
mechanism with added perpendicular drop system to scatter the nodes in all
direction and in significantly less time [13]. The pod proposed shown in fig 2
is made out of PVC foam with a density of $0.5g/cm^{3}$ to ensure floating and
cushions in landing during aerial deployment. The pod is designed to provide
stability in water by ensuring that the centre of gravity of the sensor node,
at no point is above the metacentre by placing heavy components close to the
base. The designed pod is 45 cm in height, and the maximum diameter of the
body is 30 cm. The whole pod with all the electronics weigh less than a
kilogram, making it a light and mobile sensor node so that victims can carry
it to a safe place in the locality to wait for the response team. GPS and the
activity status light would prove beneficial in retrieving sensor nodes after
use for multiple deployments.

Figure 3: Sensor Node Specifications

###  III-B Software Architecture

After the initial deployment of the WSN, only the sensors, siren, and WiFi
module are in an active state to save energy. Siren is programmed to create
loud noise at certain intervals to attract people. The operational status can
be inferred from the green LED, and the push button is placed along with it to
ensure the human presence with utmost certainty. An LCD screen and keypad are
installed on the sensor node to receive on-ground information from the victim
regarding the number of victims trapped, the severity of injury, age group,
and the presence of any medical emergency. The microcontroller has multiple
pipelines for polling sensor data, transceiving data across nodes and
processing data. These are ensured by multiple threads in the controller and
switching between the task in high frequency. A Real-Time Operating System
(RTOS) is used to ease scheduling and decrease latency. The RTOS also helps
set priority levels for different tasks, where human detection by the PIR
sensor can be prioritized more. A set of predetermined weights are associated
based on the level of emergency according to the data from the victims to
optimise the rescue operation with limited resources. With on-ground data, a
weighted graph optimisation problem can be formulated as a familiar Travelling
Salesman Problem(TSP) to find both effective and efficient routes for relief
operations.

##  IV Communication

###  IV-A Energy Radio model

The communication protocol has been designed to focus on the energy efficiency
of the networks. A dynamic energy model is considered, where the energy
dissipation depends on the number of packets being transmitted. Here, the
transmitter and receiver of the wireless sensor node spend 50nJ/bit energy. We
consider the transmitter of the amplifier requires
$\epsilon_{amp}=120PJ/bit/m^{2}$ to get the required range of Signal to Noise
Ratio (SNR).

$\displaystyle E_{Tx}(k,d)=E_{Tx-elec}(k)+E_{Tx-amp}(k,d)$ (1) $\displaystyle
E_{Tx}(k,d)=E_{elec}*k+\epsilon_{amp}*k*d^{2}$ (2) $\displaystyle
E_{Rx}(k)=E_{Rx-elec}(k)$ (3) $\displaystyle E_{Rx}(k)=E_{elec}*k$ (4)

where,
$E_{elec}$ = Minimum electronics energy required to run Tx or Rx circuitry,
$E_{Tx-elec}$ = Energy consumed due to transmit electronics,
$E_{Rx-elec}$ = Energy consumed due to receive electronics,
$\epsilon_{amp}$ = Transmit Amplifier energy consumption,
K – number of message bits,
d - distance in meters,
$E_{Tx}$ = Total Transmit energy consumed,
$E_{Rx}$ = Total Receive energy consumed,

Equation (2) and (4) give the total energy expended by transmitter and
receiver for n bits message and d meters distance. Thus, we get the total
energy consumed by wireless sensor nodes [10]. We use this radio energy model
in our MLEACH protocol.

###  IV-B Communication Protocol

The communication protocol proposed in this paper is an energy-efficient
Modified LEACH (MLEACH) protocol. The nodes in LEACH group themselves into
local clusters, with one node serving as cluster-head (CH).The usage of
clusters for data transmission minimizes the effective transmit distance
between transceivers by implementing hops among CH necessary to reach the base
station.

LEACH is a single hop routing protocol where each cluster head transmits data
directly to the base station. The limitation of a node to communicate only in
the vicinity of the base station and consumes significant amount of energy for
long distance communication making it unsuitable for disaster management
applications. To surmount this problem we propose MLEACH protocol which is a
multi-hop and layered protocol, making it apt for search and rescue operation.
In addition to this, an unique data filtering method is introduced to reduce
the data congestion and increase the energy efficiency of the networks.

The functioning of the protocol is separated into two phases: setup and data
transmission.

Set-up Phase: In this phase,The CH election, cluster formation and assignment
of the Time Division Multiple Access (TDMA) slots to each Cluster Member (CM)
is done. All nodes participate in the CH selection step by creating arbitrary
precedence value in 0-1 range. If a sensor node’s produced arbitrary number is
smaller than a threshold value $T(n)$, the node gets elected as CH. The value
of T(n) is calculated using equation (5) [7].

$T(n)=\begin{cases}\frac{P}{1-P*(rmod\frac{1}{P})}:ifn\in G\\\
0:else\end{cases}$ (5)

$P$ represents the percentage of node to emerge as CH among nodes, $r$ is the
ongoing round, and $G$ is collection of all nodes that were excluded in the CH
selection in $1/P$ rounds. This equation ensures that any node that hasn’t
been a CH has an equal probability of becoming one and energy gets distributed
equitably among the sensor nodes [9]. Once the CH is elected, it broadcasts
’hello’ messages to the nearby nodes and adds them to the cluster as CM. After
the cluster formation, the CH creates a TDMA and allocates a time slot to each
CM to transmit the sensor data to the CH. The TDMA slots prevents CM data from
colliding and allows their transmitter to be switched off until their slot is
available. This lowers each node’s energy usage and extends the network’s
lifespan.

Transmission Phase: The network is divided into the lower and upper layers for
transmission of data to the base station. Lower Layer: This layer consists of
all the CM’s. It deals with a conglomeration of sensor data from the CM on the
selected CH. Higher Layer: This layer consists of all the CH’s. Here, all the
CH’s processes the collected data from CM and filter the sensor data if a CM
has not changed beyond the threshold. Thus, only essential data is transmitted
to the base station via Dijkstra’s algorithm through other CH to achieve the
shortest route, ultimately reducing data traffic and improving energy
efficiency.

Algorithm 1 Proposed Routing Algorithm

1: begin round

2: if node!=CH then

3: while TDMA Slot exists do

4: if Node has data then

5: Transmit Data $\rightarrow$ CH

6: else

7: Send Heart Beat message $\rightarrow$ CH

8: end if

9: end while

10: else

11: Broadcast schedule $\rightarrow$ CM’s and receive data from them

12: if Received data $>$ Thresholds then

13: Send hello packet $\rightarrow$ neighbouring CH

14: Populate Routing Tables of each CH

15: initiate Dijkstra’s algorithm and send data $\rightarrow$ BS

16: End Round

17: else

18: Delete Data Packet

19: end if

20: end if

Figure 4: Routing Tree Diagram

In figure 4, we have illustrated the construction of the routing tree. Lower
layer data transmission: All the cluster members of cluster head S1, S2, S3,
S4, S5, S6, and S7 form the lower level. In the figure in the cluster radius
(Rc) of the S7 cluster head, there are 5 cluster members. They transmit data
to S7 according to their TDMA slot. Higher layer data transmission: All the
cluster heads S1, S2, S3, S4, S5, S6, and S7 form the higher layer. The base
station is in the radio range(Rr) of S1, S2, and S3 so, the data is
transmitted directly to the BS from them. For S4 CH there are two routes
$S4->S2->BS\;and\;S4->S1->BS$ . As $d(S4,S1,BS)=4.5>d(S4,S2,BS)=2.5$ , the
shortest route $S4->S2->BS$ is chosen for transmission of data. Similarly, the
routes are determined for S5, S6, and S7.

##  V Simulation

In this section, we evaluate the performance of our protocol and compare it
with another pro-active routing algorithm, Destination-Sequenced Distance
Vector (DSDV) via simulations. The simulation platform used is Network
Simulator 3 (NS-3). We assume that the BS is fixed and lies at a random
location within the sensor field. All the sensor nodes are randomly and non-
uniformly distributed in the sensor field. A random way-point mobility model
is used for simulating the movement of the sensor nodes. The IEEE 802.11ah
wireless networking protocol is used to provide internet connectivity to the
nodes. An On-Off Application was used to generate traffic to the BS in the
network. The simulation parameters are summarized in Table I.

Parameters | Values
---|---
Number of Nodes | 512
Sensor Field | 7500m*7500m
BS Location | Random (within sensor field)
Data Packet Size | 512 bytes
Initial Energy of nodes | 172800 J
Average Energy Consumption | 0.836 J
Simulation Time | 120s
TABLE I: Parameters of Simulation Figure 5: Total Energy vs Time of DSDV and
MLEACH Figure 6: Throughput comparison of DSDV and MLEACH

Figure 5 reveals a linear increase in the energy consumption as time passes.
Each node in DSDV consumes a maximum of 0.819J compared to 0.836J for MLEACH
to transmit data to the BS. Figure 6 shows the number of packets received by
the BS per second. The drop in throughput occurs when the number of available
CH reduces according to equation 5. While the MLEACH protocol consumes
slightly more energy than DSDV, the protocol’s throughput is vastly superior.
The major drawback of DSDV is the constant updating of its routing tables
resulting in more energy consumption. MLEACH can support an average of 400
packets per second in steady-state, and DSDV is only 200 packets per second
while consuming similar amounts of energy. DSDV would require double the
energy to match the throughput of MLEACH. This reduces the efficiency as
compared to MLEACH. Thus, MLEACH is more energy-efficient and can last longer
as compared to DSDV.

##  VI Conclusion

The flood is the most devastating natural disaster of all that ravages human
settlement throughout the globe every year. In India alone, close to two
thousand people lose their life due to the torrential floods that occur every
year. We have proposed a model to optimize the rescue operation with the help
of wireless sensor networks. The proposed Wireless Sensor Network uses the
custom MLEACH protocol that supports multi-hop routing and a novel cluster
head election algorithm to make the network more power-efficient. The WSN
nodes harbour essential sensors and transducers to help attract and locate the
location of victims, emergency priority and flow of the flood and transmits
relevant data to the base station to optimize the relief operation during
floods. The nodes are deployed in an aerial vehicle with a centrifugal
deployment system to ensure a wide coverage area. Sensor nodes are designed
based on landing impact, flotation stability, ease of deployment, and cost-
effectiveness to scale to the requirement. The proposed model is far better
than other communication models of a sensor network in terms of throughput.
The model has been tested in NS3 simulation, resulting in greater throughput
with the same power consumption compared to the DSDV protocol. This additional
throughput can establish communication with the victims making the rescue
process smooth and effective.

##  VII Acknowledgment

The authors are listed in ascending alphabetical order of names of the author
who have contributed equally to this paper. We thank our guide, Dr Ujjwal
Verma, for his constant support and guidance. Also, thanks to Mars Rover
Manipal for supporting and sponsoring the research.

## References

  * [1] Prakash Tripathi  “Flood Disaster in India : An Analysis of trend and Preparedness”  In _Interdisciplinary Journal of Contemporary Research_ 2.4, 2015, pp. 91–98
  * [2] Vikrant Panwar and Subir Sen  “Examining the economic impact of floods in selected Indian states”  In _Climate and Development_ 12.3 Taylor & Francis, 2020, pp. 281–296  DOI: 10.1080/17565529.2019.1614897
  * [3] Yasmin Jahir et al.  “Routing protocols and architecture for Disaster Area Network: A survey”  In _Ad Hoc Networks_ 82.August, 2019, pp. 1–14  DOI: 10.1016/j.adhoc.2018.08.005
  * [4] Mustafa Ilhan Akbaş and Damla Turgut  “APAWSAN: Actor positioning for aerial wireless sensor and actor networks”  In _Proceedings - Conference on Local Computer Networks, LCN_ , 2011, pp. 563–570  DOI: 10.1109/LCN.2011.6115518
  * [5] Bartosz Balis et al.  “A Survey on Successors of LEACH Protocol”  In _Procedia Computer Science_ 4, 2011, pp. 96–105  DOI: 10.1016/j.procs.2011.04.011
  * [6] Hossein Mojaddadi, Haleh Pradhan, Noordin Ahmad and Abdul Halim bin Ghazali  “Ensemble machine-learning-based geospatial approach for flood risk assessment using multi-sensor remote-sensing data and GIS”  In _Geomatics, Natural Hazards and Risk_ 8, 2017, pp. 1080–1102  DOI: 10.1080/19475705.2017.1294113
  * [7] Wendi Rabiner Heinzelman, Anantha handrakasan and Hari Balakrishnan  “Energy-efficient communication protocol for wireless microsensor networks”  In _Proceedings of the 33rd Hawaii International Conference on System Sciences - 2000_ , 2002  DOI: 10.1109/HICSS.2000.926982
  * [8] Maryam Sabet and Hamid Reza Naji  “A decentralized energy efficient hierarchical cluster-based routing algorithm for wireless sensor networks”  In _AEU - International Journal of Electronics and Communications_ 69, 2015, pp. 790–799  DOI: 10.1016/j.aeue.2015.01.002
  * [9] Sunil Kumar Singh, Prabhat Kumar and Jyoti Prakash Singh  “A Survey on Successors of LEACH Protocol”  In _IEEE Access_ 5, 2017, pp. 4298–4328  DOI: 10.1109/ACCESS.2017.2666082
  * [10] Ravi Kishore Kodali and Naveen Kumar Aravapalli  “Multi-level LEACH protocol using NS-3”  In _2014 International Conference on Circuits, Systems, Communication and Information Technology Applications, CSCITA 2014_ , 2014, pp. 132–137  DOI: 10.1109/CSCITA.2014.6839248
  * [11] Jan Portmann, Simon Lynen, Margarita Chli and Roland Siegwart  “People detection and tracking from aerial thermal views”  In _Proceedings - IEEE International Conference on Robotics and Automation_ , 2014, pp. 1794–1800  DOI: 10.1109/ICRA.2014.6907094
  * [12] Song Byunghun, Choi Haksoo and Su Lee Hyung  “Surveillance tracking system using passive infrared motion sensors in wireless sensor network”  In _2008 International Conference on Information Networking, ICOIN_ , 2008  DOI: 10.1109/ICOIN.2008.4472790
  * [13] Vikrant Sharma, R.. Patel, H.. Bhaduria and D. Prasad  “Policy for random aerial deployment in large scale Wireless Sensor Networks”  In _International Conference on Computing, Communication and Automation, ICCCA 2015_ , 2015, pp. 367–373  DOI: 10.1109/CCAA.2015.7148445