IEEE TRANSACTIONS ON AUTOMATIC CONTROL

1

Distributed area coverage control with imprecise
robot localization: Simulation and experimental
studies

arXiv:1612.04704v1 [] 14 Dec 2016

Sotiris Papatheodorou, Anthony Tzes Senior Member, IEEE, Konstantinos Giannousakis and Yiannis Stergiopoulos

Abstract—This article examines the area coverage problem for
a network of mobile robots with imprecise agents localization.
Each robot has uniform radial sensing ability, governed by first
order kinodynamics. The convex-space is partitioned based on
the Guaranteed Voronoi (GV) principle and each robot’s area
of responsibility corresponds to its GV-cell, bounded by hyperbolic arcs. The proposed control law is distributed, demanding
the positioning information about its GV-Delaunay neighbors.
Simulation and experimental studies are offered to highlight the
efficiency of the proposed control law.
Index Terms—Multi-agent systems, Multi-robot systems, Mobile robots, Cooperative systems, Autonomous agents.

I. I NTRODUCTION
REA coverage of a planar region by a swarm of robotic
agents is an active field of research with great interest and
potential. The task involves the dispersal inside a region of interest of autonomous, sensor–equipped robotic agents in order
to achieve sensor coverage of that region. Distributed control
algorithms are a popular choice for the task because of their
computational efficiency compared to centralized controllers,
their ability to adapt in real–time to changes in the robot
swarm and their reliance only on local information. Thus they
can be implemented on robots with low computational power
and low power radio transceivers, since each agent is required
to communicate only with its neighbors. The downside of
distributed control schemes is that they lead to local optima of
their objective function because of their lack of information
on the state of all agents.
Area coverage can be either static [1]–[3] in which case the
agents converge to some static final configuration or sweep
[4]–[6] in which the agents constantly move in order to satisfy
a constantly changing coverage objective.
Other factors that need to be taken into account are the
type of the region surveyed [7], [8], the sensing pattern and
performance of the sensors [9]–[12], or the particular dynamics
of the agents [13].
There have been varied approaches to area coverage such
as distributed optimization [14], [15], model predictive control
[16], [17], game theory [18] and optimal control [19].

A

All localization systems have an inherent uncertainty in
their measurements, thus there has been significant work on
taking this uncertainty into account either by data fusion [20],
safe trajectory planning [21] or spatial probability [22]. Our
approach creates a suitable space partioning scheme in lieu of
the agents’ positioning uncertainty and a distributed gradient–
based control law.
This article examines the problem of persistent area coverage of a convex planar region by a network of homogeneous
agents equipped with isotropic sensors. In order to take into
account the localization uncertainty, each agent is assumed
to lie somewhere inside a circular positioning uncertainty
region and Guaranteed Voronoi partitioning [23] of the region of interest is constructed from these uncertain regions.
Subsequently, a distributed control law is derived based on
the aforementioned partitioning so that a coverage objective
increases monotonously. Because of the complexity of that
law, a suboptimal one is proposed and both laws are compared
through simulations. Additionally, the suboptimal control law
is evaluated by two experiments using differential drive robots.
This article is an extension of [24] which did not contain
simulation results from the optimal control law and lacked
the experimental evaluation of the suboptimal control law.
The layout of the article is as follows. The required mathematical preliminaries are presented in Section II. The Voronoi
and the Guaranteed Voronoi space partitions are examined
in Section III. Section IV contains the definition of the
coverage objective and the derivation of the gradient ascent
based control law that maximizes it. Additionally it contains
the proposed suboptimal control law. Simulation studies are
presented in Section V in order to compare and evaluate both
control laws. Section VI contains experimental results from the
implementation of the suboptimal control law and is followed
by concluding remarks.
II. M ATHEMATICAL P RELIMINARIES
Let us assume a compact convex region Ω ∈ R2 under
surveillance and a network of n identical dimensionless mobile
agents (nodes). Each agent’s dynamics are governed by
q̇i = ui , ui ∈ R2 , qi ∈ Ω, i ∈ In

This work has received funding from the European Union’s Horizon 2020
Research and Innovation Programme under the Grant Agreement No.644128,
AEROWORKS.
The authors are with the Electrical & Computer Engineering Department,
University of Patras, Rio, Achaia 26500, Greece. Corresponding author’s
email: tzes@ece.upatras.gr

(1)

where ui is the control input for each agent and In =
{1, 2, . . . , n}.
The agents’ positions are not known precisely and thus in
the general case each agent has a positioning uncertainty of

IEEE TRANSACTIONS ON AUTOMATIC CONTROL

2

different measure. We assume that an upper bound riu for the
positioning uncertainty of each agent is known and thus its
center lies within a disk
Ciu (qi , riu ) = {q ∈ Ω : k q − qi k≤ riu } , i ∈ In

(2)

where k · k is the Euclidean metric and qi the agent’s position
as reported by its localization equipment.
All agents are equipped with identical omnidirectional sensors with circular sensing footprints and as such the sensed
region of each agent can be defined as

When constructing the Voronoi cell of a point qi , only its
Delaunay neighbors need to be considered, thus simplifying
the definition of its Voronoi cell into

Vi = q ∈ Ω : kq − qi k ≤ q − q j , ∀ j ∈ Ni , i ∈ In .
Using this property, a Voronoi diagram can be constructed
by first finding the Delaunay triangulation of the input points
and then for each cell Vi by finding the intersection of the
halfplanes defined by point qi and all points in Ni .

Cis (qi , rs ) = {q ∈ Ω : k q − qi k≤ rs } , i ∈ In
where rs is the common range of the sensors.
We also define for each agent a ‘Guaranteed Sensed Region’
(GSR) as the region that is sensed by the agent for all of its
possible positions inside Ciu . The GSRs are defined as
(
)
Cigs (Ciu ,Cis ) =

\

Cis (qi , rs ), ∀qi ∈ Ciu , i ∈ In

qi

Ciu

Cis

and since
and
are disks, the previous expression can be
further simplified into
Cigs (qi , rs − riu ) = {q ∈ Ω :k q − qi k≤ rs − riu } , i ∈ In ,

Fig. 1. Voronoi diagram (left) for 6 dimensionless nodes, with its equivalent
Guaranteed Voronoi diagram (right) in the case of disks (centered on those
points).

(3)

provided that rs ≥ riu .
When an agent’s position is known precisely, i.e. riu = 0,
its GSR is equivalent to its sensed region, whereas when an
agent’s positioning uncertainty is too large compared to its
/
sensing capabilities, i.e. riu > rs , then Cigs = 0.
III. S PACE PARTITIONING
Most distributed area coverage schemes assign an area
of responsibility to each agent based on information from
its neighboring agents. Then each agent is responsible for
maximizing the coverage solely on its own responsibility
region. The most common partitioning scheme is the Voronoi
partitioning, described in Section III-A. However in this article, because of the localization uncertainty of the agents, the
guaranteed Voronoi diagram described in Section III-B is used
for the assignment of areas of responsibility.
A. Voronoi Diagram
The Voronoi diagram for a set of points
Q = {q1 , q2 . . . qn }, qi ∈ Ω is defined as follows

Vi = q ∈ Ω : kq − qi k ≤ q − q j , ∀ j ∈ In , j 6= i , i ∈ In .
Each Vi ⊂ Ω is called the Voronoi cell of node i. A Voronoi
diagram of 6 points constrained in a subset Ω of R2 is shown
in Figure 1 (left), where the cell boundaries are shown in blue.
It is important to note that the Voronoi diagram is a complete
S
tessellation of Ω, that is i∈In Vi = Ω.
Another useful property of Voronoi diagrams is their duality to Delaunay triangulations. Let us define the Delaunay
neighbors of a point qi as

Ni = j ∈ In , j 6= i : Vi ∩V j 6= 0/ , i ∈ In ,
(4)

B. Guaranteed Voronoi Diagram
The Guaranteed Voronoi (GV) diagram is defined for a
set of uncertain regions D = {D1 , D2 , . . . , Dn }, Di ⊂ Ω. Each
uncertain region Di contains all the possible positions of a
point qi ∈ R2 whose localization cannot be measured precisely.
Thus the GV diagram is defined as
Vig

= {q ∈ Ω : max kq − qi k ≤ min q − q j
∀ j ∈ In , j 6= i, qi ∈ Di , q j ∈ D j , i ∈ In .

(5)

An equivalent definition of the GV diagram is as follows
Vig =

\

Hi j

i∈In

where

Hi j = q ∈ Ω : max kq − qi k ≤ min q − q j , ∀qi ∈ Di , ∀q j ∈ D j .

This definition is analogous to the construction of the classic
Voronoi cell by the halfplane intersection method described
previously, with the difference that the regions being intersected are no longer halfplanes.
As such, each cell Vig contains the points on the plane that
are closer to the corresponding point qi of region Di for all
possible configurations of the uncertain points. A Guaranteed
Voronoi diagram of 6 regions constrained in a subset Ω of
R2 is shown in Figure 1 (right), where similarly the cell
boundaries are shown in blue.
Remark 1. In contrast to the classic Voronoi diagram, the
GuaranteedSVoronoi diagram does not constitute a tessellation
of Ω,Ssince i∈In Vig ⊂ Ω. We define a neutral region O so that
O ∪ i∈In Vig = Ω.
Similarly to the Voronoi diagram, we can define the Guaranteed Delaunay neighbors Nig of a region Di so that in order

IEEE TRANSACTIONS ON AUTOMATIC CONTROL

3

to construct Vig , only the regions in Nig need be considered.
Thus the Guaranteed Delaunay neighbors are defined as

(6)
Nig = j ∈ In , j 6= i : ∂ Hi j ∩ ∂Vig 6= 0/ .
Having defined the Guaranteed Delaunay neighbors, the
definition of the GV diagram can be simplified into
Vig =

[

Hi j .

(7)

(a)

(b)

(c)

(d)

i∈Ni

It is important to note that when two or more uncertain
regions overlap, none of them are assigned GV cells. In
addition, if the uncertain regions degenerate into points, the
GV diagram converges to the classic Voronoi diagram.
C. Guaranteed Voronoi Diagram of Disks
The GV diagram is examined in the case when the uncertain
regions Di are disks Ciu as defined in Equation (2). In that case
it has been shown that Hi j are regions on the plane bounded
by hyperbolic branches. For any two nodes i and j, ∂ Hi j and
∂ H ji are the two branches of a hyperbola with foci at qi and
q j and semi-major axis equal to the sum of the radii of Ciu
and Cuj . The parametric equation of these hyperbola branches
can be found in the Appendix.
The dependence of the GV cells of two disks on the distance
of their centers di j = d(qi , q j ) is as follows. When the disks
Ciu ,Cuj , overlap both of the cells Vig ,V jg are empty as shown
in Fig. 2 (a). If the disks represent the possible positions of
the agents, this is an undesirable configuration as it can lead
to collisions. When the disks Ciu ,Cuj are outside tangent (i.e.,
d(qi , q j ) = riu + ruj ), the resulting cells are rays starting from
the centers of the disks qi , q j and extending along the direction
of qi q j as seen in Fig 2 (b). When the disks are disjoint, the
GV cells are bounded by the two branches of a hyperbola. As
di j increases further, the eccentricity of the hyperbola increases
and the distance of the disk centers from the hyperbola vertices
(the points of the hyperbola closest to its center) increases as
seen in Fig 2 (c), (d). The distance between the hyperbola’s
vertices remains constant at ri + r j as di j increases.
The GV cells of two disks also depend on the sum of
their radii ri + r j as seen in Figure 3. As ri + r j decreases,
the eccentricity of the hyperbola increases and the result on
the cells is the same as if the distance of the disks’ centers
increased, as explained previously. When ri + r j = 0, the GV
cells are the classic Voronoi cells.
It has been also shown that the Guaranteed Delaunay
neighbors Nig are a subset of the classic Delaunay neighbors
of the uncertain disks’ centers. The two sets Ni and Nig are
equal when the GV diagram is constructed for the whole plane,
however when the GV diagram is constrained within a subset
Ω of the plane, Nig can be a subset of Ni . Such a case is shown
in Figure 1 where nodes 1 and 5 are Delaunay neighbors
but not Guaranteed Delaunay neighbors. This way, the GV
diagram of disks can be constructed in a manner similar to
the classic Voronoi diagram. First the Delaunay triangulation
of the disks’ centers is constructed and then each cell Vig is
constructed by intersecting the regions bound by hyperbolic
arcs Hi j ∀ j ∈ Nig . This process is shown in Figure 4.

Fig. 2. Dependence of the GV cells on qi − q j .

Fig. 3. Dependence of the GV cells on ri + r j .

It is useful to define Guaranteed-Sensing Voronoi (GSV)
cells as
Vigs = Vig ∩Cigs , i ∈ In ,
(8)
which are the parts of the GV cells that are guaranteed to be
sensed by their respective nodes. It can be shown that since
GV cells are disjoint, GSV cells are also disjoint regions.
IV. P ROBLEM S TATEMENT - D ISTRIBUTED C ONTROL L AW
By taking into account the GV partitioning of the space, the
following area coverage criterion is constructed
H =

Z

Z

S

φ (q) dq =
gs

i∈In Vi

∑

gs

i∈In Vi

φ (q) dq =

∑ Hi ,

(9)

i∈In

Fig. 4. The construction of a GV cell by the hyperbolic region intersection
method.

IEEE TRANSACTIONS ON AUTOMATIC CONTROL

4

where the function φ : R2 → R+ is related to the a priori
knowledge of importance of a point q ∈ Ω indicating the
probability of an event to take place at q in a coverage
scenario. This criterion expresses the area that is guaranteed
to be closest to and covered by the nodes.
A distributed control law is designed in order to maximize
the coverage criterion (9) while taking into account the nodes’
guaranteed sensing pattern (3) and the GV partitioning (7). It
is shown that this control law leads to a monotonic increase
of the coverage criterion.

The boundary ∂ Ṽigs can be decomposed in disjoint sets as



∂Vigs = ∂Vigs ∩ ∂ Ω ∪ ∂Vigs ∩ ∂Cigs ∪ ∂Vigs ∩ ∂ O . (14)
These sets represent the parts of ∂Vigs that lie on the boundary
of Ω, the boundary of the node’s guaranteed sensing region,
and the boundary of the unassigned neutral region of Ω,
respectively. This decomposition is shown in Figure 5 with
∂Vigs ∩ ∂ Ω in green, ∂Vigs ∩ ∂Cigs in solid red and ∂Vigs ∩ ∂ O
in solid blue.

Theorem 1. For a network of mobile nodes with uncertain
localization as in (2), uniform range–limited radial sensing
performance as in (3), governed by the individual agent’s
kinodynamics described in (1), and the GV–partitioning of Ω
defined in (7), the coordination scheme
Z

ui = αi
"Z
αi

∑g

j∈Ni

gs

gs

gs

∂Vi ∩∂Ci

∂Vi ∩∂ Hi j

υii ni φ

ni φ dq +
#

Z

dq +

gs

∂V j ∩∂ H ji

υ ij n j φ

dq

(10)

where ni is the outward unit normal on ∂Vigs and αi a positive
constant, maximizes the performance criterion (9) along the
nodes’ trajectories in a monotonic manner, leading in a locally
area–optimal configuration of the network.

Fig. 5. The decomposition of ∂Vigs into disjoint sets (solid red, green and
blue).

Proof. In order to guarantee monotonous increase of the
coverage criterion (9), we evaluate its time derivative

Hence

∂ Hi
∂ Hi
∂H
=∑
q̇i = ∑
ui
∂t
∂
q
i
i∈In
i∈In ∂ qi

∂H
∂ qi

(11)

Z

gs

Vi

∂

∑g ∂ qi

j∈Ni

Z
gs

φ dq.

Vj

Considering the second summation term, infinitesimal motion
of qi may only affect ∂V jgs at ∂V jgs ∩ ∂ O, where O is the
neutral region defined in Remark 1, since both hyperbola
branches are affected by alteration of one of the foci. In
addition, only the Guaranteed Delaunay neighbors (6) of i are
considered in the summation, as a major property of GV–
partitioning. Therefore, the former expression can be written
via the generalized Leibniz integral rule [25] (by converting
surface integrals to line ones) as
∂H
=
∂ qi

Z

υ i n φ dq +
gs i i

∂Vi

Z

∑g

j∈Ni

gs

∂Vi ∩∂ O

(12)

monotonous increase of H is achieved.
We now evaluate the partial derivative of H with respect
to qi as
φ dq +
gs

gs

∂Vi ∩∂ Ω

Z

∂ Hi
ui =
, i ∈ In ,
∂ qi

Z

can be written as

=

By using the gradient–based control law

∂
∂H
=
∂ qi
∂ qi

∂H
∂ qi

∂V j ∩∂ O

υ ij n j φ dq,

where υii , υ ij stand for the transpose Jacobian matrices with
respect to qi of the points q ∈ ∂Vigs , q ∈ ∂V jgs , respectively,
i.e.
T
4 ∂q
, q ∈ ∂V jgs , i, j ∈ In .
(13)
υ ij (q) =
∂ qi

υii ni φ dq +

υii ni φ dq +

Z
gs

gs

∂Vi ∩∂Ci

Z

∑g

j∈Ni

gs

∂V j ∩∂ O

υii ni φ dq +

υ ij n j φ dq.

(15)

It is apparent that υii = 0 at q ∈ ∂Vigs ∩ Ω since all q ∈ ∂ Ω
remain unaltered by infinitesimal motions of qi , assuming no
alteration of the environment. Considering the second integral,
for any point q ∈ ∂ Ṽigs ∩ ∂Cigs it holds that υii (q) = I2 , where
I stands for the identity matrix, since they translate along the
direction of motion of qi at the same rate. As far as the sets
∂Vigs ∩ ∂ O and ∂V jgs ∩ ∂ O, j ∈ Nig are concerned, they can be
merged in pairs via utilization of the left and right hyperbolic
branches, as introduced in GV–partitioning (7), as follows
Z
gs

∂Vi ∩∂ O

υii ni φ

"Z

∑g

j∈Ni

gs

∂Vi ∩∂ Hi j

Z

dq +

∑g

gs

∂V j ∩∂ O

j∈Ni

υii ni φ

#

Z

dq +

υ ij n j φ dq =

gs

∂V j ∩∂ H ji

υ ij n j φ

dq

.

(16)

However, for any two Delaunay neighbors i, j it can be
observed that
• The hyperbolic branches ∂ Hi j and ∂ H ji are symmetric
with respect to the perpendicular bisector of qi , q j , and
are governed by the same set of parametric equations (left
and right branch).
• The vectors n j are mirrored images of the corresponding
ni with respect to the perpendicular bisector of qi , q j .

IEEE TRANSACTIONS ON AUTOMATIC CONTROL

∂H
∂ qi

Taking into account the above,
∂H
∂ qi
"Z

∑g

j∈Ni

5

can be written as

A. Constraining Nodes Inside Ω

Z

=

gs

gs
gs
∂Vi ∩∂Ci

∂Vi ∩∂ Hi j

υii ni φ

ni φ dq +
#

Z

dq +

gs

∂V j ∩∂ H ji

υ ij n j φ

dq

.

(17)

The unit normal vectors ni , n j and Jacobian matrices υii , υ ij
in the second part of (17) can be evaluated via utilization
of the parametric representations of the sets over which the
integration takes place. These sets are parts of the left (∂ Hi j )
and right (∂ H ji ) branch of the hyperbola that assigns the space
among two arbitrary nodes i, j. The second term
of the sum
S
in (17) requires knowledge of all the nodes in j∈N g N gj ⊆ In
i
and so the control law is distributed.
The decomposition of the set ∂Vigs ∩ ∂ O of node i into
mutually disjoint hyperbolic arcs ∂Vigs ∩ ∂ Hi j and ∂Vigs ∩ ∂ Hik
is shown in Figure 6. Figure 6 also shows the other domains of
integration used in (17) ∂Vigs ∩ ∂Cigs , ∂V jgs ∩ ∂ H ji and ∂Vkgs ∩
∂ Hki .

When using this control law, it is possible that a portion
of some node’s positioning uncertainty region lies outside the
region Ω, i.e. that Ciu ∩ Ω 6= Ciu for some node i. As a result,
it would be possible for that node to be outside the region Ω
given that it may be anywhere within Ciu .
In order to prevent this situation, instead of Ω, a subset
Ωsi ⊆ Ω will be used, in the general case different for each
node since their uncertainty radii are allowed to differ. This
subset is calculated as the Minkowski difference of Ω with the
disk C (riu ) = {q ∈ Ω : k q k≤ riu }.
Ωsi = {q ∈ Ω | q +C(riu ) ⊆ Ω} , i ∈ In .
Thus by constraining the center of each uncertainty disk
qi inside Ωsi , it is guaranteed that no part of the uncertainty
disk will ever be outside Ω. This is achieved by projecting
the velocity control input ui on ∂ Ωsi . This is only needed if
qi ∈ ∂ Ωsi and the contorl input ui points towards the outside
of Ωsi . As a result, the control law becomes
(
proj∂ Ωs ui ifqi ∈ ∂ Ωsi ∧ qi + ε ui ∈
/ Ωsi
i
ũi =
(18)
ui
otherwise
where ε is an infinitesimally small positive constant.
B. Suboptimal Control Law
Because of the high complexity of the integrals over Hi j and
H ji of the control law (10), the large amount of information
required in order to implement it, as described in Remark 2, as
well as the problems described in Section IV-A, a simplified
version of the control law is proposed.

Fig. 6. The domains of integration of the control law (10).

The proposed law (10) leads to a gradient flow of H along
the nodes trajectories, while H increases monotonically, since
dH
∂H
=∑
· q̇i = ∑ kui k2 ≥ 0.
dt
i∈In ∂ qi
i∈In

Conjecture 1. The computationally efficient suboptimal control law
Z
ui = αi
n φ dq
(19)
gs
gs i
∂Vi ∩∂Ci

leads the nodes to suboptimal trajectories.
The control law enhancement described in Section IV-A is
used in conjunction with this simplified control law.
V. S IMULATION S TUDIES

Remark 2. It is not immediately clear from the control law
(10) the information node i requires in order to implement
it. The integrals over ∂Vigs ∩ ∂Cigs and ∂Vigs ∩ ∂ Hi j require
knowledge of the Guaranteed Delaunay neighbors Nig in order
to compute Vigs and Hi j . However the integral over ∂V jgs ∩
∂ H ji requires knowledge of the whole cell V jgs of node i and
as such requires informations form the Guaranteed Delaunay
neighbors N gj of node j. Thus node i requires information from
all nodes in
[ g
Ñig =
Nj ,

Simulations have been conducted in order to evaluate the
efficiency of the control law as well as compare the complete
(10) with the simplified (19) law. All nodes have common
positioning uncertainty regions Ci , ∀qi ∈ Ω as well as sensing
performance Cigs . The a priori importance of points inside the
region was φ (q) = 1, ∀q ∈ Ω. In this case, the maximum
possible value for H is

H max = ∑ A Cigs ,
i∈In

g

j∈Ni

which are the Guaranteed Delaunay neighbors of all the
Guaranteed Delaunay neighbors of node i.

where A (·) is the area function, given that n guaranteed
sensing disks can be packed inside Ω. In all simulations, the
control law extension as described in Section IV-A is used.

IEEE TRANSACTIONS ON AUTOMATIC CONTROL

6

A. Case Study I
In this case study the convergence speed between the
optimal and suboptimal laws is examined. A network of three
nodes in a large region is simulated. Figure 7 shows the initial
and final states of the network, with the uncertainty disks
shown in black, the guaranteed sensing disks in red and the
GV cells in blue. In Figure 8 the trajectories of the nodes
when using the optimal control law are shown in red and
when using the suboptimal control law in blue, while the initial
positions are represented by circles and the final by squares.
It is observed that the final configurations in both cases are
very similar since the average distance between nodes’ final
states is 0.57% of the region’s diameter. This is due to the low
number of nodes in the network and is a result that can not
be generalized, as seen in the following case study. Figure 9
shows the evolution of the coverage objective H with time,
again shown in red for the optimal control law and in blue for
the suboptimal. As it is expected, the optimal law converges
faster than the suboptimal. The final value of the objective
function is the same for both control laws in this case since
the network reaches a globally optimal state in both cases.

Fig. 8. The trajectories of the nodes using the optimal (red) and suboptimal
(blue) control laws.

(a)
Fig. 9. The coverage objective H using the optimal (red) and suboptimal
(blue) control laws.

(b)

(c)

Fig. 7. Case Study I: (a) Initial state, (b) Final state using the optimal control
law, (c) Final state using the suboptimal control law.

B. Case Study II
In this case study network of ten nodes is examined. Figure
10 shows the initial and final states of the network, with the
uncertainty disks shown in black, the guaranteed sensing disks
in red and the GV cells in blue. In Figure 11 the trajectories
of the nodes when using the optimal control law are shown in
red and when using the suboptimal control law in blue, while
the initial positions are represented by circles and the final
by squares. It is observed that the final configurations in both
cases seem very similar at first glance, however the average
distance between nodes’ final states is 10.12% of the region’s

diameter. Careful examination of the trajectories shows that
not all nodes starting from the same initial position converge to
the same final position, explaining the large average distance.
. Figure 12 shows the evolution of the coverage objective H
with time, again shown in red for the optimal control law and
in blue for the suboptimal. As it is expected, the optimal law
converges faster than the suboptimal. The final value of the
objective function is the same for both control laws in this case
because of the similarity of their final configurations. However,
since both control laws converge to local maxima, there are
no guarantees that one or the other will always achieve better
performance.
This case study highlights the need for the control enhancement in Section IV-A, especially in the case of the
optimal control law. It can be seen in Figure 11 that a node
moving under the optimal law (red trajectory) was stopped
from moving outside the region on the lower left corner.
VI. E XPERIMENTAL S TUDIES
In order to evaluate the simplified control law (19), two
experiments were planned and conducted. The experiments

IEEE TRANSACTIONS ON AUTOMATIC CONTROL

7

(a)

Fig. 12. The coverage objective H using the optimal (red) and suboptimal
(blue) control laws.

(b)

(c)

Fig. 10. Case Study I: (a) Initial state, (b) Final state using the optimal control
law, (c) Final state using the suboptimal control law.

pose estimation server. ArUco [26] is a minimal library for
Augmented Reality applications and provides pose estimation
of some predefined sets of fiducial markers. Hende, by attaching a fiducial marker on top of each robot, the ArUco library
provides an estimation of each robot’s pose.
In order to estimate the positioning uncertainty of the
experimental setup, a fiducial marker was placed on the
centroid and on each of the vertices of the region Ω. A
3D triangular mesh was created from these measurements as
shown in Figure 13. The maximum uncertainty value from that
mesh rumax = 0.032 m was used for the experiments.

Fig. 11. The trajectories of the nodes using the optimal (red) and suboptimal
(blue) control laws.

consist of 3 robots, an external pose tracking system and
a computer, communicating with the robots and the pose
tracking system via a WiFi router.
The robots used are the differential-drive AmigoBots by
ActiveMedia Robotics, which provide a high level command
set, such as setting their translational and rotational velocities
directly. These robots are also equipped with encoders, used
for self-estimating their pose, but due to the drifting error,
an external tracking system was implemented. Although the
external pose estimation system is used for the control law
implementation, its performance is compared with that of
the encoders. Since the robots are not equipped with actual
sensors, circular sensing patters with a radius rs = 0.3 m were
assumed in the control law implementation.
The pose estimation system consists of a camera and an
ODROID-XU4 octa core computer, running an ArUco based

Fig. 13. Triangular mesh from positioning uncertainty measurements.

The control law was implemented algorithmically as a loop
with a period Ts = 0.1 sec. After the end of each iteration the
computer sends velocity commands to the robots. At first, the
robots’ current positions are used to create the Guaranteed
Voronoi diagram of the uncertain disks and calculate the
control inputs according to the simplified control law (19).
Subsequently, a target point qti ∈ R2 is created for each robot
based on its current position qi and its control input ũi . Once

IEEE TRANSACTIONS ON AUTOMATIC CONTROL

8

the target point for each robot has been found, and given each
robot’s position qi and orientation θi , the translational vi and
rotational ωi velocity control inputs are sent to each robot.
 t

kqi − qi k
vi = min
, vmax cos(dθi )
Ts


|dθi |
ω = min
, ωmax sin(dθi )
Ts
where dθi = ∠(qti − qi ) − θi and vmax , ωmax are constraints on
the maximum translational and rotational velocity respectively
of the Amigobot robots. Once all robots are within a predefined distance dt = 0.02 m of their respective target points, the
Guaranteed Voronoi diagram and the control law are calculated
again and new target points are produced.
The region shown in all figures regarding the experiments is
Ωs ⊆ Ω i.e. the region inside which the centers of the uncertain
disks are constrained, as explained in Section IV-A.
Both experiments are evaluated against simulations with
the same initial robot positions, positioning uncertainty and
sensing performance.
A. Experiment I
The initial and final positions of the robots for both the
experiment and the simulation are shown in Figures 14 and
15 respectively with the uncertainty disks shown in black, the
guaranteed sensing disks in red and the GV cells in blue.
Additionally, photos of the experiment initial and final configurations can be seen on Figure 16. It can be seen that in the
final configuration all guaranteed sensing disks are contained
within their respective GV cells in both the experiment and
the simulation, thus leading the network to a globally optimal
configuration. By comparing the final configurations, as well
as the robot trajectories for the experiment and simulation
shown in Figure 17, it can be seen that they differ significantly.
The mean distance between the nodes final positions in the
experiment and simulation is 8.68% of the diameter of Ωs .
Because in this current case there are multiple globally optimal
configurations the network can possibly reach and since the
implemented control law differs from the theoretical one, it
is natural that the final configurations between the experiment
and the simulation differ. Nevertheless, a globally optimal configuration was reached in both cases. The coverage objective
H did not increase monotonously due to the implemented
control law. It did increase however from 83.7% to 99.9% of
its maximum possible value. The positioning data from the
ArUco library and the robot encoders are compared in Figure
19 [Left]. As it is expected, the further a robot moves and the
more it rotates, the larger the positioning error of the encoders
grows. The mean error between the ArUco and encoder final
robot positions is 4.04% of the diameter of Ωs . Figure 19
[Right] shows the trajectories of the target points used in the
control law implementation.

Fig. 14. Experiment I: Initial configuration.

Fig. 15. Experiment I: Experiment [Left] and simulation [Right] final
configuration.

Ciu . The guaranteed sensing disks however were calculated
using only the actual uncertainty and not the one increased
by the robot radius. This was done to guarantee that the
robots would remain inside the region Ω. The initial and
final positions of the robots for both the experiment and the
simulation are shown in Figures 20 and 21 respectively with
the uncertainty disks shown in black, the guaranteed sensing
disks in red and the GV cells in blue. Additionally, photos
of the experiment initial and final configurations can be seen
on Figure 22. It can be seen that in the final configuration
all guaranteed sensing disks are almost completely contained
within their respective GV cells in both the experiment and the

Fig. 16. Experiment I: Photos from the initial [Left] and final [Right] robot
positions in the experiment.

B. Experiment II
In this experiment, half the diameter of the robots was
added to their positioning uncertainty, so that each robot was
completely contained within its positioning uncertainty disk

Fig. 17. Experiment I: Comparison of experiment (red) and simulation (green)
robot trajectories. Initial positions marked with circles and final positions with
squares.

IEEE TRANSACTIONS ON AUTOMATIC CONTROL

9

Fig. 20. Experiment II: Initial configuration.

Fig. 18. Experiment I: Coverage objective H as a function of time.

Fig. 21. Experiment II: Experiment [Left] and simulation [Right] final
configuration.

Fig. 19. Experiment I: [Left] ArUco (green) and encoder (red) reported
trajectories. [Right] Target point trajectories. Initial positions marked with
circles and final positions with squares.

simulation. By comparing the final configurations, as well as
the robot trajectories for the experiment and simulation shown
in Figure 23, it can be seen that they do not differ significantly.
The mean distance between the nodes final positions in the
experiment and simulation is 1.54% of the diameter of Ωs .
Because of the increased positioning uncertainty in this case,
which results in smaller GV cells, the number of local maxima
of the objective function is significantly smaller than in the
previous experiment. Thus it is quite possible that, despite their
differences, the theoretical and implemented control laws converge to the same final configuration. The coverage objective
H did not increase monotonously due to the implemented
control law. It did increase however from 38.8% to 98.3% of
its maximum possible value. The positioning data from the
ArUco library and the robot encoders are compared in Figure
25 [Left]. As it is expected, the further a robot moves and the
more it rotates, the larger the positioning error of the encoders
grows. The mean error between the ArUco and encoder final
robot positions is 1.54% of the diameter of Ωs . Figure 25
[Right] shows the trajectories of the target points used in the
control law implementation.
A video of the simulations and the experiments can be
found on http://anemos.ece.upatras.gr/images/stories/videos/
PTGS IEEETAC.mp4
VII. C ONCLUSION
This article examines the area coverage problem by a homogeneous team of mobile agents with imprecise localization.

A gradient ascent based control law was designed based on
a Guaranteed Voronoi partition of the region. Because of the
complexity of the resulting control law, a simpler suboptimal
control law is proposed and the performance of both is
compared in simulation studies. Additionally, two experiments
were conducted to highlight the efficiency of the suboptimal
control law.
A PPENDIX
We assume n disks Di , i ∈ In with centers qi = [xi , yi ]T ∈ R2
and radii riu , i ∈ In .

Fig. 22. Experiment II: Photos from the initial [Left] and final [Right] robot
positions in the experiment.

Fig. 23. Experiment II: Comparison of experiment (red) and simulation
(green) robot trajectories. Initial positions marked with circles and final
positions with squares.

IEEE TRANSACTIONS ON AUTOMATIC CONTROL

10

Fig. 24. Experiment II: Coverage objective H as a function of time.
Fig. 26. A hyperbola with its parameters 2a and 2c shown.

B. Hyperbola Parametric Equation
Assuming the foci of the hyperbola are qi = [−c 0]T and
q j = [c 0]T , its parametric equation is


±a cosh(t)
γ(t) =
, t ∈R
b sinh(t)
Fig. 25. Experiment II: [Left] ArUco (green) and encoder (red) reported
trajectories. [Right] Target point trajectories. Initial positions marked with
circles and final positions with squares.

A. Hyperbola Definition and Properties
We will use the GV diagram definition given in (5) for two
disks Di and D j and shown that the boundaries of their GV
cells are hyperbola branches.
Let us find the boundary of the cell of node i, Vig = Hi j ,
since we examine only two nodes. We have that

∂ Hi j = q ∈ Ω : max kq − bi k = min q − b j , ∀bi ∈ Di , ∀b j ∈ D j ,

however since Di and D j are disks, the maximum and minimum distances can be calculated and thus

∂ Hi j = q ∈ Ω : kq − qi k + riu = q − q j − riu ,

where qi and q j are the centers of Di and D j respectively. The
equation

∂ Hi j = q ∈ Ω : q − q j − kq − qi k = ri + r j ,
defines one branch of a hyperbola since it is the locus of points
whose difference of distances from two foci qi and q j is equal
to a positive constant ri +r j . The other branch of the hyperbola
is

∂ H ji = q ∈ Ω : kq − qi k − q − q j = ri + r j ,
and corresponds to the boundary of V jg .
In a hyperbola, the positive constant is called the major axis
and is usually denoted as 2a whereas the distance between
its foci is denoted as 2c. A third parameter, the minor axis
denoted as 2b is then defined by the equation a2 + b2 = c2 . As
such in this case we have that 2a = ri + r j and 2c =k qi − q j k.
Additionally, the closest points on the two branches are called
the vertices of the hyperbola and their distance is 2a.

where a, b and c are those defined in Appendix A and the
+(−) sign corresponds to the West (East) branch of the
hyperbola which is H ji (Hi j ).
By rotating the hyperbola foci around the origin by an angle
θ , the parametric equation becomes



cos(θ ) − sin(θ )
±a cosh(t)
γ(t) =
, t ∈R
sin(θ ) cos(θ )
b sinh(t)
Finally, given two foci qi = [xi yi ]T and q j = [x j y j ]T
anywhere on the plane, the general hyperbola parametric
equation is


 " xi +x j #
cos(θ ) − sin(θ )
±a cosh(t)
2
, t ∈R
γ(t) =
+ yi +y
j
sin(θ ) cos(θ )
b sinh(t)
2
(20)
y −y
where θ = arctan( x jj −xii ).
C. Outward Unit Normal Vectors
For any curve γ(t) the normal vector is defined as


γ̇
γ̇
n̂(t) = γ̈ − γ̈ ·
.
k γ̇ k k γ̇ k
The magnitude of this vector is the curvature at that particular
point on the curve and it points towards the center of curvature.
Since the hyperbolic branches define convex regions on the
plane, the outward unit normal vector is given by
n(t) = −

n̂
.
k n̂ k

(21)

However, since the resulting expressions are too large to show
here, an indicative plot of the outward unit normal vectors of
one branch of a hyperbola is shown in Figure 27.

IEEE TRANSACTIONS ON AUTOMATIC CONTROL

11

R EFERENCES

Fig. 27. The outward unit normal vectors of one branch of the hyperbola are
shown in black.

D. Jacobian Matrix
T

i
The Jacobian

matrix υ j shows the change in the curve
γ jx (t)
γ j (t) =
caused by the movement of node i, thus its
γ jy (t)
transpose is
" ∂γ
#
∂γ

υ ij =

jx

jx

∂ xi
∂ γ jy
∂ xi

∂ yi
∂ γ jy
∂ yi

(22)

The expressions for the Jacobian matrix elements are too large
to show here and as such some indicative plots of them with
respect to the parameter t are shown in Figure 28.

(a)

(b)

(c)

(d)

Fig. 28. The Jacobian matrix
t. (a)

∂ γ jx
∂ xi ,

(b)

∂ γ jx
∂ yi ,

(c)

∂ γ jy
∂ xi ,

υii
(d)

(blue) and
∂ γ jy
∂ yi .

υ ij

(red) elements as functions of

[1] H. Mahboubi and A. G. Aghdam, “Distributed deployment algorithms
for coverage improvement in a network of wireless mobile sensors:
Relocation by virtual force,” IEEE Transactions on Control of Network
Systems, vol. 61, no. 99, p. 1, 2016.
[2] C. F. van Eeden, G. Lee, S. Katsuno, and H. Yonekura, “Minimised
moving distance algorithm for geometric deployment of robot swarms,”
in Proc. 55th Annual Conf. of the Society of Instrument and Control
Engineers of Japan (SICE), Sep. 2016, pp. 1195–1200.
[3] A. Simonetto, T. Keviczky, and D. V. Dimarogonas, “Distributed solution
for a maximum variance unfolding problem with sensor and robotic
network applications,” in Proc. and Computing (Allerton) 2012 50th
Annual Allerton Conf. Communication, Control, Oct. 2012, pp. 63–70.
[4] C. Song, L. Liu, G. Feng, Y. Wang, and Q. Gao, “Persistent awareness
coverage control for mobile sensor networks,” Automatica, vol. 49, no. 6,
pp. 1867–1873, 2013.
[5] C. Franco, D. M. Stipanović, G. López-Nicolás, C. Sagüés, and
S. Llorente, “Persistent coverage control for a team of agents with
collision avoidance,” European Journal of Control, vol. 22, pp. 30–45,
2015.
[6] C. Zhai and Y. Hong, “Decentralized sweep coverage algorithm for
multi-agent systems with workload uncertainties,” Automatica, vol. 49,
no. 7, pp. 2154–2159, 2013.
[7] Y. Stergiopoulos, M. Thanou, and A. Tzes, “Distributed collaborative
coverage-control schemes for non-convex domains,” IEEE Transactions
on Automatic Control, vol. 60, no. 9, pp. 2422–2427, 2015.
[8] R. J. Alitappeh, K. Jeddisaravi, and F. G. Guimarães, “Multi-objective
multi-robot deployment in a dynamic environment,” Soft Computing, pp.
1–17, 2016.
[9] Y. Stergiopoulos and A. Tzes, “Convex Voronoi-inspired space partitioning for heterogeneous networks: A coverage-oriented approach,” IET
Control Theory and Applications, vol. 4, no. 12, pp. 2802–2812, 2010.
[10] O. Arslan and D. E. Koditschek, “Voronoi-based coverage control of heterogeneous disk-shaped robots,” in 2016 IEEE International Conference
on Robotics and Automation (ICRA). IEEE, 2016, pp. 4259–4266.
[11] L. Pimenta, V. Kumar, R. Mesquita, and G. Pereira, “Sensing and
coverage for a network of heterogeneous robots,” in Proc. of the 47th
Conference on Decision and Control, Cancun, Mexico, 2008, pp. 3947–
3952.
[12] Y. Stergiopoulos and A. Tzes, “Cooperative positioning/orientation
control of mobile heterogeneous anisotropic sensor networks for area
coverage,” in Proc. IEEE International Conference on Robotics &
Automation, Hong Kong, China, 2014, pp. 1106–1111.
[13] J. M. Luna, R. Fierro, C. Abdallah, and J. Wood, “An adaptive coverage
control algorithm for deployment of nonholonomic mobile sensors,” in
Proc. 49th IEEE Conf. Decision and Control (CDC), Dec. 2010, pp.
1250–1256.
[14] J. Cortés and F. Bullo, “Coordination and geometric optimization via
distributed dynamical systems,” SIAM Journal on Control and Optimization, vol. 44, no. 5, pp. 1543–1574, 2005.
[15] J. Cortés, S. Martinez, and F. Bullo, “Spatially-distributed coverage optimization and control with limited-range interactions,” ESAIM: Control,
Optimisation and Calculus of Variations, vol. 11, no. 4, pp. 691–719,
2005.
[16] M. T. Nguyen and C. S. Maniu, “Voronoi based decentralized coverage
problem: From optimal control to model predictive control,” in Proc.
24th Mediterranean Conf. Control and Automation (MED), Jun. 2016,
pp. 1307–1312.
[17] F. Mohseni, A. Doustmohammadi, and M. B. Menhaj, “Distributed
receding horizon coverage control for multiple mobile robots,” IEEE
Systems Journal, vol. 10, no. 1, pp. 198–207, Mar. 2016.
[18] V. Ramaswamy and J. R. Marden, “A sensor coverage game with improved efficiency guarantees,” in Proc. American Control Conf. (ACC),
Jul. 2016, pp. 6399–6404.
[19] M. T. Nguyen, L. Rodrigues, C. S. Maniu, and S. Olaru, “Discretized
optimal control approach for dynamic multi-agent decentralized coverage,” in Proc. IEEE Int. Symp. Intelligent Control (ISIC), Sep. 2016, pp.
1–6.
[20] J. A. Hage, M. E. E. Najjar, and D. Pomorski, “Fault tolerant collaborative localization for multi-robot system,” in Proc. 24th Mediterranean
Conf. Control and Automation (MED), Jun. 2016, pp. 907–913.
[21] B. Davis, I. Karamouzas, and S. J. Guy, “C-opt: Coverage-aware trajectory optimization under uncertainty,” IEEE Robotics and Automation
Letters, vol. 1, no. 2, pp. 1020–1027, Jul. 2016.

IEEE TRANSACTIONS ON AUTOMATIC CONTROL

[22] J. Habibi, H. Mahboubi, and A. G. Aghdam, “Distributed coverage
control of mobile sensor networks subject to measurement error,” IEEE
Transactions on Automatic Control, vol. 61, no. 11, pp. 3330–3343, Nov.
2016.
[23] W. Evans and J. Sember, “Guaranteed Voronoi diagrams of uncertain
sites,” in 20th Canadian Conference on Computational Geometry, 2008,
pp. 207–210.
[24] S. Papatheodorou, Y. Stergiopoulos, and A. Tzes, “Distributed area
coverage control with imprecise robot localization,” in Proc. 24th
Mediterranean Conf. Control and Automation (MED), Jun. 2016, pp.
214–219.
[25] H. Flanders, “Differentiation under the integral sign,” American Mathematical Monthly, vol. 80, no. 6, pp. 615–627, 1973.
[26] S. Garrido-Jurado, R. M. noz Salinas, F. Madrid-Cuevas, and M. Marı́nJiménez, “Automatic generation and detection of highly reliable fiducial
markers under occlusion,” Pattern Recognition, vol. 47, no. 6, pp. 2280
– 2292, 2014.

12

