1

Decentralized Algorithms for 3D Symmetric
Formations in Robotic Networks:
a Contraction Theory Approach

arXiv:1511.02547v1 [cs.SY] 9 Nov 2015

Sumeet Singh, Student Member, IEEE, Edward Schmerling, Student Member, IEEE, Marco Pavone, Member, IEEE

Abstract—This paper presents decentralized algorithms for
formation control of multiple robots in three dimensions. Specifically, we leverage the mathematical properties of cyclic pursuit
along with results from contraction and partial contraction
theory to design decentralized control algorithms that ensure
global convergence to symmetric formations. We first consider
regular polygon formations as a base case, and then extend the
results to Johnson solid and other polygonal mesh formations.
The algorithms are further augmented to allow control over
formation size and avoid collisions with other robots in the
formation. The robustness properties of the algorithms are
assessed in the presence of bounded additive disturbances and
their effect on the quality of the formation is quantified. Finally,
we present a general methodology for embedding the control
laws on complex dynamical systems, in this case, quadcopters,
and validate this approach via simulations and experiments on
a fleet of quadcopters.

I. I NTRODUCTION
Recently there has been much interest in decentralized,
multi-robot systems due to the plethora of possible applications and performance advantages. These systems have many
attractive properties such as robustness to single-point failures,
scalable implementation, and potentially lower operating costs
when compared with monolithic approaches. Furthermore,
multiple robots can together accomplish tasks which may be
impossible or very difficult for any single robot on its own
[1]–[3].
A prototypical problem in this context is decentralized
formation control, that is, the problem of ensuring convergence
to a desired formation via control algorithms amenable to a
distributed implementation [4]–[6]. Most existing results consider two-dimensional formations [7]–[9], however, a number
of robotic applications would require three-dimensional (3D)
formations. For example, it may be desired to deploy a set
of small satellites to surround a larger damaged spacecraft
to provide complete 3D reconstruction and visualization of
the external structure. Additional examples include differential
atmospheric, deep space, or underwater measurements of environmental phenomena, and filming sport and entertainment
events [10]–[12].
Accordingly, the objective of this paper is to design and
rigorously analyze decentralized control algorithms for 3D
formations. Our approach leverages the simple, yet effective
strategy of cyclic pursuit. Essentially, the cyclic pursuit strategy entails letting each robot i follow its leading neighbor
Sumeet Singh and Marco Pavone are with the Department of Aeronautics
and Astronautics; Edward Schmerling is with the Institute for Computational
and Mathematical Engineering, Stanford University, Stanford, CA, 94305,
{ssingh19, schmrlng, pavone} @stanford.edu.
This work was supported in part by the Stanford Graduate Fellowship (SGF)
and the King Abdulaziz City for Science and Technology (KACST).

i + 1 modulo n, where n is the number of robots. This
approach is attractive due to its decentralized nature and low
information requirements, namely relative position, and has
led to its popularity in recent years [13]–[15]. In this paper,
we adopt the cyclic pursuit strategy proposed in [16] (that,
in turn, generalizes earlier results in [14], [15]), where each
robot follows its neighbors on both sides, that is, robots
i + m and i − m, where m = {1, 2, . . . N } is a lookahead parameter up to a horizon N . This is referred to as the
symmetric cyclic control algorithm. To prove stability (i.e.,
convergence to a desired formation), we exploit contraction
theory. Contraction theory is a relatively recent innovation
in control system design [17], and hinges upon an exact
differential analysis of convergence. At its core, contraction
theory studies exponential convergence of pairs of system
trajectories towards each other, and by extension, to a desired
target trajectory. This represents a generalization with respect
to traditional Lyapunov analysis which studies convergence to
the origin, i.e., to a zero trajectory [18]. Partial contraction
theory extends this concept by considering the convergence
of trajectories to a “set of properties,” for example, a flowinvariant subspace [19], [20].
The contributions of this paper are as follows. First, we
study convergence of the symmetric cyclic control algorithm
to a polygon formation in 3D. Our analysis differs from
the analysis in [16] in that it considers a refined notion
of the target formation (necessary to properly constrain a
polygon formation). Our analysis (as in [16]) leverages partial
contraction theory and, in addition, yields insights into the
various tuning parameters for the controller and their impact
on performance metrics such as convergence rate and control
saturation. Second, we extend the definition of the target
formation to accommodate combinations of regular polygons
such as the set of convex regular polyhedra known as Johnson
solids, as well as more general polygonal mesh formations.
Control laws and sufficient conditions for convergence are
derived for these relatively complex polyhedron formations,
paving the way for innovative mission concepts [21]. Third,
we provide extensions to the algorithm to control the size of
the formation, and prevent collisions with other robots. Fourth,
we analyze the robustness properties of the symmetric cyclic
controller under bounded additive disturbances and provide
performance bounds with respect to the formation error. These
results provide theoretical insights into experiments performed
on board the International Space Station which highlighted
the remarkable robustness properties of cyclic controllers [15].
Finally, we present a control hierarchy for implementing the
formation control laws on a complex dynamical system, in this
case, a fleet of quadcopters.
The rest of the paper proceeds as follows: Section II

2

introduces the mathematical background for contraction and
partial contraction theory and circulant matrices which forms
the basis for all subsequent analysis. Section III formalizes the
problem, introduces the formation subspace and symmetric
cyclic control law, and provides conditions for convergence
to a regular polygon. In Section IV, we extend the results
to polyhedral formations. Section V provides extensions to
the algorithm for controlling formation size and preventing
collisions among the robots. We also assess the robustness
properties of the symmetric cyclic controller under bounded
additive disturbances. In Section VI, we demonstrate the application of our formation control laws on quadcopters, providing
verification via simulations and hardware experiments. Finally,
in Section VII, we draw our conclusions and provide avenues
for future work.
A preliminary version of this paper appeared as [22]. This
extended and revised version contains as additional contributions: (1) extensions to the control algorithms for forming
more complex polygonal mesh formations, controlling formation size, and avoiding collisions, (2) additional numerical
results, and (3) hardware experimental results on a fleet of
quadcopters.
II. M ATHEMATIC P RELIMINARIES

for system (1). Third, the matrix M (x, t) := Θ(x, t)T Θ(x, t)
is referred to as the contraction metric.
We next discuss partial contraction theory, which allows one
to address questions more general than trajectory convergence.
Consider system (1), and assume there exists a flow-invariant
linear subspace M ⊂ Rn , i.e., a linear subspace with the
property that, for all t,
x ∈ M ⇒ f (x, t) ∈ M.
Assume that the dimension of M is p and let (e1 , . . . , en )
be an orthonormal basis of Rn where the first p vectors form
a basis for M. Let Ū be a p × n matrix whose rows are
eT1 , . . . , eTp . Let V̄ be an (n − p) × n matrix whose rows
are eTp+1 , . . . , eTn . One can easily verify that matrix V̄ is subunitary and satisfies the properties V̄ T V̄ +Ū T Ū = In , V̄ V̄ T =
In−p , and x ∈ M if and only if V̄ x = 0. We will refer to V̄
as the projection matrix of M. The main theorem in partial
contraction theory can be stated as:
Theorem II.2 (Partial Contraction [20]). Consider a flowinvariant linear subspace M and its associated projection
matrix V̄ . A particular solution xp (t) of system (1) converges
exponentially to M if the auxiliary system

ẏ = V̄ f V̄ T y + Ū T Ū xp (t), t

A. Notation
Given a square matrix A, the symmetric part of A, i.e.,
(1/2)(A + AT ), is denoted by Asym . The smallest and largest
eigenvalues of Asym are denoted, respectively, by λmin (A)
and λmax (A). Accordingly, the matrix A is positive definite
(denoted A  0) if λmin (A) > 0, and negative definite
(denoted A ≺ 0) if λmax (A) < 0. Let Ik denote the k-byk identity matrix. Let A ⊗ B denote the Kronecker product
between matrices A and B. The null space of a matrix A
is denoted by N (A). Finally, the set of eigenvalues of A is
denoted by λ(A).

is contracting with respect to y. If this is true for all particular
solutions xp , all trajectories of system (1) will exponentially
converge to M from all initial conditions.
Combining Theorems and II.1 and II.2, one obtains a
powerful tool to prove convergence to a desired flow-invariant
subspace.
Corollary II.3 (Convergence to Flow-invariant Subspace
[20]). A sufficient condition for global exponential convergence to M is
∂f T
V̄ ≺ 0,
uniformly.
∂x
In this paper, the subspace M will represent a desired robot
formation.
V̄

B. Contraction and Partial Contraction Theory
This section reviews basic results in nonlinear contraction
theory [17], [23]. The basic stability result in contraction
theory reads as follows.
Theorem II.1 (Contraction [17]). Consider a general system
of the form
ẋ = f (x, t),
(1)
where x is the n × 1 state vector and f is a n × 1 nonlinear, continuously differentiable vector function with Jacobian
∂f /∂x. If there exists a square matrix Θ(x, t) such that
Θ(x, t)T Θ(x, t) is uniformly positive definite and the matrix


∂f
Θ−1
F := Θ̇ + Θ
∂x
is uniformly negative definite, then all system trajectories
converge exponentially to a single trajectory. In this case, the
system is said to be contracting.
A few remarks are in order. First, a matrix Θ(x, t) is
uniformly positive definite if there exists β > 0 such that
∀(x, t), λmin (Θ(x, t)) ≥ β, and is uniformly negative definite
if there exists β > 0 such that ∀(x, t), λmax (Θ(x, t)) ≤ −β.
Second, the matrix F is referred to as the generalized Jacobian

Remark II.4 (Non-orthonormal V ). Note that the application of partial contraction theory requires the rows of the
projection matrix V̄ to be orthonormal. However, the matrix V characterizing a subspace M may not be row-wise
orthonormal, e.g., when it is obtained by combining a set of
linearly independent equations. Nonetheless, as long as the
equations are independent, the matrix V will be full row rank
and can be transformed via an invertible matrix T into an
orthonormal counterpart V̄ which satisfies (1) V = T V̄ , (2)
T
T
V̄ x = 0 ⇔ x ∈ M, and (3) V̄ ∂f
≺ 0 ⇔ V ∂f
≺0
∂x V̄
∂x V
[20].
C. Circulant and Block-Circulant Rotational Matrices
A circulant matrix of order n is a square matrix with the
following structure:


c1 c2 · · ·
cn
cn c1 · · · cn−1 
C=
(2)
..
.. 
 ..
.
.
.
.
c2 c3 · · ·
c1

3

The elements of each row are identical to the row above, but
shifted one position to the right and wrapped around. Thus, a
circulant matrix can be compactly denoted as:
C = circ[c1 c2 . . . cn ].
A useful circulant matrix with dimensions n × n used in
this paper is defined below, parametrized by the integer
m ∈ {0, . . . , n − 1}:
Lm := circ[1, 0, . . . , 0,

−1
|{z}

, 0, . . . , 0].

(m+1)st element

D. Kronecker Product
Let A and B be matrices of size m×n and p×q respectively.
The Kronecker product of A and B is defined as:


a11 B · · · a1n B

..  ,
A ⊗ B =  ...
. 
am1 B

···

amn B

where A ⊗ B has dimensions mp × nq. If λA is an eigenvalue
of A with eigenvector vA , and similarly λB and vB are an
eigenvalue and eigenvector pair for B then the corresponding
eigenvalue and eigenvector pair for A ⊗ B is λA λB and vA ⊗
vB . We also note the following identity: (A ⊗ B)(C ⊗ D) =
AC ⊗BD, where A, B, C and D have compatible dimensions
in order for the products AC and BD to be well defined.
In the following, we will denote by CL⊗R the set of blockcirculant matrices that can be written as L ⊗ Rβ , where L is
a circulant matrix and Rβ is a rotation matrix about the axis
ez := (0, 0, 1)T with rotation angle β. From the properties
of the Kronecker product, the eigenvalues and eigenvectors
of a matrix in CL⊗R are given, respectively, by the product
of eigenvalues and Kronecker product of eigenvectors of the
circulant matrix L and the rotation matrix Rβ .
III. S YMMETRIC P LANAR F ORMATIONS
A. Problem Setup
Consider n ≥ 3 mobile robots (uniquely labelled by an
integer i ∈ {1, . . . , n}). Denote the position of robot i at time t
as xi (t), where xi (t) ∈ R3 . The overall state vector is denoted
by x = (xT1 , xT2 , . . . , xTn )T . The dynamics of each robot are
given by
ẋi = gi (x) + ui (x, t),
(3)
where ui is the control action. It is desired to design the
control actions so that (a) they drive the global state vector x
to a desired symmetric formation, and (b) they are amenable
to a decentralized implementation. As in [16], our strategy is
to “encode” a symmetric formation as a formation subspace
as discussed next. The proofs for all theorems and lemmas
introduced in this section are provided in the appendix.
B. Formation Subspace
In this section we consider regular polygons in 3D as
desired symmetric formations – the extension to non-planar
formations (which represents the main contribution of this
paper) is discussed in Section IV. Consider the case where the
direction normal to the desired formation polygon is aligned
with the vector ez (the general case can be reduced to this case

via a coordinate transformation). We encode such a formation
via the subspace:
Mn = {x ∈ R3n :
(xi+1 −xi ) = R2π/n (xi+2 −xi+1 ), i = 1, . . . , n−2, (4a)
eTz (xn −xn−1 ) = eTz (x1 −xn )},

(4b)

where the indices are considered modulo n, and R2π/n denotes
a counterclockwise rotation around ez . The n − 2 constraints
in (4a) will be referred to as rotational constraints while the
single constraint in (4b) will be referred to as the in-plane
constraint. The in-plane constraint, not considered in [16], is
needed in order to ensure that all robots lie in the same plane
(which is not ensured by the rotational constraints alone as
they permit spiral formations). The following lemma shows
that the constraints (4a) and (4b) are indeed necessary and
sufficient for the definition of a regular polygon.
Lemma III.1 (Polygon Constraints). The set of constraints
(4a) and (4b) are necessary and sufficient for the definition of
a regular polygon with normal direction ez .
Remark III.2 (Polygon Degrees of Freedom). The constraints
(4a) and (4b) together form a set of 3n−5 linearly independent
equations in 3n variables. The five missing equations correspond to five distinct degrees of freedom: three in translation,
one in scaling, and one in in-plane rotation. That is, the
polygon may be translated anywhere in space, scaled in size,
or rotated within the desired plane about the plane normal.
Both (4a) and (4b) can be compactly represented as the null
space of a certain matrix, as shown in the following lemma.
Lemma III.3 (Compact Constraints). Let Wrn
:=
[In−2 , 0(n−2)×2 ] and Wpn := [01×(n−2) , 1, 0]. Define
the 3(n − 2) + 1 × 3n matrix


Wrn ⊗ I3
(L ⊗ I3 + (L1 − L2 ) ⊗ R2π/n ) . (5)
V :=
Wpn ⊗ eTz | 1
{z
}
{z
}
|
:=Pn ∈CL⊗R
:=Wn

Then the rotational and in-plane constraints in equation (4a)
and (4b) are equivalent to the equation V x = 0.
Note that Wrn captures the rotational constraints, while
Wpn captures the in-plane constraint. The subspace Mn can
then be characterized as the null space of matrix V , i.e.,
V x = 0 ⇔ x ∈ Mn .
Our standing assumption throughout this paper is that the
robots’ internal dynamics, i.e., g(x), are flow-invariant with
respect to Mn :
Assumption 1 (Flow-invariance). The internal dynamics are
flow-invariant with respect to the desired formation, that is for
all x ∈ Mn one has V g(x) = 0.
C. Symmetric Cyclic Controller
We consider the class of symmetric cyclic controllers proposed in [16] (in turn generalizing the pursuit controllers
introduced in [14]):
ui =

N
X



T
km Rm (xi+m −xi )+Rm
(xi−m −xi ) ,

(6)

m=1

where N is the look-ahead horizon (0 < N < n − 1), km > 0
is a gain, Rm is a rotation matrix around ez with rotation angle

4

αm , and (xi+m − xi ) and (xi−m − xi ) denote the relative
coordinates among robot i and its i + m and i − m neighbors,
where m = {1, . . . , N } and summation is modulo n. Note that
the control law is spatially-decentralized, as each robot only
requires relative position information from a set of neighboring
robots.
Then, the symmetric cyclic controller can be written in
compact form as
u=−
=−

N
X
m=1
N
X

T
km [Lm ⊗ Rm + LTm ⊗ Rm
]x

km Lm x = −L x,

(7)

m=1
T
where Lm := Lm ⊗ Rm + LTm ⊗ Rm
∈ CL⊗R and
N
P
L :=
km Lm . In order to apply partial contraction theory,
m=1

the control law needs to be flow-invariant (note that, by
Assumption 1, the internal dynamics are flow invariant). The
next lemma shows that this is indeed the case.
Lemma III.4 (Flow-invariance). Subspace Mn is flowinvariant with respect to the symmetric cyclic control law given
in equation (7).
We are now in a position to apply partial contraction theory
to show that, under some assumptions, the symmetric cyclic
controller drives the system to the desired formation subspace
Mn .
As the encoded constraints are linearly independent, V is
full row rank, guaranteeing the existence of an orthonormal
counterpart V̄ (see Remark II.4). In the following, to prove
convergence to Mn , we will apply partial contraction theory
using V̄ . Let Ū be a matrix whose rows represent an orthonormal basis for the orthogonal complement of the subspace
defined by the rows of V̄ . According to Theorem II.2, we
want to show that for system
ẋ = g(x) − L x,
the associated auxiliary system


ẏ =V̄ g(V̄ T y + Ū T Ū xp ) − L (V̄ T y + Ū T Ū xp ) ,
is contracting. Note that by Assumption 1 and Lemma III.4
the closed-loop dynamics are invariant with respect to Mn .
Then, according to Corollary II.3, one requires


∂g
− L V̄ T ≺ 0,
uniformly.
V̄
∂x
By Remark II.4, since V and V̄ are related by an invertible
transformation, the above stability requirement can be reformulated in terms of V directly. By factoring V as Wn Pn
(see equation (5)) one obtains the condition


∂g
Wn Pn
− L PnT WnT ≺ 0,
uniformly.
(8)
∂x
Performing an eigenvalue analysis of (8) yields the main result
of this section. The details of the proof are provided in the
appendix.

Theorem III.5 (Polygon Convergence). Assume




N
X
∂g T
(m)
Pn − min
sup λmax Pn
km λik < 0,
1≤i≤n
∂x
x,t
m=1

(9)

k∈{−1,0,1}

where
(m)
λik

2

=
e

2π(2(i−1)+k)j
n





2πm(i − 1)
cos(kαm ) − cos kαm +
n


2
2π(i−1+k)j
2π(i−1)j
n
.
e
−1 e n
−1

Then the robots governed by system (3) under the symmetric
cyclic controller (6) globally converge to a regular polygon
formation, i.e., to the formation subspace Mn .
Note that due to the necessary inclusion of the in-plane
constraint, the eigenvalues in Theorem III.5 differ from those
in [16].
Remark III.6 (Formations of Fixed Size). Theorem III.5
provides a sufficient condition to ensure that the robots converge to a regular polygon formation. Note that depending on
the internal dynamics g(x) and the rotation angle αm , the
polygon formation will be contracting, expanding, or of fixed
size. For instance, if g(x) = 0, by setting αm < mπ/n the
polygon formation will contract towards a point; by setting
αm > mπ/n the polygon formation will expand as the robots
travel along rays emanating from the center of mass of the
formation; by setting αm = mπ/n the polygon formation
will be of fixed size, where the size depends on the initial
conditions. The freedom in scaling the formation will be
utilized later to control the formation size.
For a given number of robots, the primary design parameters
in the symmetric cyclic controller include the gains km and the
look-ahead horizon N . It is clear that increasing the gains uniformly scales all eigenvalues of the projected Jacobian, thereby
increasing the exponential convergence rate. However, this is
at the expense of a more aggressive controller which increases
the risk of control saturation. On the other hand, increasing
the look-ahead horizon, while imposing greater information
requirements for each robot, also increases the convergence
rate but at a lower risk of control saturation. This is due
to the fact that the symmetric cyclic controller leverages the
degree of asymmetry at each node of the polygon. Increasing
the look-ahead horizon therefore improves the “estimate” of
this asymmetry, rather than simply magnifying the control law
by increasing the gain. This intuition is exemplified by the
following simulation.
For simplicity, we assume that the internal dynamics g(x)
are zero. In order to compare the control signals due to
varying choices in gains and look-ahead, we maintain the
same lower bound on the closed-loop convergence rate, i.e.,
λmin V̄ LV̄ . Figure 1 shows the Euclidean norm of the net
control signal, i.e., kLx(t)k2 for six robots starting at the same
initial conditions. The control signal in Figure 1a demonstrates
a significantly lower peak value than the signal in Figure 1b
despite the fact that the two closed-loop systems possess the
same lower on bound on convergence rate.
Remark III.7 (General Plane of Convergence). For polygons
not orthogonal to ez , one can use a similarity transformed
version of the control rotation matrix Rm to obtain analogous
convergence results. Let RηT ez be the desired polygon normal

5

30

10

kuk2

kuk2

20
5

10
0

0.2

0.4

0.6 0.8
Time [s]

1

1.2

0
0

0.2

0.4
0.6
Time [s]

0.8

1

(a) Control norm with N = 2, (b) Control norm with N = 1,
km = 2.
km = 6.928.

Fig. 1: Comparison of control effort with varying gains and
look-ahead horizon.
where Rη is an arbitrary rotation matrix. Define Rms :=
RηT Rm Rη to be the similarity transform of Rm . We now
replace Rm in the individual control law given in (6) by Rms ,
and the expression in (7) becomes:
N
X
u = − (In ⊗ RηT )
km Lm (In ⊗ Rη ) x
| {z }
| {z } m=1
:=RT
η

=−

N
X

:=Rη

(10)

km Lmη x = −Lη x,

m=1

where Lmη := RTη Lm Rη and Lη :=

N
P

km Lmη . The

m=1

expression for the projection matrix V becomes:
V = Wn Pn Rη .

(11)

The convergence analysis under control law (10) is virtually
identical to the one provided above for the case where the
polygon is orthogonal to ez . The generalized control law (10)
will be leveraged in Section IV to design control laws for
polyhedral formations.
IV. P OLYHEDRAL F ORMATIONS
In this section we extend the results of Section III to
polyhedral formations. Specifically, we first focus on convex
polyhedra having regular faces and equal edge lengths, referred to as Johnson solids. Generalizations to other polygonal
mesh formations are discussed at the end of this section.
Definition IV.1 (Convex Polyhedral Solid and Polyhedral
Surface). A convex solid polyhedron, or simply polyhedron,
is defined as a body bounded by a finite number of polygons
such that it lies on one side of the plane of each polygon.
Equivalently, it is the convex hull of the polygon vertices. A
complete polyhedral surface is the set of bounding faces for
a polyhedron. A partial polyhedral surface (PPS) is a subset
of the faces of a complete polyhedral surface.
Figure 2 illustrates the distinction between a complete and
partial surface of a convex solid in the definition above. In
Figure 2a, only a portion of the surface of the Johnson solid
(cube) is specified, and thus the PPS formed by the faces
in Figure 2a contains un-matched edges. The set of these
unmatched edges (shown in dashed form in Figure 2a) is
termed the boundary of the PPS. The complete polyhedral
surface displayed in Figure 2b encompasses the entire surface
of the solid cube.
Consider now a polyhedron P with complete polyhedral surface F. As an extension to the discussion in Section III, each

(a) Partial polyhedral surface of a (b) Complete polyhedral surface
cube.
of a cube.

Fig. 2: Partial and complete representations of a bounded
polyhedron (cube) with same vertex set.
face of F can be described by a set of linearly independent
constraints on its vertices. When considering these faces in
combination, we must ensure that the full set of constraints
(i.e., due to all faces in F) remain linearly independent.
This motivates the search for a PPS representation of the
polyhedron that is minimal in the sense that it does not induce
redundant constraints yet still retains enough information to
fix the positions of each vertex (robot). We note that this
reduced set of linear constraints for the PPS will retain the
translational and scaling degrees of freedom associated with
the PPS’s constituent polygons (see Remark III.2). Thus all
further discussion of the terms “face” or “polygon” in this
section will implicitly assume translation and scale invariance.
In the following subsection, we show that a suitable minimal
PPS representation can always be constructed for any given
convex polyhedron.
A. Formulation of a Minimal Partial Polyhedral Surface
We first introduce some terminology (adapted from [24]):
Definition IV.2 (Development of a Polyhedral Surface [24]).
A development D is defined as the tuple (Q, R) where Q is
a set of polygons as defined in Section III and R is a set of
rules that identify common edges between polygons. The sets
Q and R satisfy the following properties:
1) For every pair of adjacent polygons Qi , Qj ∈ Q, there
is a rule rij ∈ R which encompasses two pieces of
information, namely, it identifies (a) the two shared
vertices between the two polygons, and (b) the dihedral
angle formed between the planes of the two polygons.
2) It is possible to create a path between any two polygons
in the collection by following the common edges (i.e., a
development cannot be split into disconnected parts).
3) Each side of every polygon is either matched with no
edge of any other polygon, or with exactly one edge of
another polygon in the collection.
Every development defines a unique (complete or partial)
polyhedral surface F, but multiple developments may correspond to the same surface. The surface in Figure 2a corresponds to a development with three squares, while the surface
in Figure 2b corresponds to a development with six squares
(this development is commonly referred to as the polyhedral
net [25]). Our goal is to be able to use the development for the
PPS in Figure 2a as a viable representation of the complete
polyhedral surface in Figure 2b by leveraging the following
result:

6

Lemma IV.3 (Completion of a Polyhedral Surface [24]).
Every bounded convex PPS F 0 gives rise to a unique complete
polyhedral surface F without adding new vertices. If A0
denotes the set of vertices for F 0 , then the completion F is
given by the surface of the convex hull of A0 .
We are now ready to present the first step in constructing
a PPS representation of a complete polyhedral surface where
the set of linear constraints induced by the faces of the PPS
are linearly independent.
Lemma IV.4 (Minimal Development). Let F be a complete
polyhedral surface and let A denote its set of vertices. Suppose
D = (Q, R) is a development for F. Choose a set of polygons
Q0 ⊂ Q that satisfies the following properties:
1) Tree: Within the dual graph of F, the subgraph induced
by Q0 is a tree,
2) Vertex span: The set of vertices spanned by Q0 is identical
to the set of vertices spanned by Q.

Fig. 3: Indexing of robots for a general polyhedron.
For each Qk ∈ Q, define the matrix Ē (k) ∈ R|Vk |×n as
(
1 if jk = i,
(k)
[Ē ]ij =
(12)
0 otherwise.

R0 = {rij ∈ R | Qi , Qj ∈ Q0 } .

Thus [Ē (k) ]ij equals one if robot j’s index within face k is
equal to i. Define x(k) ∈ R3|Vk | to be the positions of the
robots in face Qk . Using (12), we may write x(k) as


x(k) = Ē (k) ⊗ I3 x = E (k) x,
(13)

Then, the tuple D0 := (Q0 , R0 ) is a well-posed development
(essentially a spanning tree of the graph D) and therefore
defines a PPS F 0 . Furthermore, the completion of F 0 is the
original surface F.

where x ∈ R3n is the global state vector, and E (k) :=
Ē (k) ⊗ I3 . As a direct extension of Section III-B, the face Qk
(by assumption a regular polygon) is defined as the subspace
(k)
M|Vk | as follows:

Given Q0 satisfying the properties above, we now construct the
set of rules R0 ⊂ R corresponding to the polygons selected
in Q0 :

Proof. Properties 1) and 3) of Definition IV.2 hold for Q0 , R0
by construction. Property 2) holds because the ‘tree’ assumption for Q0 ensures that there exists a path between any two
polygons in Q0 by following common edges. Thus D0 is welldefined and therefore describes a PPS F 0 .
The ‘vertex span’ assumption ensures that the convex hulls
of the vertices of F and F 0 are identical. Since F, a complete
polyhedral surface, is the surface of this convex hull, by
Lemma IV.3 the completion of F 0 is F.
The key message of Lemma IV.4 is that a PPS is sufficient
to describe a polyhedron (i.e., the desired formation) provided
that it encompasses all vertices of the polyhedron and its
faces are connected. More precisely, each face within the
PPS defines a set of linear constraints on its vertices, the
combination of which completely determine the convex hull
matching the desired formation. The further restriction that
the faces induce a tree graph is an important step, as shown
next, towards ensuring that the combined linear constraints are
independent. For the rest of this section, we assume that a PPS
F described by the development D = (Q, R) has been fixed.
Denote L := |Q|, the number of polygons in F. We now make
explicit the notion of rules of a development.
Let Qk denote the k th polygon in Q. Each Qk is associated
with a set of robots Vk , and a unit outward normal nk (pointing
outward from the interior of the polyhedron). The number of
robots n is equal to the number of vertices in the polyhedron
and each robot is uniquely associated with a vertex. Each face
induces a counterclockwise (with respect to nk ) sub-indexing
of the robots, specifically, let ik denote the index of robot i in
face k, where i = 1, . . . , n and k = 1, . . . , L, see Figure 3. If
robot i does not belong to face k, then we use the convention
ik = ∞.

(k)

M|Vk | = {x(k) ∈ R3|Vk | :
(k)

(k)

(k)

(k)

Rηk (xi+1 −xi ) = R2π/n Rηk (xi+2 −xi+1 ),
i = 1, . . . , |Vk |−2,
(k)

(k)

(k)

(k)

eTz Rηk (x|Vk | −x|Vk |−1 ) = eTz Rηk (x1 −x|Vk | )},

(14a)
(14b)

where Rηk is a rotation matrix defined such that RηTk ez = nk .
The desired Johnson polyhedral formation is the convex hull
of a vertex set which can be represented as a subspace Mn ⊂
R3n defined by the polygonal constraints that stem from each
face in Q, that is
\
(k)
Mn =
M|Vk | .
(15)
k:Qk ∈Q

By Lemma III.3, the constraints in (14a) and (14b) can be
compactly represented as the null space of a 3(|Vk | − 2) + 1 ×
(k)
3n matrix, denoted by Ve (k) . In other words, x ∈ M|Vk | ⇔
Ve (k) x = 0. Using (11) and (13), one can write Ve (k) as
Ve (k) = W|Vk | P|Vk | Rηk E (k) .

(16)

Thus, the subspace Mn in (15) may be equivalently re-stated
as the intersection of null spaces:
\
Mn =
N (Ve (k) ).
k:Qk ∈Q

Note that the rules R are implicitly encoded in the indexing
used to write down the constraint matrices Ve (k) and their
nullspace intersection. Given our representation of the PPS
in terms of the set Q, we can now combine (15) and (16) to
form the global constraint matrix Ve :
h
iT
T
T
T
Ve := Ve (1)
.
(17)
Ve (2)
· · · Ve (L)

7

The global formation subspace Mn is then defined to be the
null space of Ve , i.e.,
x ∈ Mn ⇔ Ve x = 0.
As constructed, Ve is not yet full row rank. A geometric explanation for this linear dependence is that once the orientations
of two adjacent PPS faces are fixed, the in-plane rotational
degrees of freedom for the remaining faces vanish. This notion
is formalized in the following (relatively straightforward)
corollary of Lemma III.1, the proof of which is provided in
the appendix.
Corollary IV.5 (Reduced Polygon Constraints). Assume two
neighboring robots j and j + 1 of face Qk are constrained to
lie in a plane normal to nk , where the indices {j, j + 1} ∈
{1, . . . , |Vk |} are modulo |Vk |. Then the rotational constraints
given in (14a) are necessary and sufficient for the definition
of a regular polygon with normal direction nk .
Since stability analysis using partial contraction is contingent upon a full row rank constraint matrix, we present a reduction of Ve that discards redundant constraints. Consider the
following reduced constraint matrix V with similar structural
properties as Ve :
h
iT
T
T
T
V := V (1)
,
(18)
V (2)
· · · V (L)
(
Ve (k)
if k = 1, 2,
V (k) =
(19)
(Wr|Vk | ⊗ I3 )P|Vk | Rηk E (k) k = 3, . . . , L,
where, as defined in Lemma III.3, Wr|Vk |
=
[I|Vk |−2 , 0(|Vk |−2)×2 ]. That is, in-plane constraints have
been removed from all but two (adjacent) faces within Q,
specifically, Q1 and Q2 . The following lemma shows that the
reduced constraints that stem from (19) are equivalent to the
original set of constraints using Ve (k) , and that the resulting
global constraint matrix V is full row rank.
Lemma IV.6 (Minimal Representation of V ). Denote the set
of equations Ve (k) x = 0 where Ve (k) has the form given in (16)
for k = 1, . . . , L as the full-constraint set. Similarly, denote
the set of equations V (k) x = 0 where V (k) has the form given
in (19) as the reduced-constraint set. Then, the solutions to the
two sets of equations are identical. That is,
x ∈ Mn ⇔ V x = 0.
Furthermore, V is full row rank.
Proof. We assume without loss of generality that the faces
Qk (and thus the corresponding constraint matrices V (k) )
are ordered so that for any j = 1, ..., L all faces Qi
along the unique path between Q1 and Qj satisfy i ≤ j.
Let Qk = {Q1 , . . . , Qk } denote a partial list of faces,
constrained by the set of equations [V (k) ]x = 0, where
h
iT
T
T
, k ≤ L. Similarly, let
[V (k) ] = V (1)
· · · V (k)
h
iT
T
T
[Ve (k) ] = Ve (1)
. We proceed by induction.
· · · Ve (k)

Base case k = 2: It is clear that N [V (2) ] = N ([Ve (2) ]) as
[V (2) ] = [Ve (2) ]. To see that [V (2) ] is full row rank, we note
that V (1) is full row rank (Lemma III.1) and each successive
rotational constraint within V (2) (counting around the polygon
starting from the shared edge) involves at least one robot that

does not belong to any previous constraints in V (1) or V (2) .
The in-plane constraint for V (2) further reduces the dimension
of the solution set by constraining the rotational degree of
freedom for Q1 : there is only one orientation of Q1 (up to
reflection, i.e. negative scaling) under which the edge shared
by the two faces is orthogonal to n2 . Therefore all constraints
in [V (2) ] are independent.
Inductive Hypothesis: Assume that the set of equations
[Ve (k) ]x = 0 and [V (k) ]x = 0 are equivalent and that [V (k) ]
is full row rank.
We now prove that [Ve (k+1) ]x = 0 and [V (k+1) ]x = 0 are
also equivalent and that [V (k+1) ] is full row rank. The inductive hypothesis [V (k) ]x = 0 encodes rules which assemble the
faces of Qk (which form a tree by the ordering assumption)
into a unique PPS. Now, from our ordering assumption, Qk+1
must share exactly one edge (i.e., exactly 2 robots) with a face
within Qk , and that the edge must lie in a plane normal to
nk+1 . The tree structure of Qk prevents any additional shared
edges. From Corollary IV.5, the rotational constraints encoded
in V (k+1) are indeed sufficient to fully constrain Qk+1 . Thus,
[V (k+1) ]x = 0 and [Ve (k+1) ]x = 0 are equivalent.
To see that [V (k+1) ] is full row rank, we note that each
successive rotational constraint within V (k+1) involves at least
one robot not represented by any of the existing constraints
within [V (k) ] or V (k+1) . Then, by similar reasoning as for the
base case, all constraints in [V (k+1) ] are independent.
By induction, both claims are proven.
Remark IV.7 (Polyhedron
Degrees of Freedom). Note that
 L
P
V has a total of 3
|Vk | − 6L + 2 linearly independent
k=1
 L

P
equations for a total of 3
|Vk | − 6L + 6 variables.
k=1

Three degrees of freedom correspond to the three translational
degrees of freedom. The final degree of freedom corresponds
to scaling the polyhedron by some constant α ∈ R. Applying
a negative scaling factor may be interpreted as inverting the
entire formation.

B. Polyhedral Formation Control Law
Having characterized the desired formation as the null space
of a full row rank matrix, we turn our attention to the control
laws. The control law for each robot is given by the sum
of a number of contributions, one for each face in which
(k)
the robot represents a vertex. Specifically, let ui denote
the contribution to the control law for robot i due to face
k, obtained by generalizing the symmetric cyclic controller
presented in Section III-C, namely:
(k)

ui =

Nk
X
m=1

 (k) (k)
(k)
(k)
(k) 
(k)
(k)T
km
Rms (xi+m − xi ) + Rm
(xi−m − xi ) ,
s

(20)
where Nk < |Vk | − 1 is the look-ahead horizon for face
(k)
(k)
(k)
(k)
k, km > 0 is a gain, Rms = RηTk Rm Rηk with Rm =
(k)
(k)
Rmπ/|Vk | , and xi+m and xi−m are neighboring robots within
(k)
face k (modulo |Vk |). The choice Rm = Rmπ/|Vk | stems from
the requirement that robots converge to a polygon of fixed size
within each face (see Remark III.6).

8

The net control for each robot is then given by the superposition of the contributions from each face in which the robot
is a vertex, i.e.,
ui =

L
X

(k)

ui ,

i = 1, . . . , n.

(21)

k:ik 6=∞

The control laws in (20) and (21) highlight the information
requirements for each robot in the formation. In particular,
we assume that each robot knows: (a) what face(s) it belongs
to, (b) relative position with respect to its forward and rear
neighbors within each face, (c) the number of robots in these
faces, and (d) the rotation matrix Rη that accompanies the
relative position measurements. In a scenario where there
are a large number of robots in the overall formation, these
information requirements are fairly minimal, emphasizing the
decentralized framework that underpins this work.
Using (20) and (10), the overall control vector stemming
from face k is given by
u(k) = −

Nk
X

(9), required for convergence to a single regular polygon. In
particular, we now have the addition of cross-terms in the
projected Jacobian due to the interactions between the various
faces. Additionally, when working with the auxiliary system,
we did not leverage the results of flow-invariant subspaces
as they were introduced in Section III. This is because if the
global state vector x converges to one of the subspaces defined
by V (k) , then the global dynamics will not necessarily be flowinvariant with respect to the given face or any other faces in
Mn . This is due to two reasons: (a) only the rotational set
of constraints are used to define faces other than the ones
indexed by k = 1 and k = 2, which we know are insufficient
on their own to describe a regular planar polygon, and (b) the
symmetric cyclic control law for each robot as defined in (21)
introduces coupling between various faces. Thus, convergence
is predicated on the decay of the symmetric cyclic controller
to zero for all faces.
Figure 4 shows simulation results demonstrating convergence to example 3D formations.

(k) (k) (k)
(k)
km
Lmη x := −L(k)
x,
η E

m=1

where
T

(k)
Lmη

(k)

(k)

(k)

= RTηk Lm Rηk and Lm = Lm ⊗ Rm + LTm ⊗

(k)

Rm . The overall closed-loop dynamics are then given by
ẋ = u =

L
X

T

(k)
−E (k) L(k)
x.
η E

(22)

k=1

For simplicity, we are considering zero internal dynamics. A
sufficient condition for convergence to the desired formation,
i.e., to the subspace Mn , is provided by the following theorem.
Theorem IV.8 (Polyhedron Convergence). Assume
  L


P
T
(k)
J= V
−E (k) Lη E (k) V T ≺ 0.

(23)

Then, the closed-loop dynamics (22) globally converge to a
Johnson polyhedral formation, i.e., to the formation subspace
Mn .
Proof. Let V̄ represent the orthonormal counterpart of V ,
whose existence is guaranteed given the results of Lemma
IV.6. First, we need to show that the dynamics in (22) are flow
invariant. Indeed, as an immediate consequence of Remark
III.6 and the fact that the rotation angle is set equal to
mπ/|Vk |, one has
for all k = 1, . . . , L,

(b) Octahedron (6 robots)

Fig. 4: Convergence to polyhedra. Faces with bold edges
correspond to the set Q.

k=1

x ∈ Mn ⇒ u(k) (x) = 0,

(a) Hexagonal Box (12 robots)

(24)

and hence x ∈ Mn ⇒ V̄ u = 0, i.e., the dynamics in (22) are
flow invariant.
Consider, then, the following auxiliary system for (22):
!
L
X

(k)T (k) (k)
T
T
ẏ = V̄
−E
Lη E
V̄ y + Ū Ū xp .
k=1

Since the dynamics are flow invariant, and an invertible
transformation exists between V and V̄ , then by Remark II.4
and by applying Corollary II.3, one obtains the claim.
Remark IV.9 (Conditions for Polyhedron Convergence). Note
that the condition given in (23) is markedly more complex than

Remark IV.10 (Nonsymmetric Tree-Based Formations).
Lemma IV.4 sets forth the essential properties of polyhedral
surface developments required for the preceding control law
and convergence analysis to apply. That is, viable targets under this theory consist of any formation of robots expressable
as the vertices of spatially-oriented regular polygons arranged
according to a tree structure. Johnson solid formations are
a particularly motivating example, given their symmetry, but
other grid and mesh formations are also possible. For instance,
Figure 5 illustrates an example where a curved dome is
represented by the vertices of a tree of squares (Q), and the
symmetric cyclic controller is defined and assembled using
equations (20) and (21). Note that although there are 20 robots
in this formation, each robot only requires relative position
measurements with respect to three other robots (the robots in
the top square require relative measurements with respect to
four other robots) in the group. Thus, in addition to highlighting the applicability of our approach to the formation control
of nonsymmetric formations, this example also underlines the
sparsity of the decentralized formation control laws.
In future work, we intend to investigate strategies for
accommodating local robot failures. For instance, in the regular polygon case, a single point of failure simply requires
the robots to reinitialize their indices within the group and
subsequently converge to a smaller polygon. In the polyhedron
case this can be problematic if the failed robot is coupled

9

to the symmetric cyclic control law introduced in (6), now
expressed as a hybrid dynamical system:
ui (t) =

N
X



km Rm (x)(xi+m −xi )+Rm (x)T (xi−m −xi ) ,

m=1

t ∈ [kτ, (k + 1)τ ), k ∈ N,
(25)

(a) Dome shape formed by the (b) Tree of polygons (squares) Q
convex hull of the robots.
used to derive the formation control laws.

Fig. 5: Dome-shaped formation with 20 robots. The formation
may be extended to collectively film sports or entertainment
events.
to multiple faces. In this scenario, a global restart may be
required with a new configuration that fits the reduced set of
robots.

V. E XTENSIONS
In this section we present a number of extensions to the
symmetric cyclic controller discussed in Sections III and IV.
In Subsection V-A, we provide extensions to control the
formation size, and investigate convergence properties under
delays. Subsection V-B details a heuristic control law for
preventing collisions through use of Repulsive Potential Functions. Finally, in Subsection V-C, we study the disturbance
rejection properties of the symmetric cyclic controller given
in (7) with respect to bounded additive disturbances.

A. Control Over Formation Size
We first address control over the size of a regular polygon
formation. By setting αm = mπ/n in (6), the formation
converges to a regular polygon of fixed, but uncontrolled size
(in particular, the size depends on the initial conditions). In
[16], the authors introduce an additional term to (6) to control
the size of the formation, however, the auxiliary system used
to study convergence using partial contraction does not admit
the necessary stationary zero solution. We present a different
approach to the size control problem by simply modifying the
rotation angle αm embedded within the cyclic controller.
Let ρ be the desired polygon side-length for the formation
and define pi (x) := 1 − kxi+1 − xi k/ρ to be the normalized
error in the inter-robot distance between robots i and i + 1,
where k · k denotes the Euclidean norm. Let p̄(x) be the
arithmetic average over all pi (x) (hereby referred to as the
average inter-robot distance error), and define fs (p̄) to be a
continuous odd scalar function that satisfies the properties: (a)
p̄fs (p̄) > 0 for p̄ 6= 0, (b) |fs (p̄)| ≤ 1, and (c) ∃a ∈ R>0
such that fs0 (0) > a. The last (technical) property is required
to ensure sufficient control authority near p̄ = 0 and will be
formalized in Theorem V.2. Some examples for fs include
the hyperbolic tangent function, or the saturation function
(|fs (pi )| = min{|pi |, 1}). Consider the following modification

where Rm (x) is a rotation matrix around ez with rotation
angle mπ/n+αs (x), where αs (x) := αs0 fs (p̄(x((k −1)τ ))),
αs0 ∈ R>0 is a gain, and τ ∈ R>0 is a constant time lag.
In equation (25) the rotation angle has been modified
equally for all robots by the amount αs , which is a function of
the average inter-robot distance error p̄. The time lag parameter
τ accounts for the finite time required to compute p̄, e.g.,
by using discrete average consensus techniques [26]. That
is, the angle αs used in the time interval [kτ, (k + 1)τ ) is
computed using the average inter-robot distance error at time
(k − 1)τ . To simplify notation, we will use (x, τ ) to reference
the continuous/discrete hybrid dynamics expressed in (25).
As the rotation angle modification is the same for all robots,
the net control law can still be expressed using circulant
matrices:
N
X
ẋ = u = −
km [Lm ⊗ Rm (x, τ ) + LTm ⊗ Rm (x, τ )T ] x.
{z
}
|
m=1
:=Lm (x,τ )

|

{z

}

:=L(x,τ )

(26)
To analyze convergence to the desired formation, first define
Mnρ to be a subset of the original invariant subspace Mn ,
corresponding to the space of regular polygon formations with
inter-robot neighbor distance equal to ρ. Thus,
Mnρ = {x : x ∈ Mn ∧ kxi+1 − xi k = ρ, ∀i}.

(27)

To prove convergence to this subspace, we first start with a
lemma to prove that the subspace Mn is indeed flow-invariant
with respect to (26).
Lemma V.1 (Flow-invariance). The subspace Mn is flowinvariant with respect to the dynamics given in (26).
Proof. Given that the symmetric cyclic controller can still be
expressed using circulant matrices as in (26), by Lemma III.4,
the claim follows.
Having shown that Mn (which contains Mnρ ) is flowinvariant with respect to our modified dynamics, we now
present sufficient conditions to ensure that x converges to the
manifold Mnρ . The proof for the following theorem relies on
several intermediate results which are detailed in the appendix.
Theorem V.2 (Polygon Convergence with Desired Size). Let
C = 2Γ, where
p
N

2(1 − cos(2π/n)) P


km if n is even,

π
sin( n )
m=1
Γ= p
N

2(1 − cos(2π/n)) P


km if n is odd.

π
2 sin( 2n )
m=1
Assume

inf

αm ∈[mπ/n−αs0 ,mπ/n+αs0 ]

min

N
X

1≤i≤n
k∈{−1,0,1} m=1

(m)
km λik


> 0,
(28)

10

(m)

where λik has the form given in Theorem III.5, and suppose
the time lag parameter τ satisfies the following bound:


1
1
,
,
(29)
τ < min
C 8CT
where T denotes a positive constant such that (1/2)T |p̄| ≤
| sin (αs0 fs (p̄)) | ≤ T |p̄| for |p̄| < 1 (if fs is the saturation
function, T = αs0 ).
Then the robots converge to a regular polygon with the
desired side-length ρ. That is, x converges to the manifold
Mnρ .
Remark V.3 (Time Lag Bound). Note that the bound for the
time lag parameter τ given in (29) can be relaxed by using
a value for Γ that is strictly smaller than the one given in
Theorem V.2. The derivation of this relaxed bound is detailed
in the appendix.
Remark V.4 (Control Over Formation Center). The control
laws may be trivially extended to also permit control over
the geometric center of the formation. For instance,
xc
Plet
n
represent the desired center and let x0 := (1/n) i=1 xi
be the instantaneous geometric center. In similar spirit to the
continuous/discrete dynamics in (25), the center control law
uc for robot i may be expressed as:
uc (t) = kc (xc − x0 ((k − 1)τ )) , t ∈ [kτ, (k + 1)τ ), k ∈ N,
where kc ∈ R>0 is a gain and the time lag parameter
τ accounts for the finite time required to compute x0 in
a decentralized fashion (e.g., using message passing). By
leveraging analysis similar to that presented in the proof for
Theorem V.2, it is possible to derive a bound on the steady
state error kxc − x0 k.
Remark V.5 (Size Control Extension to Polyhedra). The
size controller discussed in this subsection may be trivially
extended to the formation control of polyhedra by decoupling
the polygon Q1 from all other polygons Q in the PPS. The net
control for all robots in Q1 is given by the modified symmetric
cyclic controller in (26) plus a “rotational controller” that
rotates the robots within the plane such that the first face has
an orientation compatible with adjacent faces within the PPS.
The control for all other robots in the polyhedron is given
by the standard symmetric cyclic controller as given in (20)
and (21). The proof for convergence is then also decoupled
by treating Q1 as an independent polygon with respect to the
remaining PPS. An example rotational controller along with
a proof of correctness is provided in the appendix.
B. Collision Avoidance
For collision avoidance, one must account for the size of
each robot in the formation. In particular, we assume that each
robot is completely contained within a sphere of radius r1 .
To design a collision avoidance controller, we make use of
Repulsive-Potential-Functions (RPF) [27].
Consider a “collision-detection” zone around each robot
parametrized by two variables r1 < r2 , where r2 denotes the
radius of the boundary of a robot’s detection region. The RPF
between any two robots i and j can be written as follows:
(r2 − r1 )2
+ 2(r2 − r1 ) log(r1 − dij ) − dij
r1 − dij
− v(r1 , r2 ), if r1 < dij ≤ r2 ,
(30)

Vij (dij ) := −

and equals zero for dij > r2 , where dij = kxj − xi k and
v(r1 , r2 ) := 2(r2 − r1 ) log(r1 − r2 ) − r1 is a constant. The
normalized force magnitude can be derived from the RPF by
taking the derivative with respect to dij giving:

undefined
if dij ≤ r1 ,


 (d − r )2
ij
2
fcij (dij ) =
(31)
if r1 < dij ≤ r2 ,

d (d − r1 )2

 ij ij
0
if dij > r2 .
Notice that both Vij and fcij asymptotically approach
as dij → r1 , and are C 1 continuous for dij > r1 . We
formulate the collision avoidance controller for robot
follows:
n
X
ucolli = −
fcij (dij )(xj − xi ).

+∞
now
i as
(32)

j6=i

Thus the collision avoidance controller adds a non-zero (flowinvariant with respect to Mn ) “escape” velocity to the existing
symmetric cyclic controller for robot i if and only if there
exists another robot j within the ith robot’s collision avoidance
zone.
Notice that fcij is unbounded as dij → r1 which makes
convergence analysis with respect to the subspace Mn intractable. In particular, it becomes difficult to upper bound
the maximum eigenvalue of the Jacobian of the collision
avoidance controller. In lieu of this analysis, we present Monte
Carlo simulation results in Section VI for six quadcopters
attempting to converge to a polygon while implementing the
collision avoidance controller above.
C. Robustness Analysis
To conclude our analysis of the symmetric cyclic control
algorithm, we investigate the robustness properties of the
closed-loop dynamics given in (7) in the presence of additive disturbances. We quantify robust performance by the
Euclidean norm of the deviation of the perturbed system
trajectory from the nominal contracting trajectory with respect
to the metric V̄ T V̄ , where V̄ is the orthonormal counterpart of
the nominal projection matrix as defined in (5). The analysis
presented in this section is restricted to zero internal dynamics
(i.e., g(x) = 0) and to convergence to regular polygons. An
analogous proof can be derived for the Johnson polyhedron
case.
Let z(t) = V̄ x(t) represent the dynamics of the unperturbed system under the action of the symmetric cyclic
controller (7). Then,
ż(t) = −V̄ Lx(t).

(33)

Consider the perturbed dynamics zd (t) under an additive state
and/or time-dependant disturbance:
żd (t) = −V̄ Lxd (t) + d(xd , t).

(34)

Here the disturbance is measured in the transformed z coordinates. Let δ(t) represent the difference between the nominal
and perturbed trajectories at time t. Specifically, let δz (t) :=
zd (t) − z(t) and δx (t) := xd (t) − x(t). Note that both δz and
δx equal 0 at t = 0 as the trajectories are assumed to start
from the same initial conditions. Then,
δ˙z (t) = −V̄ L δx (t) + d(xd , t).

(35)

11

Let R̄2 (t) := δzT (t)δz (t); R̄ represents the Euclidean distance
of the perturbed trajectory from the nominal one in the
transformed coordinates, referred to as the formation error.
The following theorem characterizes the robustness of the
symmetric cyclic controller.
Theorem V.6 (Cyclic Control Robustness). Assume that the
nominal system is contracting with contraction rate Λ =
−λmin V̄ LV̄ T < 0 and that the disturbance is norm
¯ Then the formation error in
bounded with upper bound d.
the transformed z coordinates, that is R̄, is upper bounded as
d¯
R̄(t) ≤ (eΛt − 1).
(36)
Λ
Proof. Starting from (35), δx = (V̄ T V̄ + Ū T Ū )δx = V̄ T δz +
Ū T Ū δx . Now, Ū T Ū δx ∈ Mn , so V̄ L (Ū T Ū δx ) = 0. Thus,
δ˙z (t) = −V̄ LV̄ T δz + d(xd , t).
Pre-multiplying by 2δz (t) we obtain:

d(δzT δz )
= −2δzT V̄ LV̄ T δz + 2δzT d(xd , t).
2δzT δ˙z =
dt
Letting R̄2 = δzT δz , one can write

dR̄2
dR̄
= 2R̄
= −2δzT V̄ LV̄ T δz + 2δzT d(xd , t).
dt
dt
Noting that:


−δzT V̄ LV̄ T δz ≤ −λmin V̄ LV̄ T δzT δz = ΛR̄2 ,
and
δzT d(xd , t) ≤ kδz kkd(xd , t)k = R̄kd(xd , t)k,
we get:

R̄R̄˙ ≤ ΛR̄2 + R̄kd(xd , t)k.

(a) Comparison for kδz k.

(b) Bound tightness.

Fig. 6: Tightness for robustness bounds. Parameters: n =
6, N = 2, Λ = −6.928, d¯ = 0.065.
VI. F ORMATION C ONTROL OF Q UADCOPTERS
In this section we demonstrate the use of the symmetric
cyclic control laws developed in Sections III, IV and V-A on
a fleet of quadcopters.
A. Quadcopter Dynamics and Control
Let xi and vi denote the position and inertial velocity
respectively of the ith quadcopter, resolved in an inertial frame
with a downward pointing Z axis. Denote qi to be the attitude
quaternion of the ith quadcopter with associated Euler angles:
Yaw (ψi ), Pitch (θi ), and Roll (φi ) with rotation order ZY X.
The angular velocity of the ith quadcopter with respect to the
inertial frame, resolved in the quadcopter body frame, is given
i
by ωi . Let Rib
denote the rotation matrix to transform vectors
in the inertial frame to the ith quadcopter body frame. Finally
let mqi and Ii denote the mass and moment of inertia tensor,
respectively, for the ith quadcopter.
The translational dynamics of the quadcopter are given by:
ẋi = vi ,

Thus,

  
0
0
 + 0 ,

0
g
−Ti /mqi


R̄˙ ≤ ΛR̄ + kd(xd , t)k,
≤ ΛR̄ + d¯

T

(37a)
(37b)

By the comparison theorem [28], the claim follows.
Remark V.7. The bound R̄(t) describes the time-varying
deviation from a nominal contracting trajectory. In the limit
¯
t → ∞, the steady state bound R̄ss := d/|Λ|
provides a
measure of the steady state formation error. Physically, it
represents a bound on the degree of asymmetry at each node
of the formation, as measured by the set of linear constraints
used to define the subspace, that is, (4a) and (4b). The inverse
proportionality between R̄ss and |Λ| highlights a trade-off
between increased robustness, control saturation, and sensor
requirements.
Figure 6 provides insights into the tightness of the bound.
Here, disturbances were introduced into the system via error
in the control rotation angles. The rotation angle in Rm for
robot i is equal to α̃i (t) = α0 + δαi (t) where α0 = mπ/n is
the nominal control angle and δαi (t) is a perturbation, sampled
randomly from a uniform distribution over the range [δ αi , δ̄αi ]
where δ αi , δ̄αi ∈ [−1, 1]◦ . The plotted curves in Figure 6a
correspond to the actual formation error R̄, the bound in (36),
and the numeric integration of the expression in (37a). Figure
6b quantifies the tightness of (36) and (37a) by plotting the
percentage difference with respect to the actual deviation.

i
v̇i = Rib

(38)

where Ti is the net thrust produced by all four propellers and
g is the gravitational acceleration. Note that the direction of
positive thrust is opposite the body z−axis. The rotational
dynamics of the quadcopter are given by:
 
1 0
q̇i =
qi ,
2 ωi


(39)
−1
ω̇i = Ii
Mi − ωi × Ii ωi ,
where
denotes the quaternion product and Mi is the net
torque on the quadrotor as a result of differential thrust forces
and reaction moments generated by the motor/propeller pairs.
As the dynamics of the motors are significantly faster than
those of the quadcopter as a rigid body, we assume that the
net thrust force and moment generated by the propellers can
be commanded instantaneously.
i
Notice that the rotation matrix Rib
couples the attitude
and translational dynamics, thereby resulting in a fairly complex nonlinear dynamical system. In order to incorporate
the symmetric cyclic control laws for formation control of
quadcopters, we employ a two-level control hierarchy where,
at the high level, the symmetric cyclic control laws in (6),
(21) or (25) are used to generate a desired inertial velocity

12

reference for each quadcopter. The lower level control system
then simply tracks this reference, and incorporates both thrust
and attitude control. Several velocity tracking controllers for
quadcopters have already been presented in the literature [29],
[30]. In the appendix we present an example controller used
for numerical simulations.
B. Numerical Simulations and Experiments
We now present simulation and experimental results demonstrating the application of the above two-level control hierarchy on a group of six quadcopters. The quadcopters were
initialized randomly in hover mode (vi = 0) around the point
[0, 0, −2]T . The size control parameters were: desired interquadcopter distance ρ = 2.0m, function fs (p̄) = tanh(p̄),
and gain αs0 = 5π/180. The desired formation center xc was
selected as the point [0, 0, −10]T m. The symmetric cyclic controller used a look-ahead horizon N = 1 and gain k1 = 0.5.
The desired plane of convergence was tilted 42o with respect to
the horizontal. We assumed a time lag τ = 0.1s. The collision
avoidance parameters were r1 = 0.4m, r2 = 1.2m, and the
controller was slightly modified to take advantage of relative
velocity measurements:
ucolli

n
X
(dij − r2 )2
(xj − xi ),
=−
dij (dij − r1 )2
j6=i

position of the quadcopters within a ball of radius 5m centered
at xc (all simulations initialized the quadcopters in hover
mode). Additionally, upon initialization, the quadcopters were
assigned indices such that their projection onto the desired
plane of convergence formed a non-intersecting polygon. This
is an example of a systematic (and decentralizable) method
for the quadcopters to identify their nearest neighbors for
symmetric cyclic control. In spite of the same speed constraint
of 3m/s, no collisions were recorded while the formation
consistently converged in terms of shape, size and geometric
center.
For experimental validation, we implemented the control
algorithms on a group of four quadcopters, each equipped with
a PX4 Pixhawk autopilot [31]. The low level velocity tracking
controller detailed in the appendix is virtually identical to
the on-board Pixhawk controller. Relative position measurements between neighboring quadcopters were provided by a
Vicon motion capture system. The desired inter-robot neighbor
separation ρ was set to 1.65m. Due to limitations in testing
space, a less aggressive collision avoidance controller was
implemented:
ucolli = −

n
X

kcoll tanh(ρ − dij )(xj − xi ), if dij ≤ r2 ,

j6=i

(40)

if r1 < dij ≤ r2 and vsij < 0,
where vsij = (vj − vi ) · (xj − xi )/dij is the line-of-sight
velocity between robots i and j. Additionally, the desired
velocity (given by the sum of the cyclic, center and collision
avoidance controllers) was subject to saturation constraints
to prevent aggressive maneuvering and increased collision
risk. In particular, each quadcopter was limited to a speed
of 3m/s. Figure 7a shows the resulting 3D trajectories as
well as a shaded plane that outlines the final positions of
the quadcopters. Figure 7b shows the evolution of the interquadcopter distances with respect to quadcopter 1. Notice how
two neighboring quadcopters entered the collision detection
region for quadcopter 1 before converging to the desired
separation.

where kcoll is a gain with value 1.2, and r2 = 0.9m. In
addition, the speed constraint was further restricted to 0.7m/s,
while the cyclic (k1 ) and center (kc ) control gains were set to
0.25 and 0.15 respectively. Note that in all of the experiments
where the quadcopters were initialized outside of the collision
detection regions, the collision avoidance controller did not
become active.
Figures 8a and 8b show the quadcopters in initial and final
(tetrahedron) configurations respectively. Figure 8c shows the
time evolution of the formation error kV̄ xk over all three
phases of flight – autonomous takeoff, formation control, and
landing. Figures 9a and 9b show the time evolution of the
inter-quadcopter distances and the error in geometric center.
The quadcopters demonstrated relatively quick convergence in
terms of formation shape, size and center, thereby validating
the two-level formation control methodology 1 .
VII. C ONCLUSIONS

(a) Quadcopter 3D trajectories (b) Inter-quadcopter distances (in
with size and center control. metres) for Quadcopter 1.

Fig. 7: Planar formation simulation.
To test the collision avoidance controller, we ran a Monte
Carlo simulation with 100 samples randomizing the initial

In this paper we studied decentralized control algorithms
for 3D formations such as regular polygons, Johnson solids,
and more generally, connected regular polygonal meshes. The
algorithms are desirable for their simplicity in that they only
require relative position measurements between neighbors and
an agreement between the robots on the desired orientation in
space. Our approach leveraged the mathematical properties of
cyclic pursuit along with results from contraction and partial
contraction theory. In addition, we provided extensions to the
algorithm for controlling the formation size and preventing
collisions, and quantified the robustness of the symmetric
cyclic controller under bounded additive disturbances. Finally,
we validated our algorithms via numerical simulations and
hardware experiments on a fleet of quadcopters.
This work allows several possible avenues for extension.
First, we would like to study (sub-)optimality properties of
1 Videos of the experiments are available at https://www.youtube.com/
playlist?list=PL8-2mtIlFIJr4fVjJeyQqymbkKfhsKKYX

13

(a) Quadcopters upon initialization in hover mode.

(a) Error in geometric center.

(b) Quadcopters in final tetrahedron configuration.

(b) Inter-quadcopter neighbor separations. ρ = 1.65m.

Fig. 9: Formation accuracy.

(c) Tetrahedron formation error. The three shaded regions indicate the flight modes (from left to right: takeoff, formation
control, landing).

Fig. 8: Quadcopter tetrahedron formation.
the symmetric cyclic controller by introducing an objective
cost function to be minimized (e.g., time or control effort),
and using it to guide gain selection and/or look-ahead horizon.
Second, we plan to investigate various examples for internal
dynamics g(x) and their interplay with the symmetric cyclic
controller when analyzing convergence. Third, the modification of the symmetric cyclic controller for formation size is
a little cumbersome in that it requires all robots to compute
a shared quantity. Consequently, we would like to investigate
size controllers that do not rely upon consensus techniques.
R EFERENCES
[1] W. Sheng, Q. Yang, J. Tan, and N. Xi, “Distributed Multi-Robot
Coordination in Area Exploration,” Robotics and Autonomous Systems,
vol. 54, no. 12, pp. 945 – 955, 2006.
[2] R. Smith and F. Hadaegh, “Control topologies for deep space formation
flying spacecraft,” in American Control Conference, vol. 4, 2002, pp.
2836–2841.
[3] F. Arrichiello, H. Heidarsson, S. Chiaverini, and G. Sukhatme, “Cooperative caging using autonomous aquatic surface vehicles,” in Proc. IEEE
Conf. on Robotics and Automation, May 2010, pp. 4763–4769.
[4] A. Jadbabaie, J. Lin, and A. Morse, “Coordination of groups of mobile
autonomous agents using nearest neighbor rules,” IEEE Transactions on
Automatic Control, vol. 48, no. 6, pp. 988–1001, June 2003.

[5] J. Lawton, R. Beard, and B. Young, “A decentralized approach to
formation maneuvers,” IEEE Transactions on Robotics and Automation,
vol. 19, no. 6, pp. 933–941, Dec 2003.
[6] W. Ren and R. Beard, “Decentralized Scheme for Spacecraft Formation
Flying via the Virtual Structure Approach,” AIAA Journal of Guidance,
Control, and Dynamics, vol. 27, no. 1, pp. 73–82, Jan-Feb 2004.
[7] J. Marshall, “Coordinated Autonomy: Pursuit Formations of Multivehicle Systems,” Ph.D. dissertation, University of Toronto, 2005.
[8] J. Chen, D. Sun, J. Yang, and H. Chen, “Leader-Follower Formation
Control of Multiple Non-holonomic Mobile Robots Incorporating a
Receding-horizon Scheme,” International Journal of Robotics Research,
vol. 29, no. 6, pp. 727–747, 2010.
[9] L. Consolini, F. Morbidi, D. Prattichizzo, and M. Tosques, “Leader
follower formation control of nonholonomic mobile robots with input
constraints,” Automatica, vol. 44, no. 5, pp. 1343 – 1349, 2008.
[10] D. P. Scharf, F. Y. Hadaegh, and S. R. Ploen, “A Survey of Spacecraft
Formation Flying Guidance and Control, Part II: Control,” in American
Control Conference, vol. 4, Boston, MA, Jun. 2004, pp. 2976–2985.
[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, vol. 9, no. 1, pp. 427–438, Feb
2013.
[12] F. Guerriero, R. Surace, V. Loscri, and E. Natalizio, “A multi-objective
approach for unmanned aerial vehicle routing problem with soft time
windows constraints,” Applied Mathematical Modelling, vol. 38, no. 3,
pp. 839–852, 2014.
[13] J. Marshall, M. Broucke, and B. Francis, “Formations of vehicles in
cyclic pursuit,” IEEE Transactions on Automatic Control, vol. 49, no. 11,
pp. 1963–1974, Nov 2004.
[14] M. Pavone and E. Frazzoli, “Decentralized Policies for Geometric
Pattern Formation and Path Coverage,” ASME Journal on Dynamic
Systems, Measurement, and Control, vol. 129, no. 5, pp. 633–643, Oct.
2007. [Online]. Available: http://dynamicsystems.asmedigitalcollection.
asme.org/article.aspx?articleid=1412612
[15] J. L. Ramirez, M. Pavone, E. Frazzoli, and D. W. Miller, “Distributed
Control of Spacecraft Formations via Cyclic Pursuit: Theory and
Experiments,” AIAA Journal of Guidance, Control, and Dynamics,

14

[16]
[17]
[18]
[19]
[20]
[21]

[22]

[23]
[24]
[25]
[26]
[27]

[28]
[29]
[30]
[31]
[32]

vol. 33, no. 5, pp. 1655–1669, Sep. 2010. [Online]. Available:
http://arc.aiaa.org/doi/abs/10.2514/1.46511
J. Ramirez-Riberos and J. Slotine, “Contraction theory approach to generalized decentralized cyclic algorithms for global formation acquisition
and control,” in Proc. IEEE Conf. on Decision and Control, Dec 2012.
W. Lohmiller and J. Slotine, “On Contraction Analysis for Non-linear
Systems,” Automatica, vol. 34, no. 6, pp. 683 – 696, 1998.
Q. Pham, N. Tabareau, and J. Slotine, “A Contraction Theory Approach
to Stochastic Incremental Stability,” IEEE Transactions on Automatic
Control, vol. 54, no. 4, pp. 816–820, April 2009.
W. Wang and J. Slotine, “On partial contraction analysis for coupled
nonlinear oscillators,” Biological Cybernetics, vol. 92, no. 1, pp. 38–53,
2005.
Q. Pham and J. Slotine, “Stable concurrent synchronization in dynamic
system networks,” Neural Networks, pp. 62–77, 2007.
O. Wallner, K. Ergenzinger, and R. Flatscher, “X-array aperture configuration in planar or non-planar spacecraft formation for Darwin/TPF-I
candidate architectures,” in Proceedings of SPIE 6693, Techniques and
Instrumentation for Detection of Exoplanets III, Sep 2007.
S. Singh, E. Schmerling, and M. Pavone, “Decentralized Algorithms for
3D Symmetric Formations in Robotic Networks - A Contraction Theory
Approach,” in Proc. IEEE Conf. on Robotics and Automation, Seattle,
WA, May 2015.
J. Jouffroy and J. Slotine, “Methodological remarks on contraction
theory,” in Proc. IEEE Conf. on Decision and Control, vol. 3, Dec 2004,
pp. 2537–2543 Vol.3.
A. Alexandrov, Convex Polyhedra, 1st ed.
Springer-Verlag Berlin
Heidelberg, 2005, english translation with additional material of the
1950 Russian monograph.
M. J. Wenninger, Polyhedron models. Cambridge University Press,
1974.
N. A. Lynch, Distributed Algorithms. San Francisco, CA, USA: Morgan
Kaufmann Publishers Inc., 1996.
D. H. Kim, H. Wang, and S. Shin, “Decentralized control of autonomous
swarm systems using artificial potential functions: Analytical design
guidelines,” Journal of Intelligent and Robotic Systems, vol. 45, no. 4,
pp. 369–394, 2006.
H. K. Khalil, Nonlinear Systems, 3rd ed. Prentice Hall, 2002.
R. Mahony, V. Kumar, and P. Corke, “Multirotor aerial vehicles:
Modeling, estimation, and control of quadrotor,” IEEE Robotics and
Automation Magazine, no. 19, pp. 20–32, 2012.
A. Kushleyev, D. Mellinger, C. Powers, and V. Kumar, “Towards a
swarm of agile micro quadrotors,” Autonomous Robots, vol. 35, no. 4,
pp. 287–300, 2013.
Pixhawk. Px4 autopilot project. [Online]. Available: http://pixhawk.org
R. Gray, “Toeplitz and Circulant Matrices: A review,” Tech. Rep., 2001.

Sumeet Singh is is a Ph.D. candidate in the Autonomous Systems Laboratory at Stanford University. He received a B.Eng. in Mechanical Engineering and a Diploma of Music (Performance)
from the University of Melbourne in 2012, and a
M.Sc. in Aeronautics and Astronautics from Stanford University in 2015. Prior to joining Stanford,
Sumeet worked in the Berkeley Micromechanical
Analysis and Design lab at University of California
Berkeley in 2011, and the Aeromechanics Branch
at NASA Ames in 2013. Sumeet’s current research
interests include decentralized algorithms for formation control, and stochastic
and robust nonlinear model predictive control, with specific applications to
spacecraft control.

Edward Schmerling received the B.Sc. degree in
mathematics and physics from Stanford University, Stanford, CA, USA, in 2010. He is currently
working towards the Ph.D. degree at the Institute
for Computational and Mathematical Engineering,
Stanford University. His research interests include
developing algorithms for robotic motion planning
subject to general dynamic and cost constraints.

Dr. Marco Pavone is an Assistant Professor of
Aeronautics and Astronautics at Stanford University,
where he is the Director of the Autonomous Systems
Laboratory. Before joining Stanford, he was a Research Technologist within the Robotics Section at
the NASA Jet Propulsion Laboratory. He received a
Ph.D. degree in Aeronautics and Astronautics from
the Massachusetts Institute of Technology in 2010.
Dr. Pavone’s areas of expertise lie in the fields of
controls and robotics. His main research interests
are in the development of methodologies for the
analysis, design, and control of autonomous systems, with an emphasis on
autonomous aerospace vehicles and large-scale robotic networks. He is a
recipient of an NSF CAREER Award, a NASA Early Career Faculty Award,
a Hellman Faculty Scholar Award, and was named NASA NIAC Fellow in
2011. He is currently serving as an Associate Editor for the IEEE Control
Systems Magazine.

15

A PPENDIX
VIII. M ATHEMATIC P RELIMINARIES
In this section we present the eigendecomposition of circulant matrices. The following theorem will be instrumental to
the proof for Theorem III.5.
Theorem VIII.1 (Adapted from Theorem 3.1 in [32]). The
eigenvalues of a circulant matrix are the Discrete Fourier
Transform of the first row of the matrix. Specifically, every
circulant matrix C = circ[c1 c2 . . . cn ] has eigenvectors:
1
vk = √ (1, e−2πj(k−1)/n , . . . , e−2πj(k−1)(n−1)/n )T ,
n
(41)
k = 1, 2, . . . , n,
and corresponding eigenvalues:
λk =

n
X

cp e2πj(k−1)(p−1)/n ,

(42)

p=1

and can be expressed in the form C = U ΛU ∗ where U
is a unitary matrix with the k th column equal to the k th
eigenvector and Λ is a diagonal matrix of corresponding
eigenvalues. Note that the same matrix U diagonalizes all circulant matrices. Thus if C and B are n×n circulant matrices,
with eigenvalues {λB,k }nk=1 and {λC,k }nk=1 respectively, then:
1) C and B commute (CB = BC) and CB is also a circulant matrix with eigenvalues equal to {λB,k λC,k }nk=1
2) C+B is also circulant with eigenvalues {λB,k +λC,k }nk=1
IX. P ROOFS
A. Proofs for Section III
Proof of Lemma III.1. Necessity is trivial. We then consider
sufficiency. First, we prove that constraints (4a) and (4b)
ensure that all robots lie in a common plane, i.e.,
eTz xi = a,

to ensure that the robots do not form a spiral formation out
of the desired plane, thereby reducing the dimensionality of
the solution subspace. Thus, the in-plane constraint must be
independent to the rotational constraints.
Proof of Lemma III.3. Equations (4a) and (4b) can be written
in matrix form as V x = 0, where

I3
−(I3 + R2π/n )
R2π/n
 03×3
I
−(I
+ R2π/n )
3
3

...
V =
..
..
..

.
.
.
eTz R2π/n
01×3
01×3

03×3
··· ···
03×3

..

R2π/n 03×3 · · ·
.
.

..
..
..
..

.
.
.
.
···
···
eTz −eTz (I3 + R2π/n )
(45)
Simplifying the above expression, one obtains
V =Wn circ[I3 , −I3 , 03×3 , . . . , 03×3 ]+
(46)
Wn circ[03×3 , −R2π/n , R2π/n , 03×3 , . . . , 03×3 ]
=Wn (L1 ⊗ I3 ) + Wn ((L1 − L2 ) ⊗ R2π/n I3 )

=Wn L1 In ⊗ I3 I3 + (L1 − L2 )In ⊗ R2π/n I3
=Wn (L1 ⊗ I3 + (L1 − L2 ) ⊗ R2π/n )(In ⊗ I3 ),
where the last equality follows from the properties of the
Kronecker product. The claim then follows.
Proof of Lemma III.4. Note that the rotational constraints in
(4a) and the in-plane constraint (4b) define the subspace Mn
to be a regular polygon with normal ez . The following analysis
therefore assumes that all robots lie in a single plane and
satisfy the rotational constraints. Then,

(43)

for all i ∈ {1, . . . , n} and some a ∈ R. Since ez is a left
eigenvector of R2π/n , with eigenvalue equal to 1, the rotational
constraints imply eTz (xi+1 −xi ) = eTz (xi+2 −xi+1 ), for i =
1, . . . , n−2. Combining this set of equations with the in-plane
constraint, one can write

ẋi+1 − ẋi =

eTz xi+2

=

eTz x2

+

i eTz (x2

− x1 ),

(44)

for i = 1, . . . , n − 1 and where the indices are modulo n.
Hence, for i = n − 1, one obtains


T
+ Rm
(x, t)(xi+1−m − xi+1 )
−

which implies that
− x1 ) = 0. Hence, by setting
eTz x2 := a, one immediately obtains equation (43). We have
then proven that all robots lie in a common plane. Now, for
n points lying in a common plane, the rotational constraints
represent the definition of a regular polygon.
For independence, we note that each successive rotational
constraint involves at least one robot that does not belong to
any of the previous constraints, proving that the rotational constraints are independent. The in-plane constraint is necessary

N
X


km (x, t) Rm (x, t)(xi+m − xi )

m=1

T
+ Rm
(x, t)(xi−m − xi ) .

(47)

Now by the constraints described by (4a), we can state:
xi+1+m − xi+1 = R2π/n (xi+2+m − xi+2 ), and
xi+1−m − xi+1 = R2π/n (xi+2−m − xi+2 ).
Thus, the first term in (47) can be written as:

eTz x1 = eTz x2 + (n − 1) eTz (x2 − x1 ),
eTz (x2


km (x, t) Rm (x, t)(xi+1+m − xi+1 )

m=1

eTz (xi+1 −xi ) = eTz (xi+2 −xi+1 ),
for i = 1, . . . , n − 1 and where the indices are modulo n. One
can then readily show that

N
X

ẋi+1 =

N
X


km (x, t) Rm R2π/n (x, t)(xi+2+m − xi+2 )

m=1


T
+ Rm
R2π/n (x, t)(xi+2−m − xi+2 )
=R2π/n

N
X


km (x, t) Rm (x, t)(xi+2+m − xi+2 )

m=1
T
Rm (x, t)(xi+2−m

+
=R2π/n ẋi+2 ,

− xi+2 )


(48)

16

where the order of Rm and R2π/n can be swapped since they
are rotation matrices about the same axis. Similarly, it can be
shown that ẋi = R2π/n ẋi+1 . Thus, rewriting (47):
ẋi+1 − ẋi = R2π/n (ẋi+2 − ẋi+1 ).
To prove the arbitrary normal case, one need only replace
Rm and R2π/n with their similarity transformed versions with
respect to Rη 6= I3 .
Proof of Theorem III.5. We provide the proof for the general
case where the desired formation plane is not necessarily the
horizontal plane. The proof leverages the following Lemma.
Lemma IX.1. Given a symmetric matrix X such that X ≺ 0,
then Wn X WnT ≺ 0.
Proof. If X ≺ 0 then λmax (X ) < 0. Noting that
Wn is sub-unitary, by the Cauchy Interlacing Theorem,
λmax (Wn X WnT ) ≤ λmax (X ) < 0.
Having proved Lemma IX.1, we are now in a position
to prove Theorem III.5. Specifically, the strategy is to only
consider the uniform negative-definiteness of the inner term
(not involving Wn ) in (8).
For the scenario where the desired formation plane is not
the horizontal plane, by the discussion in Remark III.7, we
must use the similarity transformed versions of the projection
matrix V and the symmetric cyclic controller. From (10) and
(11), note that Pn Rη Lmη RTη PnT = Pn Lm PnT since Lmη =
RTη Lm Rη . Thus, the analysis does not depend on the arbitrary
choice of the orientation of the desired plane normal in the
global coordinate system.
Now, given the matrices Pn and Lm are ∈ CL⊗R , they have
the same set of eigenvectors. Then, for any eigenvector vi in
this set, we have the following relation for the corresponding
eigenvalues:

λi Pn Lm PnT = λi (Pn )λi (PnT )λi (Lm ).
To derive the eigenvalues of Pn , we note that Pn =
L1 ⊗ I3 + (L1 − L2 ) ⊗ R2π/n which is the sum of two CL⊗R
matrices and thus the eigenvalues of Pn must be the sum of
the eigenvalues of L1 ⊗ I3 and (L1 − L2 ) ⊗ R2π/n . From
Theorem VIII.1, we have the general result:
λi (Lm ) = 1 − e

2mπ
n (i−1)j

, i = 1, . . . , n.

(49)

Additionally, from Theorem VIII.1, the eigenvalues of L1 −L2
are the sum of eigenvalues of L1 and −L2 . The eigenvalues of
2kπ
2π
R2π/n are {1, e± n j } = e n j where k ∈ {−1, 0, 1}. Thus,
the eigenvalues of Pn and PnT are the set i ∈ {1, . . . , n}, k ∈
{−1, 0, 1}:


2π
λ(Pn , PnT ) = 1 − e± n (i−1)j
 4π
 2kπ
2π
+ e± n (i−1)j − e± n (i−1)j e± n j , (50)
where the positive and negative signs are for Pn and PnT
respectively. In a similar fashion, the eigenvalues for Lm are:


2mπ
λ(Lm ) = 1 − e n (i−1)j ekαm j


−2mπ
+ 1 − e n (i−1)j e−kαm j



2πm(i − 1)
=2 cos(kαm ) − cos kαm +
.
n
(51)

Multiplying (50) and (51) gives the expression for
(m)
λ Pn Lm PnT = λik (x, t). To obtain the summation form
in(9), note that PnLm PnT ∈ CL⊗R for all m. Thus
N
N


P
P
λ Pn
km Lm PnT =
km λ Pn Lm PnT . Having
m=1

m=1

obtained the summation form, (9) follows directly from trying
to show (8) given Lemma IX.1.

B. Proofs for Section IV
Proof of Corollary IV.5. Necessity is straightforward. We then
consider sufficiency. Let us assume without loss of generality
that nk = ez ; the general case follows through a rotation
under Rηk . To prove the claim, we need to show that the
rotational constraints are sufficient to ensure that all robots
lie in the desired plane for j = 1, . . . , |Vk |. The claim then
follows directly. Now, the rotational constraints imply
eTz (xi+1 −xi ) = eTz (xi+2 −xi+1 ),

(52)

for i = 1, . . . , |Vk |−2. For two robots with indices j and j +1
that lie in the desired plane, we know that eTz xj = eTz xj+1 =
a ∈ R. If j ∈ {1, . . . , |Vk | − 1}, we can recursively use (52)
to show
eTz xi = eTz xj ,
(53)
for i = 1, . . . , |Vk |. For the case where j = |Vk |, we have:
eTz xn = eTz x1 .
However, we note that (44) still holds for i = 1, . . . , |Vk | − 2
using the rotational constraints alone. Then, setting i = |Vk | −
2, we obtain
(eTz xn −eTz x1 )+(|Vk |−1)eTz (x2 −x1 ) = (|Vk |−1)eTz (x2 −x1 ),
which gives eTz x2 = eTz x1 . We can now again use (52)
recursively to show (53), completing the proof.
C. Proofs for Section V
1) Control Over Formation Size: In order to prove Theorem
V.2, we first require several intermediate results detailed in
the following lemmas. We first prove that x converges to the
subspace Mn .
Lemma IX.2 (Convergence to Mn with Size Control). Assume


N
X
(m)
inf
min
km λik
> 0,
αm ∈[mπ/n−αs0 ,mπ/n+αs0 ]

1≤i≤n
k∈{−1,0,1} m=1

(54)
(m)

where λik has the form given in Theorem III.5. Then the
system ẋ = u under the controller (26) globally converges to
Mn .
Proof. To show convergence to Mn , consider the following
auxiliary system for the transformed closed-loop dynamics
żp = V̄ x˙p = V̄ u where u is as given in (26):
ẏ = −V̄ L(xp , τ ) (V̄ T y + Ū T Ū xp ),

17

Note that y = 0 is a solution of this auxiliary system
since V̄ L(xp , τ ) (Ū T Ū xp ) = 0, ∀xp . Then to show that the
auxiliary system is contracting in y, we require:


V̄ L(xp , τ ) V̄ T  0,
uniformly, ∀xp , τ.
By Remark II.4, the above stability requirement can be reformulated in terms of V directly, i.e.,


Wn Pn L(xp , τ ) PnT WnT  0,

uniformly, ∀xp , τ. (55)

The above condition is equivalent to (8) with g = 0,
L = L(xp , τ ) and the additional stipulation that (55) holds
for all possible xp and time lags τ . Temporarily neglecting
this additional condition and applying Theorem III.5 with
g = 0 gives the expression within the infimum in (54). Then
to account for all possible xp and time lags τ , we note that the
dependance of L on xp and τ stems from the rotation angle
αm = mπ/n + αs (xp , τ ), which changes every τ seconds. As
|fs (p̄)| ≤ 1, it follows that αm ∈ [mπ/n − αs0 , mπ/n + αs0 ],
for all time (despite the discrete switches every τ seconds).
Choosing the worst case value for αm within this interval
then completes the expression given in (54).
Remark IX.3. Note that the continuous/discrete hybrid dynamics induced by the time lag parameter τ do not affect the
convergence analysis for Mn . By taking the infimum in (54),
we have that for any modification of the standard rotation
angle mπ/n by the angle αs , condition (55) is still satisfied.
We now study the behavior of the system within the
invariant subspace Mn , with respect to the desired subset
Mnρ . Here, the time lag τ will play a crucial role in governing
the convergence of the inter-robot separation to the desired
distance.
Lemma IX.4 (Formation Size Discrete p
Dynamics). Assume
x lies in the subspace Mn . Let β = 2(1 − cos(2π/n)),
N
P
Γ=
km γm where γm are fixed constants derivable from
m=1

geometry, specifically,

γ1 = 1, and γm



1


 1,
sin( nπ ) 
∈ 
1


 1,
π
2 sin( 2n
)

if n is even,

ẋi = 2 sin(αs )
ẋi+1 = 2 sin(αs )
= 2 sin(αs )

N
X
m=1
N
X
m=1
N
X

!
km kxi+m − xi k x̂oi
!
km kxi+1+m − xi+1 k x̂oi+1
!
km kxi+m − xi k x̂oi+1 ,

m=1

where x̂oi is the unit vector pointing outward from the
instantaneous geometric center x0 towards robot i. Exploiting
symmetry, we can re-write kxi+m − xi k as γm kxi+1 − xi k,
where γ1 = 1 and γm for m = 2, . . . , N may be derived from
geometry. Using these simplifications, we obtain



:=Γ>0

for m = 2, . . . , N , and let C = 2βΓ. Then, under the action
of the modified cyclic controller given in (26), each inter-robot
distance error pi admits the following recursive solution:
pik+1 = (pik − 1)e

To evaluate (ẋi+1 − ẋi ), consider Figure 10.
From Remark III.6, we see the geometric realization of
the forward and rear components of the cyclic controller
cancelling for x ∈ Mn , when αm = mπ/n. We can
also see that the vector Rm (xi+m − xi ) where Rm is a
rotation matrix with angle mπ/n + αs , is simply a rotation of
Rmπ/n (xi+m − xi ) about the plane normal by the angle αs .
Thus, we can write:

X

 N


(ẋi+1 − ẋi ) = 2 sin(αs ) 
km γm 
 kxi+1 − xi k
m=1

| {z }

if n is odd,

C sin(αs (pik−1 ))τ

Fig. 10: Cyclic control with αs

+ 1,

(56)

where pik = pi (kτ ).
Proof. First note that by Lemma V.1, Mn is flow-invariant
with respect to the control law given in (26). Thus, we will
exploit the symmetry properties attributed to regular polygons.
In particular, recall the definition pi := 1 − kxi+1 − xi k/ρ.
We note that for x ∈ Mn , p1 = p2 = . . . = pn = p̄. Consider
the time-derivative of pi within the time-range [kτ, (k + 1)τ ):
(xi+1 − xi )T (ẋi+1 − ẋi )
ṗi = −
.
ρkxi+1 − xi k

×

p
|

(xi+1 − xi )
2(1 − cos(2π/n))
.
{z
} kxi+1 − xi k
:=β>0

We now recall that αs = αs (p̄k−1 ), and is constant over the
time range [kτ, (k + 1)τ ). Thus, we can write ṗi as:
 kxi+1 − xi k
ṗi = −2βΓ sin αs (p̄k−1 )
ρ

= C sin αs (pik−1 ) (pi − 1).

(57)

Solving the above equation exactly gives the recursive formula
given in (56).
We are now ready to derive sufficient conditions for the time
lag parameter τ to ensure that the recursive solution given in
(56) converges to 0.

18

Lemma IX.5 (Convergence to Mnρ within Mn ). Assume x
lies in the subspace Mn . Let T denote a positive constant
such that (1/2)T |pi | ≤ | sin (αs0 fs (pi )) | ≤ T |pi | for |pi | < 1
(if fs is the saturation function, T = αs0 ). Suppose the time
lag parameter τ satisfies the following bound:


1
1
τ < min
,
.
(58)
C 8CT

By the comparison theorem [28], we obtain for t ∈ [k 0 τ, (k 0 +
1)τ ]:

Then the neighbor inter-robot distance governed by the discrete dynamics given in (56) converges to the desired distance
ρ. That is, x converges to the subspace Mnρ .

pik0 +1 < ξ((k 0 + 1)τ ).

Proof. The layout of the proof is as follows: we first show
that there exists some time index k 0 such that the inter-robot
distance error pik has the same sign in two consecutive time
steps. Next, we prove that if pik has the same sign for all k >
k 0 , then pik converges to zero. Finally, we provide sufficient
conditions to ensure that either (a) pik possesses the same sign
for all k > k 0 , or (b) if there exists a time index k 00 where pik
changes sign, then it never crosses 0 again. We begin with the
following claim:
Claim IX.6. Given the inter-robot distance error discrete
dynamics given in (56), there exists a time index k 0 such that
pik0 −1 pik0 > 0. That is, the distance error has the same sign
in two consecutive time steps.
Proof. Assume by way of contradiction that pik−1 pik < 0,
pik−1 < 0 and 1 > pik > 0. Then, sin(αs (pik−1 )) < 0.
From (56), we have 1 > pik+1 > pik > 0. Alternatively,
assume that pik−1 pik < 0, pik−1 > 0 and pik < 0. Then,
sin(αs (pik−1 )) > 0. Using (56), we have pik+1 < pik < 0.
Setting k 0 = k + 1 gives the desired contradiction.
Given there exists a time index k 0 , such that pik0 −1 pik0 > 0,
consider the following result:
Claim IX.7. Let k 0 be a time index such that pik0 −1 pik0 > 0.
If
1) pik0 > 0 and pik > 0 for all k > k 0 , or
2) pik0 < 0 and pik < 0 for all k > k 0 ,
then pik converges to 0.
Proof. Consider the first case. Without loss of generality, assume pik0 −1 > pik0 . If this is not true, then take k 00 to be k 0 +1.
Then given (56), we have pik00 = pik0 +1 < pik0 = pik00 −1 , and
we proceed with our analysis for k > k 00 .
Now if pik > 0 for all k > k 0 , and pik pik−1 > 0 for all k >
0
k , then given (56), we deduce that pik is a strictly decreasing
monotonic sequence lower-bounded by 0, and therefore has a
limit l ≥ 0. We will show that l = 0.
Consider the continuous time differential equation in (57).
Given the recursive solution in (56), we know that pi (t) <
pik0 −1 (and consequently αs (pi (t)) < αs (pik0 −1 )) for all t >
(k 0 − 1)τ . Furthermore, since pi (t) − 1 < 0, we obtain the
following differential inequality for t ∈ [k 0 τ, (k 0 + 1)τ ]:


p˙i (t) = C sin αs (pik0 −1 ) (pi (t) − 1)
(59)
< C sin (αs (pi (t))) (pi (t) − 1).
We now use the lower bound: sin (αs (pi (t))) > (1/2)T pi (t).
It follows then:
CT
p˙i (t) <
pi (t)(pi (t) − 1).
2

1

pi (t) <


1+

1
p ik 0


:= ξ(t).
CT
0
− 1 e 2 (t−k τ )

In particular, we note that

Proceeding inductively, it follows that
pik < ξ(kτ ), ∀k > k 0 .
That is, the recursive solution is bounded above by ξ(t), a
Class-L function [28], for all k > k 0 . Consequently, the limit
l for pik must equal 0.
Consider now the second case. Without loss of generality,
assume pik0 > pik0 −1 . If this is not true, then take k 00 to be
k 0 + 1. Then given (56), we have pik00 = pik0 +1 > pik0 =
pik00 −1 , and we proceed with our analysis for k > k 00 .
Now if pik < 0 for all k > k 0 and pik pik−1 > 0 for all k >
k 0 , then given (56), we deduce that pik is a strictly increasing
monotonic sequence upper-bounded by 0, and therefore has a
limit −l ≤ 0. We will show that l = 0.
We first prove that l < 1. That is, there exists a time index
k ∗ such that for all k > k ∗ , |pik | < 1. By way of contradiction,
assume that l ≥ 1 and suppose pik = −l −  where  > 0
is an infinitesimally small positive number. Given (56), we
know that pik−1 < pik < −1. Then since Cτ > 0, we know
that C sin(αs (pik−1 ))τ < C sin(αs0 fs (−1))τ < 0, and thus
eC sin(αs (pik−1 ))τ < 1. Since the sequence pik is increasing
monotonically, for the limit to be −l, pik < −l for all k. To
obtain the desired contradiction, we will show that there exists
a small enough  such that pik+1 > −l. Using (56), we seek
an  > 0 such that:
(−l −  − 1)eC sin(αs (pik−1 ))τ > −l − 1.
Re-arranging the above expression, we require:
eC sin(αs (pik−1 ))τ <

l+1
,
l+1+

Since eC sin(αs (pik−1 ))τ < eC sin(αs0 fs (−1))τ < 1, there exists
an  > 0 small enough such that the above condition will be
satisfied. We now proceed with our analysis for k > k ∗ . For
t > k ∗ τ , we have that |pi (t)| < |pik∗ |. Then, using (57) and
that |pi (t)| < 1, we obtain the following differential inequality
for t ∈ [(k ∗ + 1)τ, (k ∗ + 2)τ ]:

p˙i (t) = −C sin αs (pik∗ ) (pi (t) − 1)
> −C |sin (αs (pi (t)))| (pi (t) − 1).
Using the lower bound: |sin (αs (pi (t)))| > (1/2)T |pi (t)| for
t > k ∗ τ , it follows:
CT
|pi (t)|(pi (t) − 1).
2
Once again, applying the Comparison Theorem [28], we obtain
for t ∈ [(k ∗ + 1)τ, (k ∗ + 2)τ ]:
p˙i (t) > −

|pi (t)| <
1+



1
pik∗ +1

1
 CT
:= η(t).
∗
− 1 e 2 (t−(k +1)τ )

19

that xM → Mnρ . To prove this result, let us examine the
dynamics of xM :

We note that:
|pik∗ +2 | < η((k ∗ + 2)τ ).

ẋM = Ū T Ū u = −Ū T Ū L(x, τ )x

Proceeding inductively, it follows that

= −L(x, τ )Ū T Ū x
= −L(x, τ )xM
= −L(xM , τ )xM − (L(x, τ ) − L(xM , τ )) xM ,
{z
}
|

|pik | < η(kτ ), ∀k > k ∗ + 1.
Thus, |pik | is bounded above by η(t), a Class-L function [28],
for all k > k ∗ +1. Thus, the limit −l for pik must equal 0.
Now suppose pik0 −1 pik0 > 0 and pik0 > 0. We now derive
sufficient conditions to ensure that pik > 0 for all k > k 0 . Let
k = k 0 . Since τ < 1/C, and 0 < sin(αs (pik−1 )) < 1, we have
0 < C sin(αs (pik−1 ))τ < 1. Using the inequality eζ ≤ 1 + 2ζ
for ζ ∈ [0, 1], we may deduce:
eC sin(αs (pik−1 ))τ ≤ 1 + 2C sin(αs (pik−1 ))τ.

Additionally, since | sin αs0 fs (pik−1 ) | ≤ T |pik−1 |, we have:
1 + 2C sin(αs (pik−1 ))τ ≤ 1 + 2CT
| {z τ} pik−1 .
:=A>0

Since, pik − 1 < 0, it follows

pik+1 ≥ (pik −1) 1 + Apik−1 +1 = pik +Apik−1 pik −Apik−1 .
Equivalently,
pi
pik+1
≥ 1 + Apik−1 − A k−1 .
pik
p ik
Note that pik /pik−1 > 1 > 1/2 (i.e. pik−1 /pik < 2) and
Apik−1 > 0. Then, if Apik−1 /pik < 1/2 (i.e. A < 1/4), then,
0 < (1/2)pik < pik+1 < pik .
By induction, we conclude that pik > 0 for all k > k 0 . By
Claim IX.7, pik converges to 0. The requirement A < 1/4 is
guaranteed by (58).
Finally, we address the scenario where pik0 −1 pik0 > 0 and
pik0 < 0. There are two sub-cases:
1) For all k > k 0 , pik < 0. Then by Claim IX.7 above, pik
converges to 0.
2) There exists a time index k 00 such that pik00 −1 pik00 < 0
(i.e. pik00 > 0).
Analyzing the second sub-case, using (56), we have that
pik00 +1 > pik00 > 0. Then given A < 1/4, pik > 0 for all
k > k 00 . Thus, pik does not change sign again, completing the
proof.
We are now ready to prove Theorem V.2 which links
convergence analysis both on and off the subspace Mn .
Proof of Theorem V.2. By Lemmas V.1, IX.4 and IX.5, we
know that if x(0) lies in the flow-invariant subspace Mn ,
then x converges to the manifold Mnρ .
We now address the scenario where x(0) ∈ R3n \ Mn . The
state x may be decomposed into projections both on and off
Mn as follows:
T
x = Ū
| {zŪ x} +
:=xM ∈Mn

T
V̄
| {zV̄ x}

.

:=xM ⊥ ∈M⊥
n

By Lemma IX.2, we have that kxM ⊥ k is decreasing exponentially to 0. Thus x → xM . All that remains is to show

:=up

where the order of Ū T Ū and L can be swapped since the
order of cyclic control and projection onto the subspace Mn
is irrelevant. Note that from the analysis conducted in Lemma
IX.4, the center control plays no part in the convergence
analysis to Mnρ . Thus we see that the dynamics for xM are
perturbed by the term up due to formation error. However
since kx − xM k is decreasing exponentially to 0, we have
that up is bounded and also decreasing exponentially to 0.
Thus, there exists some time t? following which the conditions
given in Lemmas IX.4 and IX.5 are sufficient to ensure that
xM → Mnρ .
Remark IX.8 (Time Lag Parameter Bound). Note that the
bound for the time lag parameter used in the proof above
(derived in the proofs for Lemmas IX.4 and IX.5) is less
restrictive than the bound given in the statement of Theorem
V.2 in the main body of the text. Thus, satisfaction of the
constraint in (29) is sufficient for the convergence results to
hold. We present the more restrictive condition in the main
body of the text to simplify the exposition.
Finally, we present the decoupled framework for controlling
the formation size of a Johnson polyhedron by decoupling Q1
from the set Q. Suppose the tuple (n1 , Rη1 ) describes the
plane of the first polygon Q1 and let (n2 , Rη2 ) define the
plane of an adjacent polygon Q2 with robots j and j + 1
forming the shared edge. Define the variables νf , νb and kr as
follows:
 
0
kr = (xj+1 − xj ) · RηT2 0
(60)
1
(
1 if kr > 0,
νf =
(61)
0 else.
(
1 if kr < 0,
νb =
(62)
0 else.
Then, the control law for robot i in Q1 is given by:
(1)

ui =

N1
X

 (1) (1)
(1)
(1)
(1) 
(1)
(1)T
Rms (xi+m − xi ) + Rm
(xi−m − xi )
km
s

m=1


T
(1)
(1)
(1)
(1) 
+ |kr | νf Rr(1)
(xi+1 − xi ) + νb Rr(1)
(xi−1 − xi ) ,
s
s
(63)
where Rms = RηT1 Rmπ/|V1 |+αs (x(1) ) Rη1 , and Rrs =
RηT1 Rπ/|V1 | Rη1 . Thus, the control for all robots in Q1 is
augmented by a “one-sided” cyclic controller (that is flow(1)
invariant with respect to the polygon M|V1 | ) that aligns Q1
with its neighboring face Q2 . Note that only one of νf or νb
will be equal to 1 if Q1 does not possess the correct in-plane
rotation.

20

(1)

The proof of convergence for x(1) to the subspace M|V1 |
is a trivial extension of Lemma IX.4 by noting that the net
control for Q1 may be written as:

where T is the desired net thrust, and [adx , ady , adz ] are
the components of ad . Given a desired yaw angle ψ, we can
solve the above system of equations for the required thrust
(T ), roll (φ) and pitch (θ). The roll, pitch and yaw angles are
(1)
ẋ(1) = u(1) = − L(1)
, τ )x(1)
then converted into a desired attitude quaternion qd which is
η (x

− |kr (x(1) )| νf L1 ⊗ Rrs + νb LT1 ⊗ RrTs x(1) . fed into the attitude control loop.
The attitude control loop uses the desired attitude qd and the
A suitable auxiliary system for the above dynamics is given current measured attitude q to compute a desired correction
by:
quaternion ∆q given by:

∆q = qd q −1
(1)
ẏ (1) = −V̄ (1) L(1)
η (xp , τ )

T
= cos(β/2) sin(β/2)n̂x sin(β/2)n̂y sin(β/2)n̂z .


T
T
+|kr (x(1)
The desired moment is then given by the following
p )| νf L1 ⊗ Rrs + νb L1 ⊗ Rrs


proportional-derivative control law:
T
T
× V̄ (1) y (1) + Ū (1) Ū (1) x(1)
,
p
M = Kp β n̂ − Kd ω,
(1)
(1)
where V̄
is the orthonormal counterpart of Ve without the where Kp and Kd are positive semi-definite gain matrices.
matrix E (1) . Application of Corollary II.3 then requires:
Finally, the desired thrust and moment are fed into the quadcopter motor controller to compute the desired motor speeds.

1
(1)  (1)T
(1)
(1)
V̄ (1) L(1)
 0, uniformly, A block diagram for this controller is shown below in Figure
η (xp , τ ) + |kr (xp )|L1η V̄
2
11.
∀xp , τ.
(1)

(1)

T

Note that V̄ (1) |kr (xp )|L1η V̄ (1)  0, ∀xp . Thus, the
above condition reduces to (55) in Lemma IX.2. The proof of
convergence for polygon size is unchanged since the rotational
controller introduced in (63) does not affect the expression for
ṗi given in (57).
The control law for the remaining robots in the polyhedron
is given as before by (20) and (21). Conditions for convergence
for the reduced PPS (defined by the set of polygons Q\Q1 ) can
be derived by eliminating the first block in the net constraint
matrix given in (17). The robots in V1 converge to the polygon
Q1 with the desired size and an orientation consistent with the
reduced PPS (as dictated by Q2 ). In turn, all polygons in the
reduced PPS will then possess the desired size ρ as dictated
by the polygon Q1 .
X. Q UADCOPTER V ELOCITY T RACKING C ONTROL
In this section we present the quadcopter velocity tracking
controller used for the numerical simulations presented in the
main body of the paper. For notational convenience we drop
the i subscript notation referring to the ith quadcopter. In
this section we assume that all quantities refer to a single
quadcopter.
Let vd denote the desired velocity for a quadcopter, as
dictated by the cyclic control laws, and let v be the current
velocity. A simple proportional-integral-derivative controller
converts the velocity error ev given by vd − v into a desired acceleration ad . To achieve this desired acceleration, a
reference attitude quaternion qd must also be generated since
the quadcopter only posseses thrust capability in one direction
(opposite its body z axis). To do this, we equate the desired
acceleration to the actual acceleration dynamics given in (38)
and linearize the resulting equation assuming small pitch and
roll angles. This yields the following system of equations:

 
 
a
−T sin(ψ) −T cos(ψ) φ
= mq dx
ady
T cos(ψ) −T sin(ψ) θ
−T + mq g = madz ,

(64a)
(64b)

Fig. 11: Quadcopter Velocity Tracking Controller.

