1

Bio-Inspired Local Information-Based Control for
Probabilistic Swarm Distribution Guidance

arXiv:1711.06869v1 [cs.MA] 18 Nov 2017

Inmo Jang, Hyo-Sang Shin, and Antonios Tsourdos

Abstract—This paper addresses a task allocation problem for
a large-scale robotic swarm, namely swarm distribution guidance
problem. Unlike most of the existing frameworks handling this
problem, the proposed framework suggests utilising local information available to generate its time-varying stochastic policies.
As each agent requires only local consistency on information
with neighbouring agents, rather than the global consistency, the
proposed framework offers various advantages, e.g., a shorter
timescale for using new information and potential to incorporate
an asynchronous decision-making process. We perform theoretical analysis on the properties of the proposed framework. From
the analysis, it is proved that the framework can guarantee the
convergence to the desired density distribution even using local
information while maintaining advantages of global-informationbased approaches. The design requirements for these advantages
are explicitly listed in this paper. This paper also provides specific
examples of how to implement the framework developed. The
results of numerical experiments confirm the effectiveness and
comparability of the proposed framework, compared with the
global-information-based framework.
Index Terms—Swarm robotics, Distributed robot systems,
Networked robots, Markov chains.

I. I NTRODUCTION
HIS paper addresses a task allocation problem for a largescale multiple-robot system, called a robotic swarm.
Robotic swarms have attracted lots of attention because they
are regarded as promising solutions to handle complicated
missions that other systems may not be able to manage [1], [2].
Agents in a swarm are assumed to be homogeneous because
the swarm is usually realised through mass production [3]. In
this context, the task allocation problem can be reduced to a
problem of how to distribute a swarm of agents into given
tasks (or bins), satisfying the desired population fraction (or
swarm density) for each task. This problem is known as the
swarm distribution guidance problem [4]–[6].
For a large number of agents, probabilistic approaches based
on Markov chains [4]–[13] or differential equations [14]–[18]
have been widely utilised. Since these approaches focus not on
individual agents but instead on the ensemble dynamics, they
are also called Eulerian [11]–[13] or macroscopic frameworks
[18], [19]. In these approaches, swarm densities for each bin
are represented as system states, and a state-transition matrix
describes stochastic decision policies, i.e., the probabilities
that agents in a bin switch to another. Individual agents in
the swarm make decisions based on these policies, but in a
random, independent, and memoryless manner.

T

Inmo Jang, Hyo-Sang Shin, and Antonios Tsourdos are with the Centre for
Autonomous and Cyber-Physical Systems, Cranfield University, MK43 0AL,
United Kingdom (e-mail: inmo.jang@cranfield.ac.uk; h.shin@cranfield.ac.uk;
a.tsourdos@cranfield.ac.uk).

Initially, open-loop-type frameworks have been proposed
[4]–[7], [16]–[18]. Agents under these frameworks are controlled by time-invariant stochastic decision policies. The policies, which make a swarm converge to a desired distribution,
are pre-determined by a central controller and broadcasted
to each agent before executing the mission. Communication
between agents is hardly required during the mission, so
that it can reduce communication complexity under these
frameworks. However, the agents only have to follow the given
policies without incorporating any feedbacks, and thus there
still remain some agents who unnecessarily and continuously
switch bins even after the swarm reaches the desired distribution. This gives rise to a trade-off between convergence rate
and long-term system efficiency [17].
There have been also some other works, called closed-looptype frameworks [9]–[15]. This type of frameworks allows
agents to adaptively construct their own stochastic decision
policies at the expense of sensing the concurrent swarm
status through interactions with other agents. Based on such
information, agents can synthesise time-inhomogeneous transition matrices to achieve certain objectives and requirements:
for example, maximising convergence rates [10], minimising
travelling costs [13], and temporarily adjusting given policies
when bins are more overpopulated or underpopulated than
certain levels [14], [15]. In particular, Bandyopadhyay et. al.
[13] recently proposed a closed-loop-type algorithm that exhibits faster convergence as well as less undesirable transition
behaviours, compared with an open-loop-type algorithm. This
algorithm is expected to mitigate the trade-off raised in openloop-type frameworks.
To the best of our knowledge, most of the existing closedloop-type algorithms are based on Global Information Consistency Assumption (GICA) [20]. GICA implies that necessary
information is required to be consistently known by entire
agents. We refer to such information as global, because
achieving information consistency needs agents to somehow
interact with all the others through a multi-hop fashion and
thus it “happens on a global communication timescale” [20].
This paper proposes a framework that requires Local Information Consistency Assumption (LICA) [20]. Unlike GICAbased algorithms, the proposed framework require only local
consistency on information with neighbouring agents, not
the global consistency. LICA can provide various alternative
advantages to the proposed framework, compared with GICA.
Firstly, it “provides a much shorter timescale for using new
information because agents are not required to ensure that this
information has propagated to the entire team before using it”
[20]. Secondly, LICA enables a foundation on which an asyn-

2

chronous decentralised decision-making process can be developed. Note that the timescales for achieving the information
consistency between the agents can be different depending on
their local circumstances. Considering any possibly-extrinsic
heterogeneity of agents (e.g., different sensing frequency due
to local communication delays), an asynchronous algorithm is
regarded as more realistic in coordinating a robotic swarm,
so increasing its system efficiency [21]–[23]. Finally, LICA
makes the proposed approach additionally robust against dynamical changes in bins and those in agents. Given that
inclusions or exclusions of bins are perceived by neighbouring agents, the proposed approach works well even without
requiring other far-away agents to know the changes.
The LICA-based framework developed in this paper utilises
local information as its feedback gains, which is motivated
from the recent GICA-based work in [13]. This framework
is inspired by the mechanism of decision-making in a fish
swarm, in which each of them adjusts its individual behaviour
based on those of neighbours [24]–[27]. Similarly, each agent
in the framework developed uses its local status, i.e. the current
density of its associated bin relative to those of its neighbour
bins, to generate its time-varying stochastic decision policies.
The agent is not required to know any global information,
and hence the aforementioned advantages of LICA can be
exploited.
We prove that, even using local information, the proposed
framework asymptotically converges to a desired swarm distribution and it retains the advantages of existing closedloop-type approaches. This paper explicitly presents the design requirements for a time-inhomogeneous Markov chain
to achieve these desired features. It is thus expected that
the user can utilise the requirements in designing their own
algorithm. In addition, three specific examples are provided
to demonstrate how to implement the proposed framework:
1) minimising travelling cost; 2) maximising convergence rate
under upper flux bounds; and 3) generation of quorum-based
policies (similar to [14], [15]).
The rest of this paper are organised as follows. Section II introduces the desired features of a swarm distribution guidance
framework along with relevant definitions and notations. Section III proposes our framework with its design requirements,
the biological inspiration, and an analysis regarding whether
the desired features are satisfied. We provide examples of how
to implement the framework for specific problems in Section
IV, and an asynchronous implementation in Section V. The
results of numerical experiments are shown in Section VI,
followed by concluding remarks in Section VII.

Notations
∅, 0, I and 1 denote the empty set, the zero matrix of
appropriate sizes, the identity matrix of appropriate sizes, and
a row vector with all elements are equal to one, respectively.
v ∈ Pn is a stochastic (row) vector such that v ≥ 0 and
v ·1> = 1. v[i] indicates the i-th element of vector v. Prob(E)
denotes the probability that event E will happen.

TABLE I
N OMENCLATURE
Symbol

Description

Bi
A
nA
k
nk [i]
ajk
xjk
Mkj
Θ
µ?k
Ak
Ck
Nk (i)

The i-th bin amongst a set of nbin bins (Definition 1);
A set of agents (Definition 1);
The number of total agents at time instant k (Definition 1);
The number of agents at the i-th bin (Eqn. (4));
The j-th agent’s state indicator vector (Definition 2);
Stochastic state vector of the j-th agent (Definition 4);
Stochastic decision policy of the j-th agent (Definition 4);
Desired swarm distribution (Definition 5);
Current (global) swarm distribution (Definition 3);
Physical motion constraint matrix (Definition 6);
Communicational connectivity matrix (Definition 7);
A set of (communicationally-connected) neighbour bins of
the i-th bin (Definition 7);
A set of agents in Nk (i);
Current local swarm density at the i-th bin (Eqn. (4));
Estimate of µ̄?k [i] by the j-th agent;
Locally-desired swarm density at the i-th bin (Eqn. (5));
Primary guidance matrix (Eqn. (10));
Secondary guidance matrix (Eqn. (10));
Primary local-feedback gain (e.g., Eqn. (6));
Secondary local-feedback gain (Eqn. (9));

ANk (i)
µ̄?k [i]
µ̄jk [i]
Θ̄[i]
Pkj
Skj
ξ̄kj [i]
Gjk [i]

II. P RELIMINARIES
A. Definitions
This section presents necessary definitions and assumptions
for our proposed framework, which will be shown in Section
III. Since most of them are embraced from the recent existing
literature [10], [13], we here briefly provide their essential
meanings.
Definition 1 (Agents and Bins). A set of agents A are
supposed to be distributed over a prescribed region in a state
space B. The entire space is partitioned into nbin disjoint bins
bin
(subspaces) such that B = ∪ni=1
Bi and Bi ∩ Bj = ∅, ∀i 6= j.
We also regard B = {B1 , ..., Bnbin } as the set of all the bins.
Each bin Bi represents a predefined range of an agent’s state,
e.g., position. The number of the entire agents is time-varying,
and its value at time instant k is denoted by nA
k = |A|. Note
that we do not assume that the agents keep track of nA
k.
Definition 2 (Agent’s state). Let ajk ∈ {0, 1}nbin be the state
indicator vector of agent j ∈ A at time instant k. If the agent’s
state belongs to bin Bi , then ajk [i] = 1, otherwise 0.
Definition 3 (Current (global) swarm distribution). The current (global) swarm distribution µ?k ∈ Pnbin is a row-stochastic
vector such that each element µ?k [i] is the population fraction
(swarm density) of A in bin Bi at time instant k:
µ?k :=

1 X j
ak .
|A|

(1)

∀j∈A

Definition 4 (Agent’s stochastic state and decision policy).
Agent j’s stochastic state is a row-stochastic vector xjk ∈

3

Pnbin in which each element xjk [i] gives the probability that
the agent’s state belongs to bin Bi at time instant k:
xjk [i] := Prob(ajk [i] = 1).

(2)

The probability that agent j in bin Bi at time instant k will
transition to bin Bl before the next time instant is called its
stochastic decision policy, denoted as:
Mkj [i, l] := Prob(ajk+1 [l] = 1|ajk [i] = 1).

(3)

Note that Mkj ∈ Pnbin ×nbin is a row-stochastic matrix such
that Mkj ≥ 0 and Mkj · 1> = 1> , and will be referred as
Markov matrix.
Definition 5 (Desired swarm distribution). The desired swarm
distribution Θ ∈ Pnbin is a row-stochastic vector such that
each element Θ[i] indicates the desired swarm density for bin
Bi .
Assumption 1. For ease of description for this paper, we
assume that Θ[i] > 0, ∀i ∈ {1, ..., nbin }. Obviously, in practice, there may exist some bins whose desired swarm densities
are zero. These bins can be accommodated by adopting any
subroutines ensuring that all agents eventually move to and
remain in any of the positive-desired-density bins, for example,
an escaping algorithm in [13, Section III.C].
Assumption 2 (The number of agents [6], [10], [13], [15],
[17]). It is assumed that nA
k  nbin so that the time evolution of the swarm distribution is governed by the stochastic
decision policy in Equation (3). Although the finite cardinality
of the agents normally cause a residual convergence error, a
lower bound on nA
k that probabilistically guarantees a desired
convergence error is analysed in [13, Theorem 6] by exploiting
Chebyshev’s equality. Note that this theorem is generally
appliable and thus is also valid for our work.
Definition 6 (Physical motion constraint [6], [10], [13]).
Motion constraints of agents are denoted by the matrix Ak ∈
{0, 1}nbin ×nbin , where Ak [i, l] = 1 if agents in bin Bi at
time instant k are allowed to transition to bin Bl by the next
time instant; Ak [i, l] = 0, otherwise. It is assumed that Ak
is symmetric and irreducible (i.e., strongly-connected); and
Ak [i, i] = 1 for all agents, bins, and time instants.
Definition 7 (Communicationally-connected). Bins Bi and Bl
are said to be communicationally-connected, if there exists at
least one agent in bin Bi who can directly communicate with
some agents in bin Bl , and vice versa. This communicational
connectivity over all the bins at time instant k is defined by the
matrix Ck ∈ {0, 1}nbin ×nbin , where Ck [i, l] = 1 indicates that
bins Bi and Bl are communicationally-connected. Note that
Ck is symmetric and all its diagonal entries are set to be one.
For each bin Bi , we define the set of its (communicationallyconnected) neighbour bins as Nk (i) = {∀Bl ∈ B | Ck [i, l] =
1}. The set of agents in any of bins in Nk (i) is denoted by
ANk (i) = {∀j ∈ A | ajk [l] = 1, ∀l : Bl ∈ Nk (i)}.
Assumption 3 (Communicational connectivity over bins). The
physical motion constraint of a robot is, in general, more
stringent than its communicational constraint. From this, it can

be assumed that if the transition of agents between bin Bi and
Bl is allowed within a unit time instant, then the both bins
are communicationally-connected, i.e., if Ak [i, l] = 1 then
Ck [i, l] = 1. Note that we set Ck [i, l] = 0 if Ak [i, l] = 0.
This implies that the matrix Ck is irreducible, as is Ak .
The communication network over the agents is assumed to
be strongly-connected [10], [13]. Using distributed consensus
algorithms [10], [11], [28], each agent can access necessary
local information in its neighbour bins.
Assumption 4 (Pre-known Information [13]). The desired
swarm distribution Θ, the motion constraint matrix Ak (also
Ck ), and other pre-determined values such as variables regarding objective functions and user-design parameters (which will
be introduced later) are known by all the agents before they
begin a mission.
Assumption 5 (Agent’s capability [10], [13]). Each agent can
determine the bin to which it belongs, and know the locations
of neighbour bins so that it can navigate toward any of these
bins. The agent is capable of collision avoidance behaviours
against other agents or obstacles.
B. Problem Statement
The objective of the swarm distribution guidance problem
considered in this paper is to distribute a set of agents A over
a set of bins B by the Markov matrix Mkj in a manner that
holds the following desired features:
Desired Feature 1. The swarm distribution µ?k asymptotically
converges to the desired swarm distribution Θ as time instant
k goes to infinity.
Desired Feature 2. Transitions of the agents between the bins
are controlled in a way that Mkj becomes close to I as µ?k
converges to Θ. This implies that the agents are settled down
after Θ is achieved, and thus unnecessary transitions can be
reduced. Moreover, the agents identify and compensate any
partial loss or failure of the swarm distribution.
Desired Feature 3. For each agent in bin Bi , the information required for generating time-varying stochastic decision policies is not global information (e.g., µ?k ) but locally
available information within ANk (i) . Thereby, the resultant
time-inhomogeneous Markov process is based on LICA, and
has benefits such as a shorter timescale for obtaining new
information (than GICA), the potential for an asynchronous
process, etc.
Remark 1. One of our main contributions is to provide
Desired Feature 3 as well as to retain Desired Features 1 and 2
by additionally adopting Assumption 3, which can be elicited
from other assumptions in the existing literature.
III. A C LOSED - LOOP - TYPE F RAMEWORK USING L OCAL
I NFORMATION
This section proposes a LICA-based framework for the
swarm distribution guidance problem. The framework is different from the recent closed-loop-type algorithms in [10], [13]
in the sense that they utilise the global information (e.g., the

4

current swarm distribution in Equation (1)) for constructing a
time-inhomogeneous Markov matrix, whereas ours uses the
local information in Equation (7). We present, in spite of
using such relatively insufficient information, how the desired
features described in the previous section can be achieved
in the proposed framework. Before that, we introduce the
biological idea, which is about decision-making mechanisms
of a fish swarm, that inspires this framework to particularly
attain Desired Feature 3. In addition, we explicitly provide
the design requirements for a Markov matrix in order for
prospective users to easily incorporate their own specific
objectives into this framework.

1

2

3

4

…

10

11

12

13

14

…

20

21

22

23

24

…

30

31

32

33

34

…

40

...

…

…

…

…

…

91

92

93

94

… 100

Fig. 1. Examples of how to calculate µ̄?k [i]: for bin B23 , µ̄?k [23] =
nk [23]/(nk [13] + nk [22] + nk [23] + nk [24] + nk [33]). In the proposed
framework, agents in the bin only need to obtain the local information
from other agents in its neighbour bins (shaded). Note that each square
indicates each bin, and the red arrow between two bins Bi and Bl means
that Ak [i, l] = 1.

A. The Biological Inspiration
For a swarm of fishes, it has commonly been assumed
that their crowdedness limits their perception ranges over
other members, and their cardinality restricts the capacity for
individual recognition [25]. How fishes end up with collective
behaviours is different from the ways of other social species
such as bees and ants, which are known to use recruitment
signals for the guidance of the entire swarm [29], [30]. Thus,
in biology domain, a question naturally has arisen about
the mechanism of fishes’ decision-making in an environment
where local information is only available and information
transfer between members does not explicitly happen [24]–
[27], [31], [32].
It has been experimentally shown that fishes’ swimming
activities vary depending on their perceivable neighbours.
According to [31], fishes have the tendency to maintain their
statuses (e.g., position, speed, and heading angle) relative to
those of other nearby fishes, which results in their organised
formation structures. In addition, it is presented in [32] that
spatial density of fishes has influences on both the minimum
distances between them and the primary orientation of the fish
school.
Based on this knowledge, the works in [24]–[27] suggest
individual-based models to further understand the collective
behavioural mechanisms of fishes: for example, their repelling,
attracting, and orientating behaviours [24], [26]; how the density of informed fishes affects the elongation of the formation
structure [25]; and group-size choices [27]. The common and
fundamental characteristic of these models is that every agent
maintains or adjusts its personal status with consideration of
those of other individuals within its limited perception range.
As inspired by the understanding of fishes, we believe
that there must be an enhanced swarm distribution guidance
approach in which each agent only needs to keep its relative
status by using local information available from its nearby
neighbours. In this approach, a global information is not
necessary to be known by agents, and thereby the corresponding requirement of extensive information sharing over all the
agents can be alleviated.
B. Fundamental Idea of the Proposed Approach
Suppose that each agent in bin Bi is required to keep its
local status µ̄?k [i], which we referred to as the current local
swarm density at bin Bi , at the value of the corresponding

locally-desired swarm density Θ̄[i]. They are respectively
defined as follows:
nk [i]
µ̄?k [i] := P
,
(4)
∀l:Bl ∈Nk (i) nk [l]
where nk [i] is the number of agents such that ajk [i] = 1; and
Θ̄[i] := P

Θ[i]
∀l:Bl ∈Nk (i)

Θ[l]

.

(5)

We use the term µ̄jk [i] as an estimate of µ̄?k [i] by agent j, which
can be obtained through a distributed information consensus
algorithm [10], [11], [28].
The fundamental idea of the proposed approach is to make
each agent j in bin Bi :
(i) only need to estimate the difference of µ̄?k [i] and
Θ̄[i], which are both locally-available information within
Nk (i); and
(ii) more reluctant to deviate from the current bin as the
difference becomes smaller (i.e., Mkj [i, i] → 1 and
Mkj [i, l] → 0 ∀l as µ̄jk [i] → Θ̄[i]).
Our proposed framework utilises the difference between
µ̄jk [i] and Θ̄[i] as a local-information-based feedback gain,
denoted by ξ¯kj [i], which is a scalar in (0, 1] that monotonically
decreases as µ̄jk [i] converges to Θ̄[i]. For instance, this paper
uses

|Θ̄[i]−µ̄jk [i]| α

if |Θ̄[i] − µ̄jk [i]| ≤ Θ̄[i]

( Θ̄[i] )
|Θ̄[i]−µ̄jk [i]| α
(6)
ξ¯kj [i] := ξ
if
(
) < ξ
Θ̄[i]



1
otherwise
where α > 0 and ξ > 0 are design parameters. We call
this gain primary local-feedback gain because it is utilised to
control the primary guidance matrix Pkj (shown in the next
subsection).
Remark 2. Equation (4) is equivalent to the i-th element of
the following vector:
X
1
µ̄k (i) =
ajk .
(7)
|ANk (i) |
∀j∈ANk (i)

µ̄?k [l]

Namely,
= µ̄k (i)[l] if l = i. Here, we intentionally
introduce Equation (7) for ease of comparison with the information required for feedback gains in the existing literature

5

(e.g., Equation (1)). From this, it is implied that, in order for
each agent in bin Bi to estimate µ̄k (i)[i] (i.e., the current
local swarm density µ̄?k [i]), the set of other agents whose
information is necessary is restricted within ANk (i) . That is,
each agent needs to have neither a large perception radius
nor an extensive information consensus process over the entire
agents.

This subsection presents our closed-loop-type framework
based on locally-available information feedbacks. The basic
form of the stochastic decision policy for agent j in bin Bi is
such that
(
(1 − ωkj [i])Pkj [i, l] + ωkj [i]Skj [i, l] if l = i
j
Mk [i, l] :=
(1 − ωkj [i])Pkj [i, l]
∀l 6= i.
(8)
Here, ωkj [i] ∈ [0, 1) is the weighting factor to have different
weights on the agent’s primary decision policy Pkj [i, l] ∈ P
and secondary decision policy Skj [i, l] ∈ P. It is defined as
(9)

where τ j is a design parameter; and Gjk [i] ∈ [0, 1] is secondary
local-feedback gain, which is based on the difference between
µ̄jk [i] and Θ̄[i]. Note that ωkj [i] is mainly affected by Gjk [i],
while diminishing as time instant k goes to infinity.
Equation (8) can be represented in matrix form as
Mkj = (I − Wkj )Pkj + Wkj Skj ,
Pkj

Requirement 4. Pkj is irreducible such that
Pkj [i, l] > 0

if Ck [i, l] = 1.

Pkj [i, l]

otherwise.

=0

(R4)

Note that Ck is already assumed to be irreducible in Assumption 3.

C. A LICA-based Closed-loop-type Framework

ωkj [i] := exp(−τ j k) · Gjk [i]

With consideration of (R1), this requirement can be fulfilled
by Θ[i]Pkj [i, l] = Θ[l]Pkj [l, i], ∀i. A Markov process satisfying
this property is said to be reversible.

(10)

Skj

where
∈ Pnbin ×nbin and
∈ Pnbin ×nbin are rowstochastic matrices, called primary guidance matrix and secondary guidance matrix, respectively. Wkj ∈ Rnbin ×nbin is a
diagonal matrix such that diag(Wkj ) = (ωkj [1], ..., ωkj [nbin ]).
The stochastic state vector of agent j is governed by the
Markov process:
xjk+1 = xjk Mkj .
(11)
For now, we claim that, in order for this Markov system to
achieve Desired Features 1-3, Pkj must satisfy the following
requirements.
Requirement 1. Pkj is a matrix with row sums equal to one,
i.e.,
n
bin
X
Pkj [i, l] = 1, ∀i.
(R1)

Requirement 5. Pkj becomes close to I as µ̄jk converges to
Θ̄, i.e.,
Pkj [i, i] → 1 as µ̄jk [i] → Θ̄[i] (or ξ¯kj [i] → 0), ∀i.
Pkj ,

Skj ,

(R5)
ξ¯kj [i]

Depending on the objectives of a user,
and
Gjk [i] can be designed differently under given specific constraints. As long as Pkj holds (R1)-(R5) for all time instant
k and all agent j ∈ A, the aforementioned desired features
are achieved. Note that (R1)-(R4) are associated with Desired
Feature 1, whereas (R5) is with Desired Feature 2. The detailed
analysis will be described in the next subsection.
Every agent executes the following algorithm at every time
instant. The detail regarding Line 6-8 will be presented in
Section IV, which shows examples of how to implement this
framework.
Algorithm 1 Decision making of agent j at time instant k
// Obtain the local information
1: Identify the current bin Bi ;
2: Identify neighbour bins Nk (i) (and Ck [i, l] ∀l);
3: Compute Θ̄[i] using (5);
j
4: Obtain µ̄k [i];
// Generate the stochastic decision policy
j
5: Compute ξ¯k [i] (using (6));
j
6: Compute Pk [i, l] ∀l;
j
7: Compute Sk [i, l] ∀l;
j
8: Compute Gk [i];
j
9: Compute ωk [i] using (9);
j
10: Compute Mk [i, l] ∀l using (8);
// Individually behave based on the policy
11: Generate a random number z ∈ unif[0, 1];
12: Select bin Bq such that
Pq−1 j
Pq
j
l=1 Mk [i, l];
l=1 Mk [i, l] ≤ z <
13: Move to the selected bin;

l=1

Pkj

In fact,
needs to be row-stochastic, for which it should
further hold that Pkj [i, l] ≥ 0, ∀i, l. Note that this constraint is
implied by (R4), which will be introduced later.
Requirement 2. All diagonal elements are positive, i.e.,
Pkj [i, i] > 0, ∀i.

(R2)

Requirement 3. The stationary distribution of Pkj is the
desired swarm distribution Θ, i.e.,
n
bin
X
i=1

Θ[i]Pkj [i, l] = Θ[l], ∀l.

(R3)

D. Analysis
We first show that the Markov process in Equation (11)
holds Desired Feature 1 under the assumption that Pkj satisfies the requirements (R1)-(R4) for each time instant. The
stochastic state of agent j at time instant k ≥ k0 , governed by
the Markov process from an arbitrary initial state xjk0 , can be
written as:
j
xjk = xjk0 Ukj0 ,k := xjk0 Mkj0 Mkj0 +1 · · · Mk−1
.

(12)

For ease of analysis, we assume that every agent j knows any
necessary information correctly, i.e., µ̄jk [i] = µ̄?k [i].

6

Theorem 1. Provided that the requirements (R1)-(R4) are satisfied for all time instants k ≥ k0 , it holds that limk→∞ xjk =
Θ pointwise for all agents, irrespective of the initial condition.
Proof. This claim can be proved by following similar steps in
proving [13, Theorem 4]. The claim is true if limk→∞ xjk =
xjk0 · limk→∞ Ukj0 ,k = xjk0 · 1> Θ = Θ. In order for that,
the matrix product Ukj0 ,k should (i) be strongly ergodic and
(ii) have Θ as its unique limit vector, i.e., limk→∞ Ukj0 ,k =
1> Θ. We will show that the two conditions are valid under
the assumption that (R1)-(R4) are satisfied.
Lemma 5 in Appendix describes the characteristics of Mkj
and Ukj0 ,k , which will be used for the rest of this proof. From
this lemma, (a) Ukj0 ,k is primitive (thus, regular); (b) there
exists a positive lower bound γ for Mkj , ∀k; and (c) Mkj is
asymptotically homogeneous. Then, from [33, Theorem 4.15,
p.150] it follows that Ukj0 ,k is strongly ergodic, which fulfils
the condition (i).
Let ek ∈ Pnbin be the unique stationary distribution vector
corresponding to Mkj (i.e., ek Mkj = ek ). Due to the prior
condition (b) and the fact that (d) Mkj is irreducible for
∀k ≥ k0 , it follows from [33, Theorem 4.12, p.149] that
the asymptotical homogeneity of Mkj with respect to Θ (i.e.,
limk→∞ ΘMkj = Θ) is equivalent to limk→∞ ek = e and
Θ = e, where e is a limit vector. According to [33, Corollary,
p.150], under the prior conditions (b) and (d), if Ukj0 ,k is
strongly ergodic with its unique limit vector v, then v = e.
Hence, it turns out that the unique limit vector of Ukj0 ,k is Θ
(i.e, limk→∞ Ukj0 ,k = 1> Θ). Thereby, the condition (ii) is also
fulfilled.

1

2

1

2

4

3

4

3

(a) Tree-type

(b) Strongly-connected

Fig. 2. Examples of simple bin topologies to help Lemma 1 & 2: (a) treetype; (b) strongly-connected. The red line in (b) indicates a newly-added route
between bin B1 and B4 based on the topology in (a).

Proof. The matrix B ∈ Rnbin ×nbin can be linearly decomposed into ne of the same-sized matrices B(i,j) , where ne is
the number of edges in the underlying graph of Ck . Here,
B(i,j) ∈ Rnbin ×nbin is a matrix such that B(i,j) [i, i] =
−Θ[j]/Θ[i] and B(i,j) [j, j] = −Θ[i]/Θ[j]; B(i,j) [i, j] =
B(i,j) [j, i] = 1; and all the other entries are zero. For example,
consider that four bins are given and connected as shown in
Figure 2(a). Clearly, B = B(1,2) + B(2,3) + B(2,4) , where

 Θ[2]
1
0
0
− Θ[1]


 1
− Θ[1]+Θ[3]+Θ[4]
1
1 
Θ[2]
,

B=

1
− Θ[2]
0

 0
Θ[3]
Θ[2]
0
1
0
− Θ[4]


− Θ[2]
Θ[1]

B(2,3)

0
0

=
0
0


B(2,4)

0
0

=
0
0

0

Θ[1]
− Θ[2]
0
0



B(1,2) =  1
 0
0


Theorem 1 implies that the stochastic state of any agent
eventually converges to the desired swarm distribution, regardless of Skj , Gjk [i] and (R5). In other words, even if (R5) is not
satisfied, the Markov system can converge to Θ. However, the
system induces unnecessary transitions of agents even after
being close enough to the desired swarm distribution, which
means that Desired Feature 2 does not hold.
For now, we present that Desired Feature 2 can be obtained
by (R5) and Theorem 2, which will be described later. Suppose
that, for every bin Bi , µ̄?k [i] converges to and eventually
reaches Θ̄[i] at some time instant k. The following shows
that at this moment it also holds that µ?k reaches Θ. From
Equations (4)-(5) and the supposition
of µ̄?k [i] = Θ̄[i] ∀i, it
P
follows that 1/Θ̄[i] · nk [i] = ∀j:Bj ∈Nk (i) nk [j] ∀i. This can
be rearranged as:


0

0 0
,
0 0
0 0

1

0
− Θ[3]
Θ[2]

0
1

1
0

− Θ[2]
Θ[3]
0

0

0
0
0
0

Θ[4]
− Θ[2]
0
1


0
0

,
0
0
0
1
0

− Θ[2]
Θ[4]




.


where X ∈ Rnbin ×nbin is a diagonal matrix such that
diag(X) = (1/Θ̄[1], 1/Θ̄[2], ..., 1/Θ̄[nbin ]); Ck is the communicational connectivity matrix (in Definition 7); and nk ∈
Rnbin is a row vector such that the i-th element indicates nk [i],
i.e., the number of agents in bin Bi at time instant k.

It is trivial that the rank of every B(i,j) is one, and the matrix
has only one linearly independent column vector, denoted by
v(i,j) . Without loss of generality, we consider v(i,j) ∈ Rnbin
1
as a column vector such that the i-th entry is − Θ[i]
, the j-th
1
entry is Θ[j] , and the others are zero: for an instance, v(1,2) =
1
1
[− Θ[1]
, Θ[2]
, 0, 0]> .
It is obvious that v(i,j) and v(k,l) are linearly independent
when the bin pairs {i, j} and {k, l} are different. This implies
that the number of linearly independent column vectors of B
is the same as that of edges in the topology. Hence, for a treetype topology of nbin bins, since there exist nbin − 1 edges,
the rank of the corresponding matrix B is nbin − 1.

Lemma 1. Given nbin bins communicationally-connected as
a tree-type topology, the rank of its corresponding matrix B
in Equation (13) is nbin − 1.

Lemma 2. Given a strongly-connected topology of bins, the
rank of its corresponding matrix B is not affected by adding
a new edge that directly connects any two existing bins.

nk · B := nk · (Ck − X) = 0

(13)

7

Proof. We will show that this claim is valid even when a
tree-type topology is given, as it is a sufficient condition for
strong-connectivity. Given the tree-type topology in Figure
2(a), suppose that bin B1 and B4 are newly connected. Then,
the new topology becomes as shown in Figure 2(b), and it has
new corresponding matrix Bnew , where Bnew = B + B(1,4) .
As explained in the proof of Lemma 1, the rank of B(1,4) is one
and it has only linearly independent vector v(1,4) . However,
this vector can be produced as a linear combination of the
existing v vectors of B (i.e., v(1,4) = v(1,2) + v(2,4) ). Thus,
the rank of Bnew retains that of B. Without loss of generality,
this implies that the rank of B of a given strongly-connected
topology is not affected by adding a new edge that directly
connects any two existing bins.
Thanks to Lemma 1 and 2, we end up with the following
corollary and theorem:
Corollary 1. Given nbin bins that are communicationally
strongly-connected, the rank of its corresponding B is nbin −1.
Theorem 2. Given nbin bins that are communicationally
strongly-connected, convergence of µ̄?k to Θ̄ is equivalent to
convergence of µ?k to Θ.
Proof. From Equation (5), it can be said that Θ·B = 0. When
µ̄?k [i] is assumed to converge to Θ̄[i] at some time instant k for
every bin Bi , Equation (13) is valid (i.e., nk · B = 0). Since
the nullity of B is one, due to Corollary 1, there is only one
linearly-independent row-vector a ∈ Rnbin such that a·B = 0.
Hence, it is obvious that nk =  · Θ, where  is an arbitrary
scalar value. This implies that µ?k [i] = nk [i]/nA
k = Θ[i], ∀i :
Bi ∈ B. Therefore, convergence of µ?k to Θ is equivalent to
convergence of µ̄?k to Θ̄.
From this theorem and (R5), Desired Feature 2 finally holds.
Corollary 2. If Pkj satisfies (R5), it can be said from Theorem
2 that Pkj becomes I as µ?k converges to Θ. And this is also
the case for the Markov process Mkj , which satisfies Desired
Feature 2.
In order for each agent j in bin Bi to generate the timevarying stochastic decision policy Mkj [i, l] in Equation (8), the
agent only needs to obtain its local information within ANk (i) .
Therefore, Desired Feature 3 is also achieved.
Remark 3 (Robustness against dynamic changes of agents
and those of bins). The proposed framework is robust with
against dynamic changes in the number of agents and bins.
Similarly to what is claimed in [13, Remark 8], as each agent
behaves based on its current bin location and local information
in a memoryless manner, Desired Features 1-3 in the proposed
framework won’t be affected by inclusion or exclusion of
agents in a swarm. Furthermore, provided changes on bins
are perceived by at least nearby agents in the corresponding
neighbour bins, robustness against those changes can be hold
in the proposed framework. This is because agents in bin
Bi utilise only local information such as Θ̄[i] and µ̄jk [i],
and are not required to know any information from other
far-away bins. Moreover, the proposed framework does not
need to recalculate Θ, reflecting such changes on bins, so

P
that
∀i Θ[i] = 1 because computing Θ̄[i] in (5) includes
normalisation of Θ.
IV. I MPLEMENTATION E XAMPLES
A. Example I: Minimising Travelling Expenses
This section provides examples on implementations of the
framework proposed. In particular, this subsection addresses
a problem of minimising travelling expenses of agents during
convergence to a desired swarm distribution.
This problem can be defined as: given a cost matrix
Ek ∈ Rnbin ×nbin in which each element Ek [i, l] represents
the travelling expense of an agent from bin Bi to Bl , find Pkj
such that
n
bin n
bin
X
X
min
Ek [i, l]Pkj [i, l]
(P1)
i=1 l=1

subject to (R1)-(R5) and
M Θ[l]f (ξ¯kj [i], ξ¯kj [l])f (Ek [i, l]) ≤ Pkj [i, l]
if Ck [i, l] = 1, ∀i 6= l

(14)

where M ∈ (0, 1] is a design parameter. f (ξ¯kj [i], ξ¯kj [l]) ∈ (0, 1]
is set by
f (ξ¯kj [i], ξ¯kj [l]) = max(ξ¯kj [i], ξ¯kj [l])
(15)
so that the value monotonically increases with regard to
increase of either ξ¯kj [i] or ξ¯kj [l] and diminishes as ξ¯kj [i] and
ξ¯kj [l] simultaneously reduces. This value controls the lower
bound of Pkj [i, l] in Equation (14). Θ[l] enables agents in bin
Bi to be distributed over its neighbour bins in proportion to the
desired swarm distribution. f (Ek [i, l]) ∈ (0, 1] is a scalar that
monotonically decreases as Ek [i, l] increases (see Equation
(29) for instances), encouraging agents in bin Bi to avoid
spending higher transition expenses. Note that we assume that
Ek is symmetric; Ek [i, l] > 0 if Ak [i, l] = 1; and its diagonal
entries are zero.
Corollary 3. The optimal matrix Pkj of the problem (P1) is
given by: ∀i, l ∈ {1, ..., nbin } and i 6= l,
(
M Θ[l]f (ξ¯kj [i], ξ¯kj [l])f (Ek [i, l]) if Ck [i, l] = 1
j
Pk [i, l] =
0
otherwise
(16)
and ∀i = l,
X j
Pkj [i, i] = 1 −
Pk [i, l].
(17)
∀l6=i

Proof. We can prove this by following the proof of [13,
Corollary 1]. Suppose that the problem is only subject to (R4)
and (14), without (R1)-(R3) and (R5). Then, the off-diagonal
elements of an optimal matrix should be their corresponding
lower bounds in (14) if Ck [i, l] = 1. The diagonal elements of
the matrix do not affect the objective function due to the fact
that Ek [i, i] = 0, ∀i. Accordingly, the matrix Pkj that holds
(16) and (17) is also an optimal matrix for the simplified
problem.
Let us now consider (R1)-(R3) and (R5). Since M ,
fP(ξ¯kj [i], ξ¯kj [l]) and f (Ek [i, l]) are upper-bounded by 1 and
j
∀l6=i Θ[l] < 1, Pk [i, i] in (17) is always positive for all i,
which fulfils (R2). It is also obvious that (R1) is satisfied by

8

Equation (17). From Equation (16), it holds that Θ[i]Pkj [i, l] =
Θ[l]Pkj [l, i], complying with (R3). Since (R1)-(R4) are satisfied, the Markov process is converging to a desired distribution
due to Theorem 1. Noting that f (ξ¯kj [i], ξ¯kj [l]) diminishes as
ξ¯kj gets close to zero (i.e., µ̄jk → Θ̄), (R5) is also fulfilled
by Equation (17). Hence, Pkj is the optimal solution for the
problem (P1).
For reducing unnecessary transitions of agents during this
process, it is favourable that agents in bin Bi such that µ̄jk [i] ≤
Θ̄[i] (i.e., underpopulated) do not deviate. To this end, we set
Skj = I and Gjk [i] as follows [13]:
Gjk [i] :=

exp(β(Θ̄[i] − µ̄jk [i]))
exp(β|Θ̄[i] − µ̄jk [i]|)

.

(18)

The gain value is depicted in Figure 3(a) with regard to β.
Remark
of P
Convergence Rate). Due to the fact
P 4 (Increase
j
that
P
[i,
l]
≤
∀l6=i k
∀l:Bl ∈Nk (i)\Bi Θ[l] from Equation
(16), the total outflux of agents from bin Bi becomes smaller as
the bin has fewer connections with other bins. This eventually
makes the convergence rate of the Markov process slower.
Adding an additional variable into Pkj [i, l] in (16) does not
affect the obtainment of Desired Features 1-3 as long as Pkj
satisfies (R1)-(R5). Thus, in order to enhance the convergence
rate under the requirements, one can add
Θ := min{ P

1

∀s:Bs ∈Nk (i)\Bi

Θ[s]

,P

1

∀s:Bs ∈Nk (l)\Bl

}

Θ[s]
(19)

into Pkj [i, l], as follows:

¯j
¯j

Θ M Θ[l]f (ξk [i], ξk [l])f (Ek [i, l])
j
Pk [i, l] =
if Ck [i, l] = 1


0
otherwise,
(20)
which can be substituted for Equation (16).
Algorithm 2 Subroutine of Algorithm 1 (Line 6–8) for P1
Compute Pkj [i, l] ∀l using (16) (or (20)) and (17);
j
j
2: Set Sk [i, i] = 1 and Sk [i, l] = 0, ∀l 6= i
j
3: Compute Gk [i] using (18);
1:

for mitigating the trade-off between convergence rate and longterm system efficiency.
For the sake of imposing upper flux bounds during the entire
process, we consider the following one-way flux constraint:
nk [i]Pkj [i, l] ≤ c(i,l) ,

∀i, ∀l 6= i.

(21)

This means that the number of agents moving from bin Bi
to Bl is upper-bounded by c(i,l) . The bound value is assumed
to be very small with consideration of mission environments
such as the number of agents, the number of bins, and their
topology. Otherwise, all the agents can be distributed over
the bins very soon so that the upper flux bounds become
meaningless, and thus the corresponding problem can be
trivial.
Regarding the convergence rate of a Markov chain, there
are respective analytical methods depending on whether it
is time-homogeneous or time-inhomogeneous. For a timehomogeneous Markov chain, if the matrix is irreducible, the
second largest eigenvalue of the matrix is used as an index
indicating its asymptotic convergence rate [34, p.389]. In
contrast, for a time-inhomogeneous Markov chain, coefficients
of ergodicity can be utilised as a substitute for the second
largest eigenvalue, which is not useful for this case [35].
Particularly, this paper uses the following proper coefficient
of ergodicity, amongst others:
Definition 8. (Coefficient of Ergodicity [33, pp. 136–137]).
Given a stochastic matrix M ∈ Pn×n , a (proper) coefficient
of ergodicity 0 ≤ τ (M) ≤ 1 can be defined as:
τ (M) := max max |M[i, s] − M[l, s]|.
∀s

∀i,∀l

(22)

A coefficient of ergodicity is said to be proper if τ (M) = 0
if and only if M = 1> · v, where v ∈ Pn is a row-stochastic
vector.
The convergence rate of a time-inhomogeneous Markov
chain Mk ∈ Pn×n , ∀k > 1 can be maximised by minimising
τ (Mk ) at each time instant k,Qthanks to [33, Theorem 4.8,
r
p.137]: τ (M1 M2 · · · Mr ) ≤ k=1 τ (Mk ). Hence, the objective of the specific problem considered in this subsection
can be defined as: find Pkj such that
min τ (Pkj )

(23)

subject to (R1)-(R5) and (21).
B. Example II: Maximising Convergence Rate within Upper
Flux Bounds
This subsection presents an example in which the specific
objective is to maximise the convergence rate under upper
bounds regarding transitions of agents between bins, denoted
by upper flux bounds. The bounds can be interpreted as safety
constraints in terms of collision avoidance and congestion:
higher congestions may induce higher collisions amongst
agents, which may bring unfavourable effects on system
performance. A similar problem is addressed by an open-looptype algorithm in [17], where transitions of agents are limited
only at a desired swarm distribution. This restriction is not for
considering the aforementioned safety constraints, but rather

Remark 5 (Advantages of the coefficient of ergodicity in (22)).
Other proper coefficients in [33, p. 137] such as
X
τ1 (M) = 1 − min
min (M[i, s], M[l, s])
i,l

∀s

or
τ2 (M) = 1 −

X
∀s

min (M[i, s]) .
∀i

may have the trivial case such that τ1 (Pkj ) = 1 (or τ2 (Pkj ) = 1)
for some time instant k, when they are applied to this problem.
This is because, given a strongly-connected topology Ck , there
may exist a pair of bins Bi and Bl such that Pkj [i, s] = 0
or Pkj [l, s] = 0, ∀s. To avoid this trivial case, the work in

9

[13] instead utilises τ1 ((Pkj )dCk ) as the proper coefficient of
ergodicity, where dCk denotes the diameter of the underlying
graph of Ck . However, this implies that agents in bin Bi are
required to additionally access the information from other bins
beside Nk (i), causing additional communicational costs. The
coefficient of ergodicity in (22) does not suffer this issue. Note
that τ (M) ≤ τ1 (M) ≤ τ2 (M) [33, p. 137].
Finding the optimal solution for the problem (23) is another
challenging issue, which can be called fastest mixing Markov
chain problem. Since the purpose of this section is to show an
example of how to implement our proposed framework, we
heuristically address this problem at this moment.
Suppose that matrix Pkj satisfying (R1)-(R5) is given,
and the topology of bins is not fully-connected. Since the
matrix is non-negative and there exists at least one zerovalue entry in each column, the coefficient of ergodicity
can be said as τ (Pkj ) = max∀i,∀s (Pkj [i, s]). Assuming that
max∀l6=i Pkj [i, l] ≤ 1/|Nk (i)|, which is generally true due
to the smallness of c(i,l) , it turns out that each diagonal
element of Pkj is the largest value in each row. Thus, we can
say that τ (Pkj ) = max∀i Pkj [i, i]. TheP
objective function of
this problem can be said as maxmin∀i ∀l6=i Pkj [i, l] because
minimising the maximum diagonal element of a stochastic
matrix is equivalent to maximising the minimum row-sum of
its off-diagonal elements.
We turn now to the constraints (R1)-(R5) and (21). In order
to comply with (R3), we initially set Pkj [i, l] = Θ[l]Qjk [i, l],
where Qjk is a symmetric matrix that we will design now.
The constraint (21), (R4), and the symmetricity of Qk are
integrated into the following constraint: ∀i, ∀l 6= i,
c(l,i)
c(i,l)
,
) ≥ Qjk [i, l] > 0 if Ck [i, l] = 1
min(
nk [i]Θ[l] nk [l]Θ[i]
Qjk [i, l] = 0 otherwise.
(24)
For (R2) and (R5), we set the diagonal entries of Pkj as
Pkj [i, i] ≥ 1 − ξ¯kj [i],

∀i.

This
Pnbincan bej rewritten, with consideration of (R1) (i.e.,
l=1 Θ[l]Qk [i, l] = 1, ∀i), as
X
Θ[l]Qjk [i, l] ≤ ξ¯kj [i], ∀i.
(25)
∀l6=i

Then, the reduced problem can be defined as: find Qjk such
that
X
maxmin∀i
Θ[l]Qjk [i, l]
(P2)
∀l6=i

subject to (24) and (25).
The algorithm for this problem is shown in Algorithm 3.
If we neglect (25), an optimal solution can be obtained by
making Qjk [i, l] equal to its upper bound of (24) (Line 2).
However, this solution may not hold (25). Thus, we lower the
entries of Qjk to satisfy (25), while keeping them symmetric
and as higher as possible (Line 3–9). In details, Line 3 (or
Line 6) ensures the constraint (25) for each bin Bi in a way
that, if this is not the case, obtains the necessary lowering
factor ¯0Q [i] (or ¯Q [i]). In order to keep Qjk as higher as

possible, we temporarily take 0Q [i, l] as the maximum value
of {¯
0Q [i], ¯0Q [l]} (Line 4). After curtailing Qjk [i, l] by applying
0
Q [i, l], we obtain the corresponding lowering factor again
(Line 5–6). The minimum value is taken for both maintaining
Qjk symmetric and satisfying (25) (Line 7). Then, the corresponding stochastic decision policy is generated based on the
resultant Qjk (Line 9–10). Note that we set Gjk [i] = 0 for all
time instants, all bins, and all agents, so Mkj = Pkj .
Algorithm 3 Subroutine of Algorithm 1 (Line 6) for P1
// Initialise Pkj
j
1: Pk [i, l] = 0, ∀l ∈ {1, 2, ..., nbin };
// Compute Qjk satisfying (24)
c(l,i)
c(i,l)
j
2: Qk [i, l] = min( n [i]Θ[l] , n [l]Θ[i] ), ∀Bl ∈ Nk (i) \ {Bi };
k
k
j
// Lower Qk to satisfy (25)
j
ξ̄k
[i]
, 1);
j
∀l6=i Θ[l]Qk [i,l]
0
0
0
4: Q [i, l] = max(¯
Q [i], ¯Q [l]), ∀Bl ∈ Nk (i) \ {Bi };
j
j
0
5: Qk [i, l] := Q [i, l]Qk [i, l], ∀Bl ∈ Nk (i) \ {Bi };
j
ξ̄k
[i]
, 1);
6: ¯Q [i] = min( P
j
∀l6=i Θ[l]Qk [i,l]

3:

¯0Q [i] = min( P

Q [i, l] = min(¯
Q [i], ¯Q [l]), ∀Bl ∈ Nk (i) \ {Bi };
Qjk [i, l] := Q [i, l]Qjk [i, l], ∀Bl ∈ Nk (i) \ {Bi };
// Compute Pkj
j
j
9: Pk [i, l] = Θ[l]Qk [i, l], ∀Bl ∈ Nk (i) \ {Bi };
P
j
j
10: Pk [i, i] = 1 −
∀l6=i Pk [i, l];
7:
8:

C. Example III: Local-information-based Quorum Model
This subsection shows that the proposed framework is
able to incorporate a quorum model, which is introduced in
[14], [15]. In this model, if a bin is overpopulated above
a certain level of predefined threshold called quorum, the
probabilities that agents in the bin move to neighbour bins
are temporarily increased, rather than following given Pkj .
This feature eventually brings an advantage to the convergence
performance of the swarm.
To this end, we set the secondary guidance matrix Skj as
follows: ∀i, l ∈ {1, ..., nbin } and ∀j ∈ A,
(
1/|Nk (i)|
if Ck [i, l] = 1
j
Sk [i, l] :=
(26)
0
otherwise.
This matrix makes agents in a bin equally disseminated over
its neighbour bins. In addition, the secondary feedback gain
Gjk [i] is defined as
!−1

µ̄jk [i] 
j
Gk [i] := 1 + exp γ(qi −
)
,
(27)
Θ̄[i]
where γ > 0 is a design parameter, and qi > 1 is the quorum
for bin Bi . The value of the gain is shown in Figure 3(b),
varying depending on γ and qi . As µ̄jk [i]/Θ̄[i] becomes higher
than the quorum, Gjk [i] gets close to 1 (i.e., Skj [i, l] becomes
more dominant than Pkj [i, l]). The steepness of the function at
the quorum value is regulated by γ.
The existing quorum models in [14], [15] require each agent
to know µ?k [i], which implies that the total number of agents

10

1

1
=1
=5
= 10
= 100

0.9
0.8

0.9

q i = 2,

=1

q i = 2,

= 10

0.8

q i = 2,

= 100

0.7

0.7

0.6

0.6

0.5

0.5

0.4

0.4

0.3

0.3

0.2

0.2

0.1

q i = 1.5,

= 10

0.1

0
-2

0
-1.5

-1

(underpopulated)

-0.5

0

0.5

1

1.5

2

(overpopulated)

0

0.5
1
(underpopulated)

(a) P1

1.5

2

2.5

3

3.5
4
(overpopulated)

(b) the Quorum Model
Gjk [i]

Fig. 3. The secondary feedback gains
depending on the associated
design parameters: (1) for P1 (Eqn. (18)); (2) for the quorum model (Eqn.
(27))

nA
k should be tracked in real time. It could be possible that
some agents in a swarm unexpectedly become faulted by
internal or external effects during a mission, which hinders
for other alive agents from keeping track of nA
k in a timely
manner. On the contrary, this requirement is not the case for
the quorum model in this subsection, and it works by using
the local information available from ANk (i) .
Algorithm 4 Subroutine of Algorithm 1 (Line 7-8) for the
quorum-based method
Compute Skj [i, l] ∀l using (26);
j
2: Compute Gk [i] using (27);
1:

V. A SYNCHRONOUS I MPLEMENTATION
A synchronous process induces extra time delays and interagent communications to make entire agents, who may have
different timescales for obtaining new information and make
decisions, remain in sync. Such unnecessary waiting time and
communications may cause unfavourable effects on mission
performance or even may not be realisable in practice [23].
In the previous sections, it was assumed that a swarm of
agents act synchronously at every time instant. Here we show
that the proposed framework allows agents to operate in an
asynchronous manner, assuming that the union of underlying
graphs of the corresponding Markov matrices across some time
intervals is frequently and infinitely strongly-connected.
Suppose that an algorithm to compute Pkj that satisfies (R1)(R5) in a synchronous environment is given (e.g., Algorithm 2
or 3). We propose an asynchronous implementation, as shown
in Algorithm 5, which substitutes for Line 6 in Algorithm
Algorithm 5 Asynchronous Construction of Pkj [i, l] (Substitute for Line 6 of Algorithm 1)
1:
2:
3:
4:
5:
6:
7:

+
if Bi ∈ R+
k & isnonempty(Rk \ {Bi }) then
j
Compute Pk [i, l] as usual, ∀Bl ∈ R+
k \ {Bi };
−
Pkj [i, l] = 0, ∀B
∈
R
;
l
k
P
Pkj [i, i] = 1 − ∀l6=i Pkj [i, l];
else
Pkj [i, i] = 1; Pkj [i, l] = 0, ∀l 6= i;
end if

1. We refer to a set of bins where agents are ready to use
their respective local information (e.g., µ̄jk [i]) as R+
k , and a
set of the other bins as R−
k . It is assumed that each agent
j in bin Bi ∈ R+
k also knows the local information of its
neighbour bin Bl ∈ Nk (i) if Bl ∈ R+
k . As shown in Line 2, the
agent follows an existing procedure as long as all information
required to generate Pkj [i, l] is available (e.g., µ̄jk [i] and µ̄jk [l]
for Algorithm 2, and ¯0Q [i] and ¯0Q [l] for Algorithm 3). On the
contrary, if any local information of its neighbour bin Bl ∈
Nk (i) is not available, the probability to transition to the bin
is set as zero (Line 3). In the meantime, each agent for whom
necessary local information is not ready does not deviate but
remains at the bin it belonged to. Equivalently, it can be said
that Pkj [i, i] = 1 and Pkj [i, l] = 0, ∀l 6= i (Line 6).
Hereafter, for the sake of differentiation from the original
Pkj generated in a synchronous environment, let us refer to the
matrix resulted by Algorithm 5 as asynchronous primary guidance matrix, denoted by P̄kj . Accordingly, the asynchronous
Markov matrix can be defined as:
M̄kj := (I − Wkj )P̄kj + Wkj Skj .
Here, we show that this asynchronous Markov process also
converges to the desired swarm distribution.
Lemma 3. The matrix P̄k , for every time instant k, satisfies
the following properties: (1) row-stochastic; (2) all diagonal
elementsP
are positive, and all other elements are non-negative;
nbin
and (3) i=1
Θ[i]P̄kj [i, l] = Θ[l], ∀l.
Proof. The matrix P̄kj is row-stochastic because of Line 4 and
6 in Algorithm 5. Furthermore, given that Pkj satisfies (R2),
the property (2) is valid for P̄kj because P̄kj [i, i] ≥ Pkj [i, i] for
∀i.
Let us now
to the property (3). For ∀Bi ∈ R−
k , it is
Pnturn
bin
trivial that l=1 Θ[l]P̄kj [l, i] = Θ[i] because of Line 6. For
j
∀Bi ∈ R+
k , it turns out from Algorithm 5 that (i) P̄k [i, l] =
j
j
j
−
P̄k [l, i] = 0 for ∀Bl ∈ Rk ; (ii) P̄k [i, l] =PPk [i, l] for ∀Bl ∈
j
j
R+
Pkj [i, l].
k \ {Bi } and (iii) P̄k [i, i] = Pk [i, i] +
∀l:Bl ∈R−
k
We apply the findings into the following equation:
n
bin
X

Θ[l]P̄kj [l, i] =

X

Θ[l]P̄kj [l, i]

∀l:Bl ∈R−
k

l=1

+

X

Θ[l]P̄kj [l, i] + Θ[i]P̄kj [i, i].

∀l:Bl ∈R+
k \{Bi }

(28)
The first term of the right hand side becomes zero because of
j
j
(i). Due to (ii) and the fact that Θ[i]P
P k [i, l] = Θ[l]Pjk [l, i]
∀l, the second term becomes Θ[i] ∀l:Bl ∈R+ \{Bi } Pk [i, l].
P k
The last term becomes Θ[i]Pkj [i, i] + Θ[i] ∀l:Bl ∈R− Pkj [i, l]
k
because of (iii). Putting
Equation (28)
P all of them together,
is equivalent to Θ[i] ∀l:Bl ∈R+ \{Bi } Pkj [i, l] + Θ[i]Pkj [i, i] +
k
P
Pnbin j
Θ[i] ∀l:Bl ∈R− Pkj [i, l] = Θ[i] l=1
Pk [i, l] = Θ[i].
k

Lemma 4. If the union of a set of underlying graphs of
{P̄k1 , P̄k1 +1 , ..., P̄k2 −1 } is strongly-connected, then the matrix
product P̄k1 ,k2 := P̄k1 P̄k1 +1 · · · P̄k2 −1 is irreducible.

11

Proof. Thanks to Lemma 3 and 4, the matrix product P̄ki ,ki+1
for each time interval [ki , ki+1 ) satisfies (R1)-(R4). Therefore,
one can prove this theorem by similarly following the proof
of Theorem 1.
VI. N UMERICAL E XPERIMENTS
A. Effects of Primary Local-feedback Gain ξ¯kj [i]
Depending on the shape of primary feedback gain ξ¯kj [i], the
performance of the proposed framework changes, especially
with respect to convergence rate, fraction of transitioning
agents, and residual convergence error. Let us first investigate
the effect of changes in the feedback gain using Algorithm 2
with Equation (20).
We consider a scenario where a set of 2, 000 agents are
supposed to be distributed over an arena consisting of 10 × 10
bins, as depicted in Figure 1. There are vertical and horizontal
paths between adjacent bins. Note that the agents are allowed
to move at most 3 paths away within a unit time instant. All
the agents start from a bin, which reflects the fact that they
are generally deployed from a base station at the beginning
of a mission. The desired swarm distribution Θ is uniformrandomly generated at each scenario. The agents are assumed
to estimate necessary information correctly, e.g. µ̄jk [i] = µ̄?k [i].
For the rigorous validation, the performance of the proposed
algorithm will be compared with that of the GICA-based
algorithm [13]. To this end, f (Ek [i, l]) is set to be the same
as the corresponding coefficient in [13, Corollary 1]:
f (Ek [i, l]) := 1 −

Ek [i, l]
Ek,max + E

(29)

where Ek,max is the maximum element of the travelling
expense matrix Ek , and E is a user-design parameter. Ek [i, l]
is defined as a linear function based on the distance between
bin Bi and Bl :
Ek [i, l] := E1 · ∆s(i,l) + E0

25

Fraction of Transitioning Agents (%)

Primary feedback gain with varying

0.8
0.7
0.6
0.5
0.4
= 0.2
= 0.4
= 0.6
= 0.8
= 1.0
= 1.2

0.3
0.2
0.1
0

= 0.2
= 0.4
= 0.6
= 0.8
= 1.0
= 1.2

0.2

0.4

0.6

Primary feedback gain when

15

10

5

0.8

1

0

500

1000

1500

=1

2000

2500

3000

3500

4000

Time Instant

(a) Primary Local-feedback Gain

(b) Fraction
Agents

1

of

Transitioning

0.08
α
α
α
α
α
α

0.9
0.8
0.7

= 0.2
= 0.4
= 0.6
= 0.8
= 1.0
= 1.2

α
α
α
α
α
α

0.07

0.6
0.5
0.4
0.3
0.2

0.06

= 0.2
= 0.4
= 0.6
= 0.8
= 1.0
= 1.2

0.05
0.04
0.03
0.02
0.01

0.1
0

20

0
0

Convergence Error ( D H(Θ, µ *k) )

Theorem 3. Suppose that there exists an infinite sequence of non-overlapping time intervals [ki , ki+1 ), i =
0, 1, 2, ..., such that the union of underlying graphs
of {P̄ki , P̄ki +1 , ..., P̄ki+1 −1 } in each interval is stronglyconnected. Let the stochastic state of agent j at time instant k ≥ k0 , governed by the corresponding Markov process from an arbitrary state xjk0 , be xjk = xjk0 Ūkj0 ,k :=
j
xjk0 M̄kj0 M̄kj0 +1 · · · M̄k−1
. Then, it holds that limk→∞ xjk = Θ
pointwise for all agents, irrespective of the initial condition.

1
0.9

Convergence Error ( D H(Θ, µ *k) )

Proof. Since the union of a set of underlying graphs of
{P̄k1 , P̄k1 +1
P̄k2 −1 } is strongly-connected, the underlying
P,k...,
2 −1
graph of
k=k1 P̄k is also strongly-connected. Noting that
every P̄k , ∀k ∈ {k1 , k1 + 1, ..., k2 − 1} is a nonnegative
nbin × nbin matrix and its diagonal elements are positive (by
Lemma 3), it follows from [36, Lemma 2] that P̄k1 ,k2 ≥
Pk2 −1
γ k=k
P̄k , where γ > 0. This implies that the underlying
1
graph of P̄k1 ,k2 is strongly-connected, and thus the matrix
P̄k1 ,k2 is irreducible.

0

500

1000

1500

2000

2500

3000

3500

4000

0
3500

3600

3700

Time Instant

(c) Convergence Error

3800

3900

4000

Time Instant

(d) Convergence Error (zoomed-in)

Fig. 4. Performance comparison depending on the primary local-feedback
gain ξ̄kj [i] with different setting of α: (a) the value of ξ̄kj [i]; (b) the fraction
of transitioning agents; (b) the convergence performance; (d) the convergence
performance (zoomed-in for time instant between 2000 and 4000)

are assumed to follow any shortest route when they transition
between two bins. The design parameters are set as follows:
E1 = 1 and E0 = 0.5 in (30); E = 0.1 in (29); ξ = 10−9
in (6); M = 1 in (20); β = 1.8 × 105 in (18); and τ j = 10−6
in (9).
As a performance index for the closeness between the current swarm distribution µ?k and Θ, we use Hellinger Distance,
i.e.,
v
unbin 
2
q
uX p
1
?
Θ[i] − µ?k [i] ,
DH (Θ, µk ) := √ t
2 i=1
Hellinger Distance is known as a “concept of measuring
similarity between two distributions” [37] and is utilised as
a feedback gain in the existing work [13].
More importantly, to examine the effects of the shape of
ξ¯kj [i], we set α in (6) as 0.2, 0.4, 0.6, 0.8, 1 and 1.2.
Figure 4 reveals that the convergence rate can be traded off
against the fraction of transitioning agents and the residual
convergence error. As ξ¯kj [i] becomes more concave, i.e. the
value of α decreases, the summation of off-diagonal entries
of Pkj becomes higher, leading to more transitioning agents,
but a faster convergence rate. At the same time, a higher value
of ξ¯kj [i] even at a low value of |Θ̄[i]−µjk [i]| gives rise to unnecessarily higher off-diagonal entries of Pkj . Hence, the swarm
tends to be prevented from converging to the desired swarm
distribution properly, resulting in higher residual convergence
error.
B. Comparison with a GICA-based Method

(30)

where ∆s(i,l) is the minimum required number of paths from
Bi to Bl ; E1 and E0 are user-design parameters. The agents

Let us now compare the LICA-based method for (P1) with
the GICA-based method in [13]. The scenario considered is the
same as the one in the previous subsection except for α = 0.6.

12

1

25
LICA
GICA

4

LICA
GICA

0.8
0.7
0.6
0.5
0.4
0.3
0.2

LICA
GICA

20

Cumulative Travel Expenses

Fraction of Transitioning Agents (%)

Convergence Error (DH ( , *k ))

0.9

10

2.5

15

10

5

2

1.5

1

0.5

0.1
0

0
0

50

100

150

200

250

300

350

400

0
0

50

100

Time Instant

150

200

250

300

350

Time Instant

(a) Convergence Error

(b) Fraction of Transitioning Agents

400

0

50

100

150

200

250

300

350

400

Time Instant

(c) Cumulative Travel Expenses

Ratio of Cumulative Travel Expenses until DH = 0.1

Ratio of Required Time Instants for Convergence (LICA/GICA)

Fig. 5. Performance comparison between the proposed method (LICA) with the existing method (GICA): (a) the error between the current swarm status and
the desired status, shown as Hellinger Distance; (b) the fraction of agents who transitions between any two bins; (c) the cumulative travel expenses of all the
agents from the beginning.

1.4

1.3

1.2

1.1

1

0.9

0.8

method and those by the compared method are 1.72 × 104 and
1.96 × 104 , respectively, and their ratio is 0.878. This is also
confirmed by the statistical result in Figure 6(b). A possible
explanation is that when some of the bins do not meet their
desired swarm densities, the entire agents in the GICA-based
method would obtain higher feedback gains, which might lead
to unnecessary transitions. On the contrary, this is not the case
in the LICA-based method since agents are only affected by
their neighbour bins.

0.92
0.9
0.88
0.86
0.84
0.82
0.8
0.78
0.76
0.74

0.7
0.3

0.28 0.26 0.24 0.22

0.2

0.18 0.16 0.14 0.12

Convergence Error (DH ( , *k ))

(a) Converged Time Instants

0.1

1

LICA/GICA

(b) Cumulative
Expenses

Travel

Fig. 6. Performance comparison (Monte-Carlo experiments) between the
proposed method (LICA) and the existing method (GICA): (a) the required
time instants to converge to DH (Θ, µ?k ) ∈ {0.30, 0.28, ..., 0.12, 0.10} (i.e.,
convergence rate); (2) the ratio of the cumulative travel expenses by LICA to
those by GICA until converging to DH (Θ, µ?k ) = 1.

Note that Θ in Remark 4 can control convergence rate, but is
not discussed in [13]. For the fair comparison, Θ is applied
to both the methods.
We conduct 100 runs of Monte Carlo experiments. Figure
5 presents the results of one representative scenario and the
statistical results of the Monte Carlo experiments are shown
in Figure 6. According to Figure 5(a), the convergence rate
of the proposed method is slower at the initial phase, but
similar to that of the the GICA-based method as reaching
DH (Θ, µ?k ) = 0.10. This is confirmed by the statistical results
in Figure 6(a), where the ratio of the required time instants
for converging to DH (Θ, µ?k ) ∈ {0.30, 0.28, ..., 0.12, 0.10} in
the LICA-based method to those of the GICA-based method
is presented. At this point, it is worth noting that these convergence rate results are presented in respect to time instants
of each Markov process. As the LICA-based framework may
have a much shorter timescale, its convergence performance in
practice could be better than that of the GICA-based method.
Figure 5(c) shows that the cumulative travel expenses are
smaller in the proposed method. The expenses by the proposed

C. Robustness in Asynchronous Environments
This subsection investigates the effects of asynchronous
environments in the proposed LICA-based method for (P1)
and compares them with those in the GICA-based method
in [13]. Hence, a realistic scenario where an asynchronous
process is required is considered: in the scenario, it is assumed
that agents in some bins cannot communicate for some reason
(such bins are called blocked) and thus other agents in normal
bins have to perform their own process without waiting them.
The proportion of blocked bins to the entire bins is set to
be different values, i.e. 0%, 10%, 20% and 30%. At each time
instant, the corresponding proportion of bins are randomly
selected as blocked bins. For the proposed framework, the
asynchronous implementation in Section V is built upon
Algorithm 2. In the GICA-based method, for the comparison
purpose only, it is assumed that agents in normal bins obtain
µjk = µ?k without interacting with agents in the blocked bins.
The rest of scenario setting are the same as those in Section
VI-B.
Figure 7 illustrates the performance of each method: convergence rate, fraction of transitioning agents, and cumulative travel expenses. As the proportion of the blocked bins
increases, the GICA-based method tends to have faster convergence speed, whereas it loses Desired Feature 2 and thus
increases cumulative travelling expenses (as shown in Figure
7(a), 7(b), and 7(c), respectively). On the contrary, the LICAbased method shows a graceful degradation in terms of Desired
Feature 2 (as shown in Figure 7(b)). A possible explanation
for these results could be that higher feedback gains due to the
communication disconnection induce faster convergence performance in each method than the normal situation. This effect

13

30

Convergence Error (DH ( , *k ))

0.8
0.7
0.6

Fraction of Transitioning Agents (%)

LICA (Blocked 0%)
LICA (Blocked 10%)
LICA (Blocked 20%)
LICA (Blocked 30%)
GICA (Blocked 0%)
GICA (Blocked 10%)
GICA (Blocked 20%)
GICA (Blocked 30%)

0.9

0.5
0.4
0.3
0.2

10

14
LICA (Blocked 0%)
LICA (Blocked 10%)
LICA (Blocked 20%)
LICA (Blocked 30%)
GICA (Blocked 0%)
GICA (Blocked 10%)
GICA (Blocked 20%)
GICA (Blocked 30%)

25

20

4

LICA (Blocked 0%)
LICA (Blocked 10%)
LICA (Blocked 20%)
LICA (Blocked 30%)
GICA (Blocked 0%)
GICA (Blocked 10%)
GICA (Blocked 20%)
GICA (Blocked 30%)

12

Cumulative Travel Expenses

1

15

10

5

10

8

6

4

2

0.1
0

0
0

50

100

150

200

250

300

350

400

450

0
0

50

100

Time Instant

150

200

250

300

350

400

450

0

50

100

Time Instant

(a) Convergence Error

150

200

250

300

350

400

450

Time Instant

(b) Fraction of Transitioning Agents

(c) Cumulative Travel Expenses

80

Swarm Distribution (D H ( , *k ))

Fractions of Transitioning Agents (%)

P2
P2 + Quorum

0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0

P2
P2 + Quorum

70
60
50
40
30
20
10
0

0

500

1000

1500

2000

2500

3000

3500

4000

Time Instant

(a) Convergence Error

0

500

1000

1500

2000

2500

3000

3500

Time Instant

(b) Fraction of Transitioning Agents

4000

40

400
Actual (Case 1)
Probabilistic (Case 1)
Actual (Case 2)
Probabilistic (Case 2)

35

350

30

300

25

250

20

200

15

150

10

100

5

50

0
0

500

1000

1500

2000

2500

3000

3500

0
4000

Max (# Transitions of Agents for Each Path) - Case 2

1

Max (# Transitions of Agents for Each Path) - Case 1

Fig. 7. Performance comparison in communication-disconnected situations: (a) the error between the current swarm status and the desired status, shown as
Hellinger Distance; (b) the fraction of agents who transitions between any two bins; (c) the cumulative travel expenses of all the agents from the beginning.

Time Instant

(c) Maximum Number of Transitioning Agents via
Each Path (for P2)

Fig. 8. Comparison results between the method for P2 and the quorum-based model: (a) the error between the current swarm status and the desired status,
shown as Hellinger Distance, at each time instant; (b) the fraction of agents who transition between any two bins at each time instant. (c) The maximum
number of transitioning agents via each path in the method for P2 (Case 1: |A| = 10, 000 and c(i,l) = 20, ∀i 6= l; Case 2: |A| = 100, 000 and c(i,l) = 200,
∀i 6= l)

is dominant for the GICA-based method because it affects
the entire agents, who use global information. However, in
the LICA-based framework, the communication disconnection
only locally influences so that its effectiveness is relatively
modest.
D. Demonstration of Example II and III
This subsection demonstrates the LICA-based method for
(P2) (i.e., Algorithm 3) and the quorum model (i.e., Algorithm
4). For the former, we consider a scenario where 10, 000 agents
and an arena consisting of 10 × 10 bins are given. The arena
is as depicted in Figure 1, where the agents are allowed to
move only one path away within a unit time instant. For each
one-way path, the upper flux bound per time instant is set as
20 agents (i.e., c(i,l) = 20, ∀i 6= l). All the agents start from a
bin, and the desired swarm distribution is uniform-randomly
generated.
For the latter, we build the quorum model upon the LICAbased method for (P2). This can be a good strategy for a user
who wants to achieve not only faster convergence rate but
also lower unnecessary transitions after equilibrium, which are
regulated by the upper flux bounds. Thus, in the same scenario
described above, we will demonstrate the combined algorithm

that computes Skj and Gjk by Algorithm 4 and Pkj by Algorithm
3. We set qi = 1.3 and γ = 30 for (27); α = 1 and ξ = 10−9
for (6); and τ j = 10−6 for (9).
Figure 8(a) and 8(b) presents that the both approaches make
the swarm converge to the desired swarm distribution. It is
observed that the number of transitioning agents in the method
for (P2) are restricted because of the upper flux bound during
the entire process. Meanwhile, the quorum-based method very
quickly disseminates the agents, who are initially at one bin,
over other bins, and thus the fraction of transitioning agents
is very high at the initial phase. After that, the population
fraction by the quorum-based method drops and remains as
low as that by the method for (P2).
Figure 8(c) presents the maximum value amongst the number of transitioning agents via each (one-way) path. The
red line indicates the actual result by the method for (P2),
while the green line indicates the corresponding probabilistic
value (i.e., max∀i max∀l6=i nk [i]Pkj [i, l]). It is shown that the
stochastic decision policies reflect the given upper bound,
meanwhile this bound is often violated in practice due to the
finite cardinality of the agents. However, the result in the same
scenario with setting |A| = 100, 000 and c(i,l) = 200, ∀i 6= l
(denoted by Case 2), depicted by the blue and magenta lines
in Figure 8(c), suggests that such violation can be mitigated

14

as the number of given agents increases.
VII. C ONCLUSION
This paper poposed a LICA-based closed-loop-type framework for probabilistic swarm distribution guidance. Since the
feedback gains can be generated based on local information,
agents have shorter and different timescales for using new
information, and can incorporate an asynchronous decisionmaking process. Even using local information, the proposed
framework converges to a desired density distribution, while
maintaining scalability, robustness, and long-term system efficiency. Furthermore, the numerical experiments have showed
that the proposed framework is suitable for a realistic environment where communication between agents is partially and
temporarily disconnected. This paper has explicitly presented
the design requirements for the Markov matrix to hold all these
advantages, and has provided specific problem examples of
how to implement this framework.
Future works include optimisation of ξ¯kj [i], which can mitigate the trade-off between convergence rate and residual error.
In addition, it is expected that the communication cost required
for the proposed framework can be reduced by incorporating
a vision-based local density estimation [38].
ACKNOWLEDGEMENT
The authors gratefully acknowledge that this research was
supported by International Joint Research Programme with
Chungnam National University (No. EFA3004Z). Thanks to
Sangjun Bae for supportive discussions.
A PPENDIX
A. Regarding the Convergence Analysis in Theorem 1
Definition 9 (Irreducible). A matrix is reducible if and only
if its associated digraph is not strongly connected. A matrix
that is not reducible is irreducible.
Definition 10 (Primitive). A primitive matrix is a square
nonnegative matrix A such that for every i, j there exists k > 0
such that Ak [i, j] > 0.
Definition 11 (Regular). A regular matrix is a stochastic
matrix such that all the entries of some power of the matrix
are positive.
Definition 12. [33, pp.92, 149] [13] (Asymptotic Homogeneity) “A sequence of stochastic matrices Mk ∈ Mn×n , k ≥ k0 ,
is said to be asymptotically homogeneous (with respect to
d) if there exists a row-stochastic vector d ∈ Pn such that
limk→∞ dMk = d.”
Definition 13. [33, pp.92, 149] [13] (Strong Ergodicity)
The matrix product Uk0 ,k := Mk0 Mk0 +1 · · · Mk−1 , formed
from a sequence of stochastic matrices Mr ∈ Pn×n , r ≥
k0 , is said to be strongly ergodic if for each i, l, r we
get limr→∞ Uk0 ,r [i, l] = v[l], where v ∈ Pn is a rowstochastic vector. Here, v is called its unique limit vector (i.e.,
limr→∞ Uk0 ,r = 1> v).

Lemma 5. Given the requirements (R1)-(R4) are satisfied, Mkj
in Equation (11) has the following properties:
1) row-stochastic;
2) irreducible;
3) all diagonal elements are positive, and all other elements
are non-negative;
4) there is a positive lower bound γ such that 0 < γ ≤
j
+
min+
refers to the minimum
i,l Mk [i, l] (Note that min
of the positive elements);
5) asymptotically homogeneous with respect to Θ.
In addition, Ukj0 ,k in Equation (12) has the following
properties:
6) irreducible;
7) all diagonal elements are positive, and all other elements
are non-negative;
8) primitive [39, Lemma 8.5.4, p.541]
Proof. This lemma can be proved by similarly following the
mathematical development for [13, Theorem 4]. Mkj is rowstochastic because Mkj · 1> = (I − Wkj )Pkj · 1> + Wkj Skj · 1> =
(I − Wkj ) · 1> + Wkj · 1> = 1> . Pkj is irreducible and
ωkj [i] is always less than 1, thus Mkj is also irreducible
(i.e., Mkj [i, l] > 0 if Pkj [i, l] > 0). The property 3) is true
because diag(I − Wkj ) > 0, Wkj ≥ 0, and Pkj is also a nonnegative matrix such that its diagonal elements are positive.
The property 4) is implied by either the property 2) or 3). From
the definition of Wkj , it follows that limk→∞ Wkj = 0 (because
of exp(−τ j k)), and thereby limk→∞ Mkj = limk→∞ Pkj .
Hence, limk→∞ ΘMkj = limk→∞ ΘPkj = Θ, and the property
5) is valid.
Let us now turn to Ukj0 ,k . It is irreducible due to the
fact that if Mrj [i, l] > 0 for some r ∈ {k0 , ..., k − 1} and
∀i, l ∈ {1, ..., nbin }, then the corresponding element Ukj0 ,k [i, l]
is greater or equal to the product of positive diagonal elements
and Mrj [i, l] and the property 2). The property 7) is true
because Ukj0 ,k is a product of nonnegative matrices where all
diagonal entries are positive. It follows from [39, Lemma 8.5.4,
p.541] that Ukj0 ,k is primitive: “if a square matrix is irreducible,
nonnegative and all its main diagonal entries are positive, then
the matrix is primitive”.
R EFERENCES
[1] M. Brambilla, E. Ferrante, M. Birattari, and M. Dorigo, “Swarm
Robotics: a Review from the Swarm Engineering Perspective,” Swarm
Intelligence, vol. 7, no. 1, pp. 1–41, 2013.
[2] A. Kolling, P. Walker, N. Chakraborty, K. Sycara, and M. Lewis,
“Human Interaction With Robot Swarms: A Survey,” IEEE Transactions
on Human-Machine Systems, vol. 46, no. 1, pp. 9–26, 2016.
[3] E. Sahin, “Swarm Robotics: From Sources of Inspiration to Domains of
Application,” in Swarm Robotics. Berlin: Springer, 2005, pp. 10–20.
[4] B. Acikmese and D. S. Bayard, “A Markov Chain Approach to Probabilistic Swarm Guidance,” in American Control Conference, Montreal,
Canada, 2012, pp. 6300–6307.
[5] ——, “Probabilistic Swarm Guidance for Collaborative Autonomous
Agents,” in American Control Conference, Portland, OR, USA, 2014,
pp. 477–482.
[6] ——, “Markov Chain Approach to Probabilistic Guidance for Swarms
of Autonomous Agents,” Asian Journal of Control, vol. 17, no. 4, pp.
1105–1124, 2015.
[7] I. Chattopadhyay and A. Ray, “Supervised Self-Organization of Homogeneous Swarms Using Ergodic Projections of Markov Chains,” IEEE
Transactions on Systems, Man, and Cybernetics, Part B (Cybernetics),
vol. 39, no. 6, pp. 1505–1515, 2009.

15

[8] N. Demir and B. Acikmese, “Probabilistic Density Control for Swarm
of Decentralized ON-OFF Agents with Safety Constraints,” in American
Control Conference, Chicago, IL, USA, 2015, pp. 5238–5244.
[9] R. Luo, N. Chakraborty, and K. Sycara, “Supervisory Control for CostEffective Redistribution of Robotic Swarms,” in IEEE International
Conference on Systems, Man, and Cybernetics (SMC), San Diego, CA,
USA, 2014, pp. 596–601.
[10] N. Demir, U. Eren, and B. Açkmee, “Decentralized Probabilistic Density
Control of Autonomous Swarms with Safety Constraints,” Autonomous
Robots, vol. 39, no. 4, pp. 537–554, 2015.
[11] S. Bandyopadhyay, S.-J. Chung, and F. Y. Hadaegh, “Inhomogeneous
Markov Chain Approach To Probablistic Swarm Guidance Algorithms,”
in 5th Int. Conf. on Spacecraft Formation Flying Missions and Technologies (SFFMT), Munich, Germany, 2013.
[12] D. Morgan, G. P. Subramanian, S. Bandyopadhyay, S.-J. Chung, and
F. Y. Hadaegh, “Probabilistic Guidance of Distributed Systems using
Sequential Convex Programming,” in IEEE/RSJ International Conference on Intelligent Robots and Systems, Chicago, IL, USA, 2014, pp.
3850–3857.
[13] S. Bandyopadhyay, S.-J. Chung, and F. Y. Hadaegh, “Probabilistic and
Distributed Control of a Large-Scale Swarm of Autonomous Agents,”
IEEE Transactions on Robotics, vol. 33, no. 5, pp. 1103–1123, 2017.
[14] A. Halasz, M. A. Hsieh, S. Berman, and V. Kumar, “Dynamic Redistribution of a Swarm of Robots among Multiple Sites,” in IEEE/RSJ
International Conference on Intelligent Robots and Systems, San Diego,
CA, USA, 2007, pp. 2320–2325.
[15] M. A. Hsieh, A. Halasz, S. Berman, and V. Kumar, “Biologically
Inspired Redistribution of a Swarm of Robots among Multiple Sites,”
Swarm Intelligence, vol. 2, no. 2-4, pp. 121–141, 2008.
[16] S. Berman, A. Halasz, M. Ani Hsieh, and V. Kumar, “Navigation-based
Optimization of Stochastic Strategies for Allocating a Robot Swarm
among Multiple Sites,” in 47th IEEE Conference on Decision and
Control, Cancun, Mexico, 2008, pp. 4376–4381.
[17] S. Berman, A. Halasz, M. A. Hsieh, and V. Kumar, “Optimized
Stochastic Policies for Task Allocation in Swarms of Robots,” IEEE
Transactions on Robotics, vol. 25, no. 4, pp. 927–937, 2009.
[18] T. W. Mather and M. A. Hsieh, “Macroscopic Modeling of Stochastic
Deployment Policies with Time Delays for Robot Ensembles,” The
International Journal of Robotics Research, vol. 30, no. 5, pp. 590–
600, 2011.
[19] K. Lerman, A. Martinoli, and A. Galstyan, “A Review of Probabilistic
Macroscopic Models for Swarm Robotic Systems,” in Swarm Robotics.
Berlin: Springer, 2005, pp. 143–152.
[20] L. B. Johnson, H.-L. Choi, and J. P. How, “The Role of Information Assumptions in Decentralized Task Allocation: A Tutorial,” IEEE Control
Systems, vol. 36, no. 4, pp. 45–58, 2016.
[21] S. Xue, C. Sun, J. Zeng, Y. Jin, and R. Cheng, “Effect of Communication
Modes to Swarm Robotic Search,” The Open Electrical & Electronic
Engineering Journal, vol. 8, no. 1, pp. 240–244, 2014.
[22] B. I. Koh, A. D. George, R. T. Haftka, and B. J. Fregly, “Parallel
Asynchronous Particle Swarm Optimization,” International Journal for
Numerical Methods in Engineering, vol. 67, no. 4, pp. 578–595, 2006.
[23] L. Johnson, S. Ponda, H.-L. Choi, and J. How, “Asynchronous
Decentralized Task Allocation for Dynamic Environments,” in Infotech@Aerospace 2011, St.Louis, MO, USA, 2011.
[24] I. D. Couzin, J. Krause, R. James, G. D. Ruxton, and N. R. Franks,
“Collective Memory and Spatial Sorting in Animal Groups.” Journal of
theoretical biology, vol. 218, no. 1, pp. 1–11, 2002.
[25] I. D. Couzin, J. Krause, N. R. Franks, and S. A. Levin, “Effective
Leadership and Decision-making in Animal Groups on the Move,”
Nature, vol. 433, no. 7025, pp. 513–516, 2005.
[26] J. Gautrais, C. Jost, and G. Theraulaz, “Key Behavioural Factors in a
Self-Organised Fish School Model,” Annales Zoologici Fennici, vol. 45,
no. 5, pp. 415–428, 2008.
[27] D. J. Hoare, I. D. Couzin, J. G. J. Godin, and J. Krause, “Contextdependent Group Size Choice in Fish,” Animal Behaviour, vol. 67, no. 1,
pp. 155–164, 2004.
[28] S. Bandyopadhyay and S.-J. Chung, “Distributed Estimation using
Bayesian Consensus Filtering,” American Control Conference, pp. 634–
641, 2014.
[29] T. D. Seeley, The Wisdom of the Hive. Cambridge, Massachusetts:
Havard University Press, 1995.
[30] L. Keller, M. J. B. Krieger, and J.-B. Billeter, “Ant-like Task Allocation
and Recruitment in Cooperative Robots,” Nature, vol. 406, no. 6799,
pp. 992–995, 2000.
[31] B. L. Partridge, “The Structure and Function of Fish Schools,” Scientific
American, vol. 246, no. 6, pp. 114–123, 1982.

[32] C. Becco, N. Vandewalle, J. Delcourt, and P. Poncin, “Experimental
Evidences of a Structural and Dynamical Transition in Fish School,”
Physica A: Statistical Mechanics and its Applications, vol. 367, pp. 487–
493, 2006.
[33] E. Seneta, Non-negative Matrices and Markov Chains, ser. Springer
Series in Statistics. New York, NY: Springer New York, 1981.
[34] Y. Bestaoui Sebbane, Planning and Decision Making for Aerial Robots,
ser. Intelligent Systems, Control and Automation: Science and Engineering. Cham: Springer International Publishing, 2014, vol. 71.
[35] I. C. F. Ipsen and T. M. Selee, “Ergodicity Coefficients Defined by Vector
Norms,” SIAM Journal on Matrix Analysis and Applications, vol. 32,
no. 1, pp. 153–200, jan 2011.
[36] A. Jadbabaie, Jie 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, 2003.
[37] J. Chung, P. Kannappan, C. Ng, and P. Sahoo, “Measures of Distance
between Probability Distributions,” Journal of Mathematical Analysis
and Applications, vol. 138, no. 1, pp. 280–292, 1989.
[38] S. A. M. Saleh, S. A. Suandi, and H. Ibrahim, “Recent Survey on Crowd
Density Estimation and Counting for Visual Surveillance,” Engineering
Applications of Artificial Intelligence, vol. 41, pp. 103–114, 2015.
[39] R. A. Horn and C. R. Johnson, Matrix Analysis, 2nd ed. Cambridge:
Cambridge University Press, 2012.

