IEEE TRANSACTIONS ON
1
Motion Planning and Collision Avoidance
using Non-Gradient Vector Fields
Dimitra Panagou
Abstract
This paper presents a novel feedback method on the motion planning for unicycle robots in
environments with static obstacles, along with an extension to the distributed planning and coordination
in multi-robot systems. The method employs a family of 2-dimensional analytic vector fields, whose
integral curves exhibit various patterns depending on the value of a parameter λ. More specifically, for
an a priori known value of λ, the vector field has a unique singular point of dipole type and can be
used to steer the unicycle to a goal configuration. Furthermore, for the unique value of λ that the vector
field has a continuum of singular points, the integral curves are used to define flows around obstacles.
An almost global feedback motion plan can then be constructed by suitably blending attractive and
repulsive vector fields in a static obstacle environment. The method does not suffer from the appearance
of sinks (stable nodes) away from goal point. Compared to other similar methods which are free of local
minima, the proposed approach does not require any parameter tuning to render the desired convergence
properties. The paper also addresses the extension of the method to the distributed coordination and
control of multiple robots, where each robot needs to navigate to a goal configuration while avoiding
collisions with the remaining robots, and while using local information only. More specifically, based
on the results which apply to the single-robot case, a motion coordination protocol is presented which
guarantees the safety of the multi-robot system and the almost global convergence of the robots to their
goal configurations. The efficacy of the proposed methodology is demonstrated via simulation results
in static and dynamic environments.
I. INTRODUCTION
Motion planning, coordination and control for robotic systems still remains an active research
topic in many respects. The primary motivation has been the computation of safe, collision-free
Dimitra Panagou is with the Department of Aerospace Engineering, University of Michigan, Ann Arbor, MI, USA;
dpanagou@umich.edu.
October 22, 2014
DRAFT
arXiv:1404.5828v3  [cs.RO]  13 Oct 2014
IEEE TRANSACTIONS ON
2
trajectories for robotic agents, mechanisms and autonomous vehicles which operate in constrained
and/or uncertain environments. Research within the robotics community has attributed various
formulations and methodologies on the motion planning problem, often specialized based on the
control objectives and the characteristics of the problems at hand. These methodologies range
from Lyapunov-based control methods, to sampling-based planning, to combinatorial planning,
to formal methods [1]–[3]. Multi-robot systems have attracted the interest of the control systems
community as well. Emphasis has been given in consensus, flocking and formation control
problems for multiple agents [4].
Avoiding obstacles and inter-agent collisions is a requirement of highest priority in motion
planning and coordination problems. Recently, significant interest has been paid to the high-
level task planning under complex goals, where the problem for an autonomous robot has
transitioned from the classical motion planning formulation (i.e., move from point A to point
B) to the consideration of complex goals under temporal specifications; such specifications are
typically described as: “visit region A, and then visit either region B or region C”. Despite the
tremendous and elegant contributions in this area, which provide elegant solutions to the high-
level mission synthesis with rigorous guarantees under certain assumptions on the considered
environments [2], [3], the interconnection of high-level tasking with the physical layer/system
is still an open problem in many respects. One issue is the consideration of multiple agents in
dynamic environments and the associated complexity in finding provably correct solutions in the
presence of nonlinearities, arbitrary constraints, and uncertainty.
The scope of this paper is to provide a solution to the motion planning problem for single
and multiple nonholonomic agents in dynamic environments, where agents have local sensing
and communication capabilities and which may be populated by dynamic (moving) obstacles.
Our goal is to provide a feedback synthesis of low-level planning controllers along with certain
guarantees, which can later on be combined with high-level tasks, such as dynamic coverage
[5], towards provably correct feedback solutions for a specific class of dynamical systems in
dynamic environments. The technical tools which we use towards this goal are set-invariance
methods, which have been proved efficient in constrained control problems of a class of nonlinear,
under-actuated systems [6].
The spirit of the proposed solutions is similar, but not identical to, Lyapunov-like scalar
functions, such as the Avoidance Functions in [7] and the Artificial Potential Fields (APF) in
October 22, 2014
DRAFT
IEEE TRANSACTIONS ON
3
[8], [9]. More specifically: It is well-known that, although scalar functions offer the merit of
Lyapunov-based control design and analysis, yielding thus solutions in closed-form with certain
guarantees [10], they suffer from the drawback of possible local minima away from the goal
point, i.e., of points in the state space other than the desired equilibrium at which the gradient
vector vanishes; this in principle results in system trajectories which get stuck away from the
goal point. Certain forms of potential functions may overcome this limitation; namely, navigation
functions [11] and harmonic functions [12], [13], but under some cost: the caveat in the former
case is that the Morse property which guarantees the non-existence of local minima is rendered
after a tuning parameter exceeds a lower bound, which is not a priori known. In the latter
case, harmonic functions may be constructed with either discrete or continuous approaches, but
the computational cost of discrete methods is quite demanding. Continuous approaches which
employ the analogies of Laplace equation with fluid mechanics yield closed-form solutions
for certain dynamic environments [14]. Stream functions [15] combine the local-minima-free
property of harmonic functions along with hydrodynamic concepts to yield streamlines which
may be preferable for second order systems. The method of vortex fields [16] uses the anti-
gradient of a scalar function to define flows around obstacles.
Now, let us note that one common ground in this class of solutions is the resulting gradient
vector field which is employed in the control synthesis. In this respect, the idea of directly defining
vector fields encoding obstacle avoidance has been studied for robot motion planning problems.
In [17], for instance, simple smooth vector fields are locally constructed in given convex cell
decompositions of polygonal environments, so that their integral curves are by construction
collision-free and, in a sequential composition spirit, convergent to a goal point. The method,
nevertheless, presumes the existence of a high-level discrete motion plan which determines the
successive order of the cells from an initial to a final configuration. Recent work employing
vector fields for vehicles’ navigation is presented also in [18] and in [19]. The approach with
velocity vector fields in [20] is also relevant to the context. However, these contributions address
only the position control of the robot, while the orientation is not guaranteed to converge to a
desired value.
Stepping now a little further away from single-agent problems: when it comes to multiple
agents, their motion towards goal configurations defines a dynamic environment and poses
challenges to the planning, coordination and control design, even in the absence of static physical
October 22, 2014
DRAFT
IEEE TRANSACTIONS ON
4
obstacles. At the same time, limitations in the available sensing and communication platforms
impose additional constraints to the multi-agent system. Given a pair (i, j) of agents i and j,
agents typically make decisions on their actions based on available information, which can be
either locally measured using onboard sensors, or transmitted and received across the nodes of the
multi-agent system via wireless communication links. Thus, information flow between two agents
can be either bidirectional (undirected) or unidirectional (directed). During the past ten years,
research efforts have achieved the formalization of problems such as consensus and formation
control in multi-agent networks using tools and notions from graph theory, matrix theory and
Lyapunov stability theory [21]–[25]. The case of directed information exchange has recently
attracted increased interest [26]–[30], motivated in part by the fact that undirected information
flow is not always a realistic and practical assumption, due to bandwidth limitations in the
network, anisotropic sensing of the agents etc. Extending consensus algorithms to nonlinear
systems has also become popular, see for instance [31], [32].
Nevertheless, despite that consensus, flocking, and formation control algorithms achieve colli-
sion avoidance in multi-vehicle systems by carefully selecting initial conditions and controlling
relative distance and heading, they are typically not used in encoding problems such as navigation
to specific goal locations for each one of the agents. In this respect, the development of planning
and coordination algorithms for the motion of multiple agents along with safety and performance
guarantees is an open problem in many respects.
A. Overview
This paper presents a novel method on the motion planning and coordination in environments
with static and/or dynamic obstacles, which results in feedback motion plans for unicycle robots
along with collision avoidance guarantees. The method employs a family of two-dimensional
analytic vector fields, originally introduced in [33], given as:
F(r) = λ(pTr)r −p(rTr),
(1)
where λ ∈R is a parameter to be specified later on, r = [x y]T the position vector with respect
to (w.r.t.) a global cartesian frame and p = [px py]T, with p ̸= 0.1
1The role the vector p ∈R2 plays in the properties of the vector field (1) becomes evident later on in Theorem 2.
October 22, 2014
DRAFT
IEEE TRANSACTIONS ON
5
In [33] the family of vector fields (1) was employed in the control design for steering kinematic,
drift-free systems in chained form in obstacle-free environments.
In this paper we first show that, except for a known value of the parameter λ, the vector
field (1) has a unique singular point on R2. More specifically: (i) For λ > 1 the pattern of the
integral curves around the unique singular point is dipolar [34]. Such vector field can be used
for steering a unicycle to a goal configuration. (ii) For λ = 1 the vector field has a continuum
of singular points and can be used to define tangential flows around circular obstacles. (iii) For
λ < 0 the pattern of the integral curves is suitable for defining repulsive flows away from lines,
and as thus, away from polygonal obstacles. A preliminary example is given in the Appendix
of [35].
We then consider the single-agent case in a static environment of circular obstacles and propose
a blending mechanism between attractive and repulsive vector fields, which yields almost global
feedback motion plans. In other words, we construct vector fields whose integral curves are
convergent to a goal configuration, except for a set of initial conditions of Lebesgue measure
zero, and collision-free by construction. This in turn results in simple feedback control laws,
which force the system to flow along the vector field.
We finally consider the extension of the methodology to the distributed coordination and
control for multiple nonholonomic agents. Based on the results for the single-agent case in static
obstacle environments, we propose a coordination protocol for multiple agents which need to
converge to specific goal configurations, using local information only. The proposed protocol
yields collision-free and almost globally convergent trajectories for the multi-agent system.
B. Contributions and Organization
When it comes to the single-agent case, i.e., to a robot operating in a known, static environment
of circular obstacles, the proposed method does not suffer from the appearance of sinks (stable
nodes) away from goal point. Furthermore, compared to similar feedback methods which rely
on scalar (potential) functions, such as [11], the main difference and advantage of the proposed
approach is that:
(i) no parameter tuning is needed in order to render the desired convergence properties; the
values of the parameter λ of the vector field are known a priori.
Compared to similar methods which rely on vector fields, such as [17], the proposed method:
October 22, 2014
DRAFT
IEEE TRANSACTIONS ON
6
(ii) requires neither the computation of a cell decomposition of the free space, nor the existence
of a high-level discrete motion plan, and as thus it is free of any computational complexity
issues,
(iii) addresses the motion planning and collision avoidance for multiple agents in dynamic
environments, and is scalable as the number of agents increases.
Finally, compared to other similar vector field based methods, such as [18]–[20], the proposed
method:
(iv) guarantees the convergence of the orientation trajectories of the robots to any predefined
value.
Remark 1: While here we consider circular, not polygonal, obstacle environments, preliminary
results reveal that the method can be used for defining repulsions around polygonal obstacles as
well, see the Appendix in [35].
When it comes to the multi-agent case, i.e., to multiple agents moving towards goal configu-
rations while avoiding collisions, the proposed method:
(v) offers the flexibility to directly impose the minimum allowable clearance among agents,
something which typically is not the case with gradient-based solutions. This character-
istic might be desirable, for instance, when considering multi-robot systems in confined
environments.
(vi) being a non-gradient vector field approach, the technical developments are based on set
invariance concepts rather than Lyapunov-based methods. This in principle provides less
conservative solutions, while it might desirable in extending the method to more complicated
dynamical models.
Compared to our earlier work, the vector field construction presented here is not the same
with the one in [36]. Furthermore, the proposed construction, coordination protocol and technical
developments are not the same with the ones in [37]. Moreover, since it offers feedback solutions
with certain convergence guarantees, it can be used as a basis in constrained model predictive
control designs [38], which are appropriate for uncertain environments. The case of mixed
environments, i.e., of multiple agents operating among physical obstacles under uncertainty,
are not considered in this paper and this topic is left open for future research.
October 22, 2014
DRAFT
IEEE TRANSACTIONS ON
7
Part of this work has appeared in [39]. The current paper additionally includes: (i) a detailed
presentation of the overall method both for the static and the dynamic case, along with the proofs
which have been omitted in the conference version in the interest of space, (ii) more simulation
results which demonstrate the efficacy of the method in static and dynamic environments.
The paper is organized as follows: Section II includes a brief overview of the notions regarding
the topology of two-dimensional vector fields that are used throughout the paper. Section III
characterizes the singular points of our vector fields w.r.t. the parameter λ, while section IV
presents the blending mechanism among vector fields, the construction of the almost global
feedback motion plans and the underlying control design, along with simulation results in
static obstacle environments. Section V presents the extension of the method to the distributed
coordination and collision-free motion of multiple agents under various sensing/communication
patterns. Our conclusions and thoughts on future work are summarized in Section VI.
II. SINGULAR POINTS OF VECTOR FIELDS
This section provides an overview of notions from vector field topology. For more information
the reader is referred to [34], [40], [41].
Definition 1: A vector field on an open subset U ⊂Rn is a function which assigns to each
point p ∈U a vector Xp ∈Tp(Rn). A vector field on Rn is C∞(smooth) if its components
relative to the canonical basis are C∞functions on U.
Definition 2: Given a C∞vector field X on Rn, a curve t →F(t) defined on an open interval
J of R is an integral curve of X if dF
dt = XF(t) on J.
Definition 3: A point p of U at which Xp = 0 is called a singular, or critical, point of the
vector field.
Center-type and non-center type singularities: Singular points are typically distinguished
to those that are reached by no integral curve (called center type) and those that are reached
by at least two integral curves (called non-center type). In the case of a center type singularity,
one can find a neighborhood of the singular point where all integral curves are closed, inside
one another, and contain the singular point into their interior. In the case of non-center type
singularities, one has that at least two integral curves converge to the singular point. The local
structure of a non-center type singularity is analyzed by considering the behavior of all the
integral curves which pass through the neighborhood of the singular point. This neighborhood
October 22, 2014
DRAFT
IEEE TRANSACTIONS ON
8
is made of several curvilinear sectors. A curvilinear sector is defined as the region bounded by
a circle C of arbitrary small radius, and two integral curves, S and S′, which both converge
(for either t →+∞, or t →−∞) to the singular point. The integral curves passing through the
open sector g (i.e., the integral curves except for S, S′) determine the following three possible
types of curvilinear sectors [42]: (i) Elliptic sectors: all integral curves begin and end at the
critical point. (ii) Parabolic sectors: just one end of each integral curve is at the critical point.
(iii) Hyperbolic sectors: the integral curves do not reach the critical point at all. The integral
curves that separate each sector from the next are called separatrixes, see also Fig. 1.
Fig. 1.
A typical isolated critical point. Image taken from [34].
First-order and high-order singularities: A singular point p of a vector field X on R2 is
called a first-order singular point if the Jacobian matrix JX(·) of the vector field X does not
vanish (i.e., is nonsingular) on p, i.e., if: det (JX(p)) ̸= 0; otherwise the singular point is called
high-order singular point.
October 22, 2014
DRAFT
IEEE TRANSACTIONS ON
9
III. NAVIGATION VIA VECTOR FIELDS
Consider the motion of a robot with unicycle kinematics in an environment W with N static
obstacles. The equations of motion read:


˙x
˙y
˙θ

=


cos θ
0
sin θ
0
0
1



u
ω

,
(2)
where q =

rT θ
T is the configuration vector, r = [x y]T is the position and θ is the orientation
of the robot w.r.t. a global frame G, and u, ω are the linear and the angular velocity of the robot,
respectively. The robot is modeled as a closed circular disk of radius ϱ, and each obstacle Oi
is modeled as a closed circular disk of radius ϱoi centered at roi = [xoi yoi]T, i ∈{1, . . . , N}.
Denote Oi = {r ∈R2 | ∥r −roi∥≤ϱoi}.
A. A family of vector fields for robot navigation
We consider the class of vector fields F : R2 →R2 given by (1). The vector field components
Fx, Fy read:
Fx = (λ −1)pxx2 + λpyxy −pxy2,
(3a)
Fy = (λ −1)pyy2 + λpxxy −pyx2.
(3b)
Theorem 1: The origin r = 0 is the unique singular point of the vector field F (1) if and
only if λ ̸= 1.
Proof: It is straightforward to verify that r = 0 is a singular point of F. Let us write the
vector field components (3) of F in matrix form as:

Fx
Fy

=

(λ −1)x2 −y2
λxy
λxy
(λ −1)y2 −x2


|
{z
}
A(λ,r)

px
py

.
(4)
The determinant of the matrix A(λ, r) is: det(A(λ, r)) = −(λ −1)(x2 + y2)2. This implies that
A(λ, r) is nonsingular away from the origin r = 0 if and only if λ ̸= 1. Therefore, for λ ̸= 1
and r ̸= 0, one has F = 0 if and only if p = 0. Since p ̸= 0 by definition, if follows that the
vector field F is nonsingular everywhere but the origin r = 0, as long as λ ̸= 1.
October 22, 2014
DRAFT
IEEE TRANSACTIONS ON
10
Theorem 2: The line l : y = tan ϕ x, where tan ϕ ≜py
px, is an axis of reflection, or mirror
line, for F (1).
Proof: Consider two points A, B of equal distance and on opposites sides w.r.t. the line l
(Fig. 2). Their position vectors rA = [xA yA]T, rB = [xB yB]T w.r.t. G read:
xA = R cos a,
yA = R sin a,
(5a)
xB = R cos(2ϕ −a),
yB = R sin(2ϕ −a),
(5b)
where (R, a), (R, (2ϕ −a)) are the polar coordinates of A, B, respectively. We need to prove
α
β
β
φ
y = tanφ x
xA
yA
xB
yB
υο
υ’ο
υp
υp
A
B
xg
yg
xl
yl
Fig. 2.
The line l : y = tan ϕ x, where ϕ = arctan( py
px ), is a reflection (or mirror) line for the vector field F.
that the vector F(rA), denoted FA, reflects to the vector F(rB), denoted FB, w.r.t. the line
l : y = tan ϕ x. Recall that the reflection matrix about the considered line l is:
H(2ϕ) =

cos 2ϕ
sin 2ϕ
sin 2ϕ
−cos 2ϕ

.
(6)
October 22, 2014
DRAFT
IEEE TRANSACTIONS ON
11
Substituting (5a) into (4) and after some standard algebra yields:
FA = (λ −2)R2∥p∥
2

cos ϕ
sin ϕ


|
{z
}
vp
+ λR2∥p∥
2

cos(ϕ −2a)
−sin(ϕ −2a)


|
{z
}
vo
,
(7)
where ∥p∥=
p
px2 + py2. Similarly, substituting (5b) into (4) yields:
FB = (λ −2)R2∥p∥
2

cos ϕ
sin ϕ


|
{z
}
vp
+ λR2∥p∥
2

cos 2ϕ
sin 2ϕ
sin 2ϕ
−cos 2ϕ



cos(ϕ −2a)
−sin(ϕ −2a)


|
{z
}
v′o
.
(8)
One has: FA = vp + vo and FB = vp + v′
o. Out of (7), (8) one gets that v′
o = H(2ϕ)vo, i.e., v′
o
is the reflection of the vector vo about the line l. Thus, one may write
vo = vl
ox ˆxl + vl
oy ˆyl and v′
o = vl
ox ˆxl −vl
oy ˆyl,
where ˆxl, ˆyl are the unit vectors along the axes xl, yl, respectively, see Fig. 2. Furthermore, vp
is parallel to the vector p, i.e., parallel to the candidate reflection line l. Consequently, one may
write: vp = vl
px ˆxl + 0 ˆyl. It follows that:
FA = (vl
ox + vl
px)ˆxl + vl
oy ˆyl,
FB = (vl
ox + vl
px)ˆxl −vl
oy ˆyl,
i.e., that the vector FB is a reflection of vector FA about the line l. This completes the proof.
Remark 2: The Jacobian matrix of F is singular at r = 0, which implies that r = 0 is a
high-order singularity. Thus, one may expect that the pattern of the integral curves around the
singular point will be more complicated compared to those around a first-order singularity, i.e.,
around nodes, saddles, foci or centers.
Theorem 3: The equation of the integral curves of F for p = [1 0]T is given as:
(x2 + y2)
λ
2 = c y(λ−1), c ∈R .
(9)
October 22, 2014
DRAFT
IEEE TRANSACTIONS ON
12
Proof: Consider the polar coordinates (r cos φ, r sin φ) of a point (x, y) where:
r =
p
x2 + y2,
cos φ = x
r ,
sin φ = y
r.
(10)
After substituting (10) and px = 1, py = 0 into (4) the vector field components read:
Fx = r2 (λ −1) cos2 φ −sin2 φ


,
(11a)
Fy = r2 (λ cos φ sin φ) .
(11b)
An integral curve of (1) is by definition the solution of the system of ordinary differential
equations:
dx
dt = Fx
dy
dt = Fy
,
which further reads:
dx
dy = Fx
Fy
,
(12)
while the differentials between Cartesian and polar coordinates satisfy the formula:

dr
rdφ

=

cos φ
sin φ
−sin φ
cos φ



dx
dy

.
(13)
Plugging (13), (11) into (12) results in:
1
r dr = (λ −1)cos φ
sin φ dφ,
while integrating by parts yields:
ln(r) = (λ −1) ln(sin φ) + ln(c) ⇒
ln(r) = ln
c sin(λ−1) φ


⇒
r = c sin(λ−1) φ ⇒r = c y(λ−1)
r(λ−1) ⇒
rλ = c y(λ−1) ⇒(x2 + y2)
λ
2 = c y(λ−1), where c ∈R .
This completes the proof.
Remark 3: It is straightforward to verify that:
• For λ = 0, (9) reduces to y = c, i.e., the integral curves are straight lines parallel to
p = [1
0]T.
• For λ = 1, (9) reduces to
p
x2 + y2 = c, i.e., the integral curves are circles of radius √c,
where c > 0, centered at the origin (x, y) = (0, 0).
October 22, 2014
DRAFT
IEEE TRANSACTIONS ON
13
B. Attractive vector fields
Let us consider the case λ = 2. Take for simplicity p = [1 0]T and write the vector field
components as:
Fx = x2 −y2,
(14a)
Fy = 2xy
(14b)
Following [34], the singular point r = 0 of (14) is a dipole. More specifically, the vector field
(14) has two elliptic sectors, with the axis y = 0 serving as the separatrix. This implies that
all integral curves begin and end at the singular point, except for the separatrix y = 0. The
separatrix converges to r = 0 for x < 0 and diverges for x > 0 (Fig. 3). Out of Theorem 2, the
separatrix y = 0 is the reflection line for the vector field (14).
-1.0
-0.5
0.0
0.5
1.0
-1.0
-0.5
0.0
0.5
1.0
Fig. 3.
The integral curves of (1) for λ = 2, px = 1, py = 0.
Furthermore, Theorem 2 implies that the axis the vector p ̸= 0 lies on is, in general, a
October 22, 2014
DRAFT
IEEE TRANSACTIONS ON
14
reflection line for (1). This means that the resulting integral curves are symmetric w.r.t. the
vector p ∈R2. In that sense, any of the integral curves of F offers a path to r = 0, while at the
same time the direction of the vector p dictates the symmetry axis of the integral curves w.r.t.
the global frame G.
Therefore, defining a feedback motion plan for steering the unicycle to a goal configuration
qg =

rT
g
θg
T has been based in earlier work of ours’ [33] on the following simple idea: Pick
a vector field F out of (1) in terms of (r −rg),2 with λ = 2 and p = [px py]T, so that the
direction of the vector p coincides with the goal orientation: ϕ ≜arctan( py
px) = θg. Then, the
integral curves serve as a reference to steer the position trajectories r(t) to the goal position rg,
and the orientation trajectories θ(t) to the goal orientation θg.
C. Repulsive vector fields
Let us consider the case λ = 1, i.e., the case when the vector field (1) has multiple singular
points. The vector field components read:
Fx = pyxy −pxy2,
(15a)
Fy = pxxy −pyx2.
(15b)
The vector field (15) vanishes on the set V = {r ∈R2 | pyx−pxy = 0}. Out of Theorem 2, the
singularity set V coincides with the reflection line of the vector field (15). The equation of the
integral curves can be computed for pyx −pxy ̸= 0 as: dx
dy =
y
−x ⇒x2 + y2 = c2, where c ∈R,
which implies that the integral curves are circles centered at the origin r = 0, see Fig. 4.
The signum of x (in general, of piTr) dictates whether the integral curves escape the singularity
set V (see the half-plane x > 0) or converge to the singularity set V (see the half-plane x < 0).
We say that the singular point r = 0 of the vector field (15) is of center type; this means that
no integral curve reaches the singular point.3
Thus, one may employ (15) to define tangential vector fields locally around circular obstacles.
2This is to have the unique singular point of F coinciding with the desired position rg.
3Characterizing this particular singularity as of center type is slightly inconsistent with standard notation, since in this case
the singular point r = 0 is not isolated.
October 22, 2014
DRAFT
IEEE TRANSACTIONS ON
15
-1.0
-0.5
0.0
0.5
1.0
-1.0
-0.5
0.0
0.5
1.0
Fig. 4.
The vector field F for λ = 1 and px = 1, py = 0.
IV. ALMOST GLOBAL FEEDBACK MOTION PLANS
Given the class of attractive and repulsive vector fields, the idea on defining an almost global
feedback motion plan F⋆on the collision-free space F is now simple: we pursue to combine an
attractive-to-the-goal vector field Fg with (local) repulsive vector fields Foi around each obstacle
Oi, so that the integral curves of F⋆: 1) converge to the goal qg, and 2) point into the interior
of F on the boundaries of the obstacles Oi. The vector field F⋆can then serve as a feedback
motion plan on W.
Remark 4: Combining the vector fields Fg, Foi should be done carefully so that the resulting
vector field F⋆does not have any undesired singularities on F. For this reason, we consider the
October 22, 2014
DRAFT
IEEE TRANSACTIONS ON
16
normalized unit vector fields:
Fn
g =



Fg
∥Fg∥,
for r ̸= 0;
0,
for r = 0.
(16a)
Fn
oi =



Foi
∥Foi∥,
for r /∈Vi;
0,
for r ∈Vi.
(16b)
respectively, when defining the blending mechanism, see later on in Section IV-C.
A. Attractive vector field to the goal
Without loss of generality we assume that qg = 0. An attractive-to-the-goal vector field Fg
may be taken out of (1) for λ = 2, pg = [1 0]T, which yields the vector field (14). The
components of the normalized vector field Fn
g taken out of (16a) for x ̸= 0, y ̸= 0 read:
Fn
gx = x2 −y2
x2 + y2,
Fn
gy =
2xy
x2 + y2.
B. Repulsive vector field w.r.t. a circular obstacle
Consider an obstacle Oi and the region Zi :

r ∈R2 | ∥r −roi∥≤ϱZi
	
, where ϱZi = ϱoi +
ϱ + ϱε, see Fig. 5. The parameter ϱε ≥0 is the minimum distance that the robot is allowed to
keep w.r.t. the boundary of the obstacle.
A repulsive vector field w.r.t. the point roi can be picked out of (15) for pi = [pxi pyi]T,
where pxi = cos φi, pyi = sin φi, φi = atan2(−yoi, −xoi) + π as:
Foxi = pyi(x −xoi)(y −yoi) −pxi(y −yoi)2,
(17a)
Foyi = pxi(x −xoi)(y −yoi) −pyi(x −xoi)2.
(17b)
Note that the vector pi is picked such that it lies on the line connecting the center roi of the
obstacle with the goal point rg = 0. Therefore, the singularity set Vi of (17) lies by construction
on this line, which is also the reflection axis of the vector field (17). Denote Ai = {r ∈
Zi | piT(r −roi) ≥0}, Bi = {r ∈Zi | piT(r −roi) < 0}, where Zi = Ai
S Bi, and consider
the behavior of the integral curves around the singularity set Vi. The integral curves depart from
the singularity set Vi in the region Ai (see the red vectors around Vi in Fig. 5), and converge to
the singularity set Vi in the region Bi (the corresponding vectors have not been drawn in Fig.
5). The integral curves in region Ai render safe, tangential reference paths around the obstacle
October 22, 2014
DRAFT
IEEE TRANSACTIONS ON
17
o i
ρ
i
p
ρ
i
O
ε
ρ
g
x
B i
A i
λ = 1
λ = 0
Ζ i
λ = 2
V i
gy
Fig. 5.
Defining a repulsive vector field Foi around the obstacle Oi. Note that we take the vector field (1) with λ = 1 in
region Ai and with λ = 0 in region Bi.
Oi. However, their pattern in region Bi is undesirable, since it may trap the system trajectories
r(t) away from rg. To overcome this, in region Bi we define a vector field out of (1) for λ = 0
and pi as before, whose vector field components read:
Foxi = −pxi(x −xoi)2 −pxi(y −yoi)2,
(18a)
Foyi = −pyi(x −xoi)2 −pyi(y −yoi)2.
(18b)
This vector field is co-linear with pi and vanishes at the unique singular point r = roi.
Remark 5: The transition of the integral curves between regions Ai, Bi is smooth, since the
vectors at the points where piT(r −roi) = 0 coincide.
October 22, 2014
DRAFT
IEEE TRANSACTIONS ON
18
In summary, the vector field Foi around a circular obstacle Oi is picked out of the family of
vector fields (1) as:
Foi =



F(λ=1) (δri) ,
for piT(δri) ≥0;
F(λ=0) (δri) ,
for piT(δri) < 0,
(19)
where δri ≜r −roi, φi ≜atan2(−yoi, −xoi) + π, pi = [cos φi sin φi]T. The normalized vector
field then reads:
Fn
oi =









F(λ=0)(δri)
∥F(λ=0)(δri)∥,
for piT(δri) < 0;
F(λ=1)(δri)
∥F(λ=1)(δri)∥,
for piT(δri) ≥0 and r /∈Vi;
0,
for piT(δri) ≥0 and r ∈Vi;
(20)
C. Blending attractive and repulsive vector fields
Define the obstacle function βi(·) : R2 →R as:
βi(r, roi, ϱoi) = ϱoi
2 −∥r −roi∥2,
(21)
which is positive in the interior Int(Oi) of the obstacle, zero on the boundary ∂Oi of the obstacle,
and negative everywhere else. Denote the value of the constraint function βi on the boundary
∂Zi of the region Zi as βiZ = −2ϱoi (ϱ + ϱε) −(ϱ + ϱε)2.
The repulsive vector field Fn
oi is then locally defined on the set: (Zi \ Int(Oi)) = {r ∈
R2 | βiZ ≤βi(r) ≤0}. At the same time, the attractive vector field Fn
g should be defined exterior
to Zi, i.e., for βi(r) < βiZ. To encode this, define the smooth bump function σi(·) : R2 →[0, 1]:
σi =









1,
for βi(r) ≤βiF;
aβi
3 + bβi
2 + cβi + d,
for βiF < βi(r) < βiZ;
0,
for βiZ ≤βi(r);
(22)
where βiZ is the value of (21) at distance ϱZi w.r.t. roi, βiF is the value of (21) at some distance
ϱFi > ϱZi w.r.t. roi, and the coefficients a, b, c and d are computed as:
a =
2
(βiZ −βiF)3,
b = −3(βiZ + βiF)
(βiZ −βiF)3 ,
c =
6βiZβiF
(βiZ −βiF)3,
d = βiZ
2(βiZ −3βiF)
(βiZ −βiF)3
,
October 22, 2014
DRAFT
IEEE TRANSACTIONS ON
19
so that (22) is a C2 function. Having this at hand, and inspired by [17], one may now define the
vector field:
Fi = σiFn
g + (1 −σi)Fn
oi.
(23)
Lemma 1: The vector field (23) is:
(i) Attractive to the goal qg for ∥r −roi∥≥ϱFi, i.e., for βi(r) ≤βFi where σi = 1, via the
effect of Fn
g.
(ii) Repulsive w.r.t. Oi for ϱoi ≤∥r −roi∥≤ϱZi, i.e., for βZi ≤βi(r) where σi = 0, via the
effect of Fn
oi.
(iii) Nonsingular in the region ϱZi < ∥r −roi∥< ϱFi, i.e., for βFi < βi(r) < βZi where
0 < σi < 1.
(iv) Safe w.r.t. the obstacle Oi and convergent to the goal qg for almost all initial conditions.
Proof: The first two arguments have been proved in the previous section. To verify the third
argument, consider the norm of vector field Fi in the blending region Di : {r ∈R2 | ϱZi <
∥r −roi∥< ϱFi}, which reads:
∥Fi∥=
p
1 −2σi(1 −σi) + 2σi(1 −σi) cos α,
where α the angle between the vectors Fn
g, Fn
oi at some point r ∈Di. Then, for r /∈Vi one has
that ∥Fi∥vanishes at the points where σi is the solution of:
2(1 −cos α)σi
2 −2(1 −cos α)σi + 1 = 0.
The discriminant reads ∆= −4(1 −cos α)2, which implies that there are no real solutions, i.e.,
that the vector field Fi is nonsingular for r /∈Vi. Moreover, for r ∈Vi one has Fn
oi = 0, and
therefore: ∥Fi∥= σi ̸= 0.
Finally, to verify the fourth argument, consider first that the integral curves which do not intersect
with the blending region Di are convergent by construction to rg. Consider now the boundary
Si : {r ∈R2 | ∥r −roi∥2 −ϱFi
2 = 0}
of the region Di and let us analyze the behavior of the integral curves on the manifolds:
S−
i : {r ∈R2 | ∥r −roi∥= ϱFi + δϱ},
S+
i : {r ∈R2 | ∥r −roi∥= ϱFi −δϱ},
October 22, 2014
DRAFT
IEEE TRANSACTIONS ON
20
with δϱ > 0 arbitrarily small. After some calculations:
∇Si Fi = 2(r −roi)TFn
g,
∇S−
i Fi = 2(r −roi)TFn
g.
For ∇S+
i Fi, consider the following cases:
Case 1. The vector field Fn
oi satisfies: (r −roi)TFn
oi = 0, and therefore:
∇S+
i Fi = 2σi(r −roi)TFn
g.
Then: (∇S−
i Fi)(∇S+
i Fi) > 0, which implies that the integral curves cross the switching
surface Si and enter Ai. Consider now the behavior of the integral curves in Ai. Assume
that ∇S+
i Fi = 2σi(r −roi)TFn
g > 0; this would imply that ∇Si Fi > 0 as well, i.e., that
the integral curves did not cross Si, a contradiction. Then:
∇S+
i Fi = 2σi(r −roi)TFn
g < 0,
which yields that the integral curves approach the boundary
Ti : {r ∈R2 | ∥r −roi∥2 −ϱZi
2 = 0}
of the blending region Di. Denote
T −
i : {r ∈R2 | ∥r −roi∥= ϱZi + δϱ}
and note that: ∇T −
i Fi = ∇S+
i Fi < 0, and that ∇Ti Fi = 0, since on Ti one has σi = 0.
Then, Fi ̸= 0 is tangent to Ti, which means that the integral curves slide along Ti, until
reaching region Bi.
Remark 6: The integral curves are not defined on the (unique) point on Ti where Fi = 0.
This further implies that system trajectories which either start or reach this point get stuck
away from the goal configuration.
Let us now consider the pattern of the integral curves in the vicinity of the singularity and
characterize the set of initial conditions from which the system trajectories end there. It
was shown in the previous section that the integral curves around the singularity set Vi are
departing the set, except for one integral curve which converges to Vi. For this condition to
occur the goal orientation θg should be co-linear with the line the singularity set Vi lies on.
To see why, recall that the vector field in the blending region reads: Fi = σiFn
g, and that
October 22, 2014
DRAFT
IEEE TRANSACTIONS ON
21
the vector field Fn
g should point to the singularity set Vi. Consequently, this condition arises
if and only if the obstacle is positioned such that the direction of the vector pi coincides
with the direction of the vector pg. Therefore, the set of initial conditions from which the
integral curves of Fi converge to the singularity set Vi is of Lebesgue measure zero. Note
also that if the direction of pg does not coincide with the direction of pi, then the singular
points of Fi are confined in Zi on a line segment of length ϱε, correspond to the initial
conditions from which solutions are not defined, and are reached by no integral curve.
Case 2. In region Bi one may follow a similar analysis to conclude that the integral curves exit
Bi.
In summary, the vector field (23) is safe and globally convergent almost everywhere, i.e., except
for a set of initial conditions of measure zero.
D. Motion plan in static obstacle environments
Theorem 4: Assume a workspace W of N circular obstacles Oi, i ∈{1, . . . , N}, positioned
such that the inter-obstacle distances dij = ∥roi −roj∥satisfy:
dij ≥ϱZi + ϱZj, ∀(i, j), j ∈{1, . . . , N}, j ̸= i.
(24)
Then, the vector field F⋆: R2 →R2, given as:
F⋆=
N
Y
i=1
σiFg +
N
X
i=1
(1 −σi)Foi,
(25)
where Fg is the normalized attractive vector field (16a), Foi is the normalized repulsive vector
field (20) around an obstacle Oi, and σi is the bump function (22) defined in terms of the obstacle
function βi given by (21), is a safe, almost global feedback motion plan in F, except for a set
of initial conditions of measure zero.
Proof: By construction, the first term in (25) cancels the effect of the attractive vector
field Fg where at least one of the bump functions σi = 0, i.e., in the corresponding region Zi
around obstacle Oi. At the same time the second term shapes the corresponding vector field
Foi in Zi. Thus, the attractive vector field Fg is activated through (25) only when βi < βZi
∀i ∈{1, . . . , N}, i.e., outside the regions Zi. Furthermore, setting the inter-obstacle distance
dij ≥ϱZi + ϱZj implies that the repulsive flows around obstacles do not overlap, and therefore
October 22, 2014
DRAFT
IEEE TRANSACTIONS ON
22
are both safe and almost globally convergent to the goal, as proved in Lemma 1. This completes
the proof.
Remark 7: The condition (24) reads that the minimum distance among the boundaries of the
obstacles should be at least 2(ϱ+ϱε). This clearance is not conservative or restrictive in practice,
since the parameter ϱε can be chosen arbitrarily close to zero, or even equal to zero, in case the
robot is allowed to touch the obstacle.
E. Control design and simulation results
Having (25) at hand, the control design for the unicycle (2) is now straightforward. We use
the control law:
u = ku tanh
∥r −rg∥2

,
(26a)
ω = −kω(θ −ϕ) + ˙ϕ,
(26b)
where ϕ ≜arctan(
F⋆
y
F⋆
x) is the orientation of the vector field F⋆at a point (x, y), with its time
derivative reading:
˙ϕ
(2)=

∂F⋆
y
∂x cθ +
∂F⋆
y
∂y sθ

F⋆
x −

∂F⋆
x
∂x cθ + ∂F⋆
x
∂y sθ

F⋆
y

u,
with the linear velocity u given by (26a), see in [35], and kω > 0, ku > 0. Then, the orientation
θ of the unicycle is Globally Exponentially Stable (GES) to the safe orientation ϕ, and the robot
flows along the integral curves of F⋆until converging to rg.
To demonstrate the efficacy of the proposed navigation and control design we consider the
motion of a robot in an environment with N = 10 static obstacles (Fig. 6), where the goal
position is rg = [−0.1 0.08]T. The radii of the obstacles are set equal to ϱoi = 0.03. The
blending zone Di around each obstacle Oi is illustrated between the boundary surfaces Si (black
line) and Ti (red line), respectively. The resulting collision-free path under the control law (26),
with the control gains picked equal to ku = 0.075, kr = 2.5, are depicted in blue color.
Remark 8: The integral curves of Foi in the region Ai around an obstacle Oi forces the robot
to perform a sharp maneuver in order to follow the tangential direction and avoid collision. This
in practice is plausible for unicycle-type vehicles (e.g. differentially driven mobile robots), yet
it may not be desirable for input-constrained vehicles, such as car-like vehicles and aircraft. Our
current work focuses in encoding curvature constraints via (1).
October 22, 2014
DRAFT
IEEE TRANSACTIONS ON
23
−0.6
−0.5
−0.4
−0.3
−0.2
−0.1
0
−0.4
−0.2
0
0.2
1
2
3
4
5
6
7
8
9
10
Fig. 6.
The path of a unicycle in a obstacle environment.
V. EXTENSION TO DYNAMIC ENVIRONMENTS
Consider N agents i ∈{1, 2, . . . , N} of unicycle kinematics which are assigned with the task
to converge to goal configurations qgi while avoiding collisions.
Each agent i has a circular communication/sensing region Ci of radius Rc centered at ri =
h
xi
yi
iT
, denoted as:
Ci : {r ∈R2 | ∥ri −r∥≤Rc},
and can reliably exchange information with any agent j ̸= i which lies within its communication
region Ci. In other words, we say that a pair of agents (i, j) is connected, or equivalently, that
agent j is neighbor to agent i, as long as the inter-agent distance dij = ∥ri −rj∥≤Rc. Denote
the set of neighbors j ̸= i of agent i with Ni.
October 22, 2014
DRAFT
IEEE TRANSACTIONS ON
24
Agents j ̸= i serve as dynamic (moving) obstacles to agent i. Navigating safely to an assigned
goal qgi is then reduced into finding a feedback motion plan F⋆
i such that its integral curves:
(i) point into the interior of the collision-free space Fi on the boundaries of the agents j ̸= i,
and (ii) converge to the goal qgi. Towards this end, we would like to employ a vector field F⋆
i
for each agent i as:
F⋆
i =
Y
j∈Ni
σijFgi +
X
j∈Ni
(1 −σij)Fi
oj,
(27)
where the attractive term Fgi is taken out of (16a), the bump function σij is defined later on,
and the repulsive term Fi
oj around each each agent j ̸= i is replaced with a normalized repelling
node,4 given out of:
Fi
xoj =
xi −xj
p
(xi −xj)2 + (yi −yj)2,
(28a)
Fi
yoj =
yi −yj
p
(xi −xj)2 + (yi −yj)2.
(28b)
In order to utilize the (almost global) convergence and safety guarantees applying to the static
case, we need to ensure that the repulsive flows around agents do not overlap, for any pair of
agents (i, j). Recall from the static case that this condition equivalently reads as that the minimum
distance dm between any pair of (moving, in the dynamic case) agents (i, j) is dm = 2(2ϱ+ϱϵ),
or equivalently, that the minimum clearance between any pair of agents is (2ϱ + ϱϵ), where
ϱϵ > 0 arbitrarily small and ϱ is the radius of the agents.
In that respect, the bump function σij in (27) is defined as:
σij =









1,
for dm ≤dij ≤dr;
a dij
3 + b dij
2 + c dij + d,
for dr < dij < dc;
0,
for dij ≥dc;
(29)
where the coefficients a, b, c, d have been computed as:
a = −
2
(dr −dc)3,
b = 3(dr + dc)
(dr −dc)3 ,
c = −
6 drdc
(dr −dc)3,
d = dc
2(3dc −dr)
(dr −dc)3
,
4The tangential repulsive vector field (19) defined for static obstacles is not a suitable choice for the dynamic case; the reason
is that the repulsive integral curves of the vector field (25) of agent i around agent j are rendered an invariant set under the
proposed velocity coordination protocol, forcing thus the trajectories ri(t), rj(t) of a pair of agents i, j converge to undesired
locations away from the goal locations rgi, rgj, see also the analysis in [?].
October 22, 2014
DRAFT
IEEE TRANSACTIONS ON
25
so that (22) is a C2 function, dm ≥2(2ϱ + ϱϵ), dm < dr < dc.
Remark 9: The communication/sensing range of each agent i should be Rc ≥dc.
A. Control design
Each agent i moves under the control law:
ui =





max

0,
min
j∈Ni|Jj<0 ui|j

,
dm ≤dij ≤Rc,
uic,
Rc < dij;
,
(30a)
ωi = −kωi (θi −ϕi) + ˙ϕi,
(30b)
where: ϕi is the orientation of the vector field F⋆
i at a point (x, y), the vector field F⋆
i is given
by (27), ui|j is the safe velocity of agent i w.r.t. an agent j lying in the communication region
of Ci of agent i, given as:
ui|j = uic
dij −dm
Rc −dm
+ εi uis|j
Rc −dij
Rc −dm
,
(31)
with the terms in (31) defined as:
uic = kui tanh(∥ri −rgi∥),
uis|j = uj
rjiTηj
rjiTηi
,
ηi =

cos ϕi
sin ϕi

,
Jj = rji
Tηi,
rji = ri −rj,
and εi > 1, kui, kωi > 0.
Theorem 5: Consider N agents i ∈{1, . . . , N} assigned to converge to goal configurations
qgi. Then, under the control law (30) each agent safely converges to its goal configuration almost
globally, except for a set of initial conditions of measure zero.
Proof: The closed loop trajectories of each agent i are forced to flow along the vector field
(27). If dij(t) > Rc, ∀t ≥0 and ∀j ∈{1, . . . , N}, then σij(t) = 1, implying that agent i flows
safely along (16a) and converges to qgi.
Let us now assume that at some time t ≥0 the distance dij(t) between a pair of agents (i, j) is
dij(t) ≤Rc. By definition agent i lies in the sensing/communication region of agent j and vice
versa, which implies that they exchange information on their current positions ri(t), rj(t) and
October 22, 2014
DRAFT
IEEE TRANSACTIONS ON
26
velocities νi(t), νj(t). Consider the blending region Di : {rj ∈R2 | dr < ∥ri −rj∥< dc} and
the surfaces:
Si(t) : {ri(t), rj(t) ∈R2 | ∥ri(t) −rj(t)∥−dc = 0},
Ti(t) : {ri(t), rj(t) ∈R2 | ∥ri(t) −rj(t)∥−dr = 0}.
Lemma 2: Agent i avoids collision with any of its neighbor agents j ∈Ni.
Proof: Collision-free motion is realized as ensuring that dij(t) ≥2ϱ, ∀t ≥0, for any pair
(i, j). Let us consider the time derivative of inter-agent distance function, which after some
calculations reads:
d
dt dij = (xi −xj)( ˙xi −˙xj)
dij
+ (yi −yj)( ˙yi −˙yj)
dij
(2)= ui rjiTηi −uj rjiTηj
dij
.
(32)
The control law (30) renders the value of the time derivative (32) positive when the value of the
distance function is dij = dm, implying thus that the inter-agent distance is forced to increase.
Since dm > 2ϱ, this further implies that collisions are avoided.
In order to draw conclusions about the convergence of the agents’ trajectories to their goal
configurations we need to examine the behavior of the integral curves around the switching
surfaces Si(t), Ti(t). With the vector fields F⋆
i , F⋆
j well-defined everywhere in the corresponding
blending regions, and the linear velocities ui, uj vanishing only at the goal locations rgi, rgj, we
are interested in identifying conditions under which the system trajectories ri(t), rj(t) are forced
to get stuck on Si(t), or on Ti(t), for infinite amount of time. This can be seen as identifying
sufficient conditions of the appearance of (chattering) Zeno behavior, or Zeno points [43]. A
sufficient condition on the appearance of Zeno points is given in [44], Theorem 2. Based on
this result, we study under which conditions the system (i.e., agents’) trajectories converge to
a Zeno point. Consider the case with N = 2 agents. Denote the dynamics of the k-th agent as
˙qk = fk(qk), k ∈{i, j}, q =

qiT qjTT, r =

riT rjTT, and take:
∇Sif(q) = 2ui rij
T

cos θi
sin θi

−2uj rij
T

cos θj
sin θj

,
(33)
where rij = ri −rj. Note that the control law (30b) renders the orientation θk of the k-th agent
GES to the orientation ϕk of the vector field F⋆
k. Thus, the unit vector [cos θk
sin θk]T coincides
October 22, 2014
DRAFT
IEEE TRANSACTIONS ON
27
with the vector field F⋆
k(rk), evaluated at rk ∈R2. With this at hand and after some algebraic
calculations one has:
∇S−
i f(q) = 2ui rij
TFgi −2uj rij
TFgj
|
{z
}
A
,
∇S+
i f(q) = σij
2ui rij
TFgi −2uj rij
TFgj


+ (1 −σij)
2ui rij
TFi
oj −2uj rij
TFj
oi


|
{z
}
B
.
The set of Zeno points is:
Zi = {r ∈R2N | ∇S−
i f(q) = ∇S+
i f(q) = 0},
which reads:
A = σijA + (1 −σij)B = 0 ⇒A = B = 0 ⇒
ui rij
T(Fgi −Fi
oj) = uj rij
T(Fgj −Fj
oi).
Not surprisingly, the set Zi is depended on the current positions ri, rj and the goal locations
qgi, qgj. The Zeno condition reduces to (Fgi + Fgj) = 0, which corresponds to current positions
ri(t), rj(t) and goal locations rgi, rgj lying on the same line. Then, the set of initial conditions
(positions) from which agents’ trajectories converge to the set Zi is confined on R, i.e., on a
lower dimensional manifold, and as thus is of measure zero. The same analysis holds along
the switching surface Ti(t), yielding exactly the same condition as before regarding on the
appearance of Zeno points.
The case of N > 2 agents can be treated accordingly. Consider an agent i lying at distance
dim ≤dc w.r.t. M ≤(N −1) agents m ̸= i. The vector field F⋆
i includes the repulsive effect
Fi
om of all M connected agents. To check whether undesired singularities appear, one needs to
consider the norm ∥Fi∥in the blending region Di. The analytical expression is more involved
compared to the N = 2 case. To make the probability of more than one agents lying in the
blending region Di as low as possible, one can define the width dc −dr →0. Define also:
Sim(t) : {ri(t), rm(t) ∈R2 | ∥ri(t) −rm(t)∥−dc = 0}
the M ≤(N −1) switching surfaces of agent i w.r.t. its neighbors m. The conditions on the
appearance of Zeno points around each switching surface read:
∇S−
imf(ri, rm) = ∇S+
imf(ri, rm) = 0, ∀m ∈{1, . . . , N}.
October 22, 2014
DRAFT
IEEE TRANSACTIONS ON
28
This results in NM
2
switching surfaces, since for any pair of agents (i, m) it holds that: Sim = Smi,
and NM Zeno conditions.
Now, note that the NM
2
Zeno conditions S−
imf(ri, rm) = 0 introduce 2 NM
2
= NM unknown
terms of the form rimTFgi, rimTFgm.
In the same spirit, the NM
2
Zeno conditions S+
imf(ri, rm) = 0 additionally introduce NM
2 (M −
1)2 = NM(M −1) unknown terms of the form rimTFi
ok, rimTFm
ok.
In total, one has NM 2 unknown terms and NM equations. Given that the N goal locations
rgi, rgm, . . . , are known, the number of unknown terms reduces to NM 2−NM = NM(M −1).
To have as many equations as unknown terms, it should hold that M −1 = 1 ⇒M = 2. This
implies that each agent i ∈{1, . . . , N} is connected with at most M = 2 agents; note that this
is irrespective of the total number of agents N. Then, the geometric conditions which result in
Zeno points are given as the solutions of the resulting linear system; these solutions express NM
relations of the form rT
imFi
ok, rT
imFm
ok, which dictate the Zeno points, i.e., the Zeno positions
among the N agents. Then, the set of initial conditions from which the agents converge to these
Zeno positions are confined to a lower dimensional manifold, since they correspond to initial
positions confined on a line, and to a specific initial orientation for each agent, and as thus are
of measure zero.
Finally, let us note that the case of M > 2 neighbors is not of interest for the proposed algorithm,
as each agent i makes the avoidance decision w.r.t. the worst-case neighbor agent, i.e., w.r.t. the
agent which is more susceptible to collision. This is realized via considering the safe velocity
ui|m w.r.t. each neighbor agent m and taking the minimum over safe velocities in the definition
of the linear velocity control law (30a). The maximum function is defined to ensure that each
agent i will never be forced to move with negative linear velocity, i.e., backwards; this is to
ensure that there is no possibility of back-to-back colliding agents.
In summary, the motion of each agent i remains collision-free w.r.t. its neighbor agents j ∈Ni
under the control law (30), and each agent i converges to its goal location qgi almost globally,
except for a set of initial configurations of measure zero. This completes the proof.
Remark 10: Theorem 5 justifies that the set of initial conditions for which the multi-robot
system exhibits Zeno trajectories (chattering across a switching surface for infinite amount of
time) which result in robots getting stuck away from their goals, is of measure zero. To avoid
sliding along a switching surface, which can be seen as “finite-time chattering”, one can employ
October 22, 2014
DRAFT
IEEE TRANSACTIONS ON
29
hysteresis logics [45].
B. Simulation Results
We consider N = 30 agents which are moving towards their goal locations (depicted with
square markers) starting from goal positions (depicted with cross markers) while avoiding colli-
sions, see the resulting paths in Fig. ??. The goal locations are defined sufficiently far apart so
that the communication regions do not overlap when agents lie on their goal locations.
VI. CONCLUSIONS
This paper presented a novel methodology for the motion planning of unicycle robots in
environments with obstacles, with extensions to the collision avoidance in multi-agent systems.
The method is based on a family of vector fields whose integral curves exhibit attractive or
repulsive behavior depending on the value of a parameter. It was shown that attractive-to-the-
goal and repulsive-around-obstacles vector fields can be suitably blended in order to yield almost
global feedback motion plans in environments with circular obstacles. The case of collision
avoidance under local sensing/communication in multi-agent scenarios was also treated. No
parameter tuning is needed in order to avoid local minima, as needed in similar methods which
are based on scalar (potential) functions. Current work focuses on the definition of vector fields
encoding input constraints, such as curvature bounds, which may be more appropriate for aircraft
and car-like vehicles.
REFERENCES
[1] L. E. Parker, “Path planning and motion coordination in multiple mobile robot teams,” in Encyclopedia of Complexity and
System Science, R. A. Meyers, Ed.
Springer, 2009.
[2] H. Kress-Gazit, G. E. Fainekos, and G. J. Pappas, “Temporal-logic-based reactive mission and motion planning,” IEEE
Transactions on Robotics, vol. 25, no. 6, pp. 1370–1381, Dec. 2009.
[3] A. Bhatia, M. R. Maly, L. E. Kavraki, and M. Y. Vardi, “A multi-layered synergistic approach to motion planning with
complex goals,” IEEE Robotics and Automation Magazine, vol. 18, no. 3, pp. 55–64, 2011.
[4] W. Ren and Y. Cao, “Overview of recent research in distributed multi-agent coordination,” in Distributed Coordination of
Multi-agent Networks, ser. Communications and Control Engineering.
Springer-Verlag, 2011, ch. 2, pp. 23–41.
[5] D. Panagou, D. M. Stipanovi´c, and P. G. Voulgaris, “Vision-based dynamic coverage control for nonholonomic agents,”
in Proc. of the 53rd IEEE Conference on Decision and Control, Los Angeles, CA, Dec. 2014, p. to appear.
[6] D. Panagou and K. J. Kyriakopoulos, “Viability control for a class of underactuated systems,” Automatica, vol. 49, no. 1,
pp. 17–29, Jan. 2013.
October 22, 2014
DRAFT
IEEE TRANSACTIONS ON
30
Fig. 7.
Collision-free motion of 30 nonholonomic agents under the proposed control strategy.
October 22, 2014
DRAFT
IEEE TRANSACTIONS ON
31
Fig. 8.
Collision-free motion of 30 nonholonomic agents under the proposed control strategy (Continued).
[7] G. Leitmann and J. Skowronski, “Avoidance control,” Journal of Optimization Theory and Applications, vol. 23, pp.
581–591, Dec. 1977.
[8] O. Khatib, “Real-time obstacle avoidance for manipulators and mobile robots,” The International Journal of Robotic
Research, vol. 5, no. 1, pp. 90–98, Spring, 1986.
[9] E. G. Hernandez-Martinez and E. Aranda-Bricaire, “Convergence and collision avoidance in formation control: A survey
of the artificial potential functions approach,” in Multi-Agent Systems - Modeling, Control, Programming, Simulations and
Applications, F. Alkhateeb, E. A. Maghayreh, and I. A. Doush, Eds.
InTech, 2011, ch. 6, pp. 103–126.
[10] D. M. Stipanovi´c, C. J. Tomlin, and G. Leitmann, “Monotone approximations of minimum and maximum functions and
multi-objective problems,” Applied Mathematics and Optimization, vol. 66, pp. 455–473, 2012.
[11] E. Rimon and D. Koditschek, “Exact robot navigation using artificial potential functions,” IEEE Transactions on Robotics
and Automation, vol. 8, no. 5, pp. 501–518, Oct. 1992.
[12] C. I. Connolly, J. B. Burns, and R. Weiss, “Path planning using Laplace’s equation,” in Proc. of the 1990 IEEE Int. Conf.
on Robotics and Automation, May 1990, pp. 2102–2106.
[13] P. Szulczy´nski, D. Pazderski, and K. Kozłowski, “Real-time obstacle avoidance using harmonic potential functions,” Journal
of Automation, Mobile Robotics and Intelligent Systems, vol. 5, no. 3, pp. 59–66, 2011.
[14] H. J. S. Feder and J.-J. E. Slotine, “Real-time path planning using harmonic potentials in dynamic environments,” in Proc.
of the 1997 IEEE Int. Conf. on Robotics and Automation, Albuquerque, New Mexico, Apr. 1997, pp. 874–881.
[15] S. Waydo and R. M. Murray, “Vehicle motion planning using stream functions,” in Proc. of the 2003 IEEE Int. Conf. on
Robotics and Automation, Taipei, Taiwan, Sep., pp. 2484–2491.
[16] A. D. Luca and G. Oriolo, “Local incremental planning for nonholonomic mobile robots,” in Proc. of the 1994 IEEE Int.
Conf. on Robotics and Automation, May 1994, pp. 104–110.
[17] S. R. Lindemann and S. M. LaValle, “Simple and efficient algorithms for computing smooth, collision-free feedback laws
over given cell decompositions,” The International Journal of Robotics Research, vol. 28, no. 5, pp. 600–621, 2009.
[18] T. Liddy, T.-F. Lu, P. Lozo, and D. Harvey, “Obstacle avoidance using complex vector fields,” in Proc. of the 2008
Australasian Conference on Robotics and Automation, Canberra, Australia, Dec. 2008.
[19] E. G. Hern´andez-Mart´ınez and E. Aranda-Bricaire, “Multi-agent formation control with collision avoidance based on
discontinuous vector fields,” in Proc. of the 35th Annual Conf. of IEEE Industrial Electronics, Nov. 2009, pp. 2283 –2288.
October 22, 2014
DRAFT
IEEE TRANSACTIONS ON
32
[20] W. Dixon, T. Galluzo, G. Hu, and C. Crane, “Adaptive velocity field control of a wheeled mobile robot,” in Proc. of Fifth
Int. Workshop on Robot Motion and Control, Jun. 2005, pp. 145–150.
[21] M. Mesbahi and M. Egerstedt, Graph Theoretic Methods in Multiagent Networks.
Princeton University Press, 2010.
[22] W. Ren, R. W. Beard, and E. M. Atkins, “Information consensus in multi-vehicle cooperative control,” IEEE Control
Systems Magazine, vol. 27, no. 2, pp. 71–82, apr 2007.
[23] R. Olfati-Saber, J. A. Fax, and R. M. Murray, “Consensus and cooperation in networked multi-agent systems,” Proceedings
of the IEEE, vol. 97, no. 1, pp. 215–233, 2007.
[24] D. V. Dimarogonas and K. J. Kyriakopoulos, “Connectedness preserving distributed swarm aggregation for multiple
kinematic robots,” IEEE Transactions on Robotics, vol. 24, no. 5, pp. 1213–1223, Oct. 2008.
[25] S. G. Loizou and K. J. Kyriakopoulos, “Navigation of multiple kinematically constrained robots,” IEEE Transactions on
Robotics, vol. 24, no. 1, pp. 221–231, Feb. 2008.
[26] P. Lin, Y. Jia, and L. Li, “Distributed robust h∞consensus control in directed networks of agents with time-delay,” Systems
& Control Letters, vol. 57, pp. 643–653, 2008.
[27] J. Mei, W. Ren, and G. Ma, “Distributed containment control for lagrangian networks with parametric uncertainties under
a directed graph,” Automatica, vol. 48, pp. 653–659, 2012.
[28] Z. Qu, C. Li, and F. Lewis, “Cooperative control with distributed gain adaptation and connectivity estimation for directed
networks,” International Journal of Robust and Nonlinear Control, 2012.
[29] J. liang Zhang, D. lian Qi, and M. Yu, “A game theoretic approach for the distributed control of multi-agent systems
under directed and time-varying topology,” International Journal of Control, Automation, and Systems, vol. 12, no. 4, pp.
749–758, 2014.
[30] Z. Li, G. Wen, Z. Duan, and W. Ren, “Designing fully distributed consensus protocols for linear multi-agent
systems with directed graphs,” IEEE Transactions on Automatic Control, to appear, 2014. [Online]. Available:
http://arxiv.org/abs/1312.7377v2
[31] W. Ren, “On consensus algorithms for double-integrator dynamics,” in Proc. of the 46th Conference on Decision and
Control, New Orleans, LA, USA, Dec. 2007, pp. 2295–2300.
[32] K. Liu, G. Xie, W. Ren, and L. Wang, “Consensus for multi-agent systems with inherent nonlinear dynamics under directed
topologies,” Systems & Control Letters, vol. 62, no. 2, pp. 152–162, Feb. 2013.
[33] D. Panagou, H. G. Tanner, and K. J. Kyriakopoulos, “Control of nonholonomic systems using reference vector fields,”
in Proc. of the 50th IEEE Conf. on Decision and Control and European Control Conf., Orlando, FL, Dec. 2011, pp.
2831–2836.
[34] M. Henle, A Combinatorial Introduction to Topology.
Dover Publications, 1994.
[35] D. Panagou, “Motion planning and collision avoidance using non-gradient navigation vector fields. Technical Report.”
[Online]. Available: www.arxiv.org
[36] D. Panagou and V. Kumar, “Cooperative visibility maintenance for leader-follower formations in obstacle environments,”
IEEE Transactions on Robotics, vol. 30, no. 4, pp. 831–844, Aug. 2014.
[37] D. Panagou, D. M. Stipanovi´c, and P. G. Voulgaris, “Multi-objective control for multi-agent systems using Lyapunov-
like barrier functions,” in Proc. of the 52nd IEEE Conference on Decision and Control, Florence, Italy, Dec. 2013, pp.
1478–1483.
[38] S. Maniatopoulos, D. Panagou, and K. J. Kyriakopoulos, “A MPC scheme for the navigation of a nonholonomic vehicle with
field-of-view constraints,” in Proc. of the 2013 American Control Conf., Washington DC, USA, Jun. 2013, pp. 3967–3972.
October 22, 2014
DRAFT
IEEE TRANSACTIONS ON
33
[39] D. Panagou, “Motion planning and collision avoidance using navigation vector fields,” in Proc. of the 2014 IEEE
International Conference on Robotics and Automation, Hong Kong, China, Jun. 2014, pp. 2513–2518.
[40] W. M. Boothby, An Introduction to Differentiable Manifolds and Riemannian Geometry - 2nd Edition.
Academic Press,
1986.
[41] J. M. Lee, Introduction to Smooth Manifolds.
Springer, 2002.
[42] X. Tricoche, G. Scheuermann, and H. Hagen, “A topology simplification method for 2D vector fields,” in Proc. of
Visualization 2000, Salt Lake City, UT, USA, Oct., pp. 359–366.
[43] A. D. Ames and S. Sastry, “Characterization of Zeno behavior in hybrid systems using homological methods,” in Proc. of
the 2005 American Control Conf., Portland, OR, USA, Jun. 2005, pp. 1160–1165.
[44] F. Ceragioli, “Finite valued feedback laws and piecewise classical solutions,” Nonlinear Analysis, vol. 65, no. 5, pp.
984–998, 2006.
[45] D. Liberzon, Switching in Systems and Control.
Birkhauser Boston, 2003.
APPENDIX
Here we present some preliminary ideas on the extension of the method to polygonal envi-
ronments.
Consider the pattern of the integral curves for λ < −1, shown in Fig. 9. The repulsive nature
of the integral curves w.r.t. the axis the vector p lies on can be used to define a repulsive flow
w.r.t. each side of polygonal obstacles, as shown in Fig. 10. The effect of the repulsive flows
can be confined around the polygonal obstacle using blending mechanisms as those presented
in Section 4. Identifying sufficient minimum clearance around the obstacles which guarantees
the almost global convergence of the integral curves to a goal configuration in such a polygonal
environment is currently ongoing work, and beyond the scope and the length of the current
paper.
October 22, 2014
DRAFT
IEEE TRANSACTIONS ON
34
-3
-2
-1
0
1
2
3
-3
-2
-1
0
1
2
3
Fig. 9.
The integral curves for λ = −1, px = −1, py = 0.
−0.68
−0.66
−0.64
−0.62
−0.6
−0.58
−0.56
−0.54
−0.52
−0.18
−0.16
−0.14
−0.12
−0.1
−0.08
−0.06
−0.04
Fig. 10.
Integral curves around a polygonal obstacle.
October 22, 2014
DRAFT