1

A guiding vector field algorithm for path following
control of nonholonomic mobile robots

arXiv:1610.04391v2 [cs.SY] 7 Feb 2017

Yuri A. Kapitanyuk, Anton V. Proskurnikov and Ming Cao

Abstract—In this paper we propose an algorithm for pathfollowing control of the nonholonomic mobile robot based on the
idea of the guiding vector field (GVF). The desired path may be
an arbitrary smooth curve in its implicit form, that is, a level
set of a predefined smooth function. Using this function and the
robot’s kinematic model, we design a GVF, whose integral curves
converge to the trajectory. A nonlinear motion controller is then
proposed which steers the robot along such an integral curve,
bringing it to the desired path. We establish global convergence
conditions for our algorithm and demonstrate its applicability
and performance by experiments with real wheeled robots.
Index Terms—Path following, vector field guidance, mobile
robot, motion control, nonlinear systems

I. I NTRODUCTION
Many applications in industrial and mobile robots is built
upon the functionality of following accurately a predefined
geometric path [1], [2]. Path following is one of the central
problems in automatic guidance of mobile robots, such as
aerial vehicles [3], underwater [4] and surface [5] marine
crafts, wheeled ground robots and autonomous cars [6], [7].
Although mathematical models of robots may differ, the
principles of their path following control are similar and can
be classified in several major categories.
One approach to steer to the desired path is to fix its
time-parametrization; the path is thus treated as a function of
time or, equivalently, the trajectory of some reference point.
Path steering is then reduced to reference-point following,
which is a special case of the reference tracking problem
that has been extensively studied in control theory and solved
for a very broad class of nonlinear systems [8], [9]. An
obvious advantage of the trajectory tracking approach is its
applicability to a broad range of paths that can be e.g. nonsmooth or self-intersecting. However, this method has unavoidable fundamental limitations, especially when dealing with
general nonlinear systems. As shown in [10], unstable zero
dynamics generally make it impossible to track a reference
trajectory with arbitrary predefined accuracy; more precisely,
Yuri Kapitanyuk and Ming Cao are with the ENTEG institute at the
University of Groningen, The Netherlands. Anton V. Proskurnikov is with
the Delft Center for Systems and Control at Delft University of Technology,
Delft, The Netherlands, and also with Institute for Problems of Mechanical
Engineering of Russian Academy of Sciences (IPME RAS) and ITMO
University, St. Petersburg, Russia. E-mail: i.kapitaniuk@rug.nl,
anton.p.1982@ieee.org, m.cao@rug.nl
Y. Kapitanyuk and M. Cao acknowledge the partial support from the
European Research Council (ERCStG- 307207). The work of A. Proskurnikov
is supported by NWO TTW, the Netherlands, under the project STW#13712
“From Individual Automated Vehicles to Cooperative Traffic Management Predicting the benefits of automated driving through on-road human behaviour
assessment and traffic flow models (IAVTRM)”.

the integral tracking error is uniformly positive independent
of the controller design. Furthermore, in practice the robot’s
motion along the trajectory can be quite “irregular” due to
its oscillations around the reference point. For instance, it is
impossible to guarantee either path following at a precisely
constant speed, or even the motion with a perfectly fixed
direction. Attempting to keep close to the reference point, the
robot may “overtake” it due to unpredictable disturbances.
As has been clearly illustrated in the influential paper [11],
these performance limitations of trajectory tracking can be
removed by carefully designed path following algorithms.
Unlike the tracking approach, path following control treats
the path as a geometric curve rather than a function of
time, dealing with its implicit equations or some time-free
parametrization. The algorithm thus becomes “flexible to use
a timing law as an additional control variable” [11]; this
additional degree of freedom allows to maintain a constant
forward speed or any other desired speed profile, which is
extremely important e.g. for aerial vehicles, where the lifting
force depends on the robot’s speed. The dynamic controller,
maintaining the longitudinal speed, is usually separated from
“geometric” controller, steering the robot to the desired path.
A widely used approach to path-following, originally proposed for car-like wheeled robots [6], [12], assumes the
existence of the projection point, that is, the closest point
on a path, and the robot’s capability to measure the distance
to it (sometimes referred to as the “cross-track error”). The
robot’s mathematical model is represented in the Serret-Frenet
frame, consisting of the tangent and normal vector to the
trajectory at the projection point. This representation allows
to design efficient path following controllers for autonomous
wheeled vehicles [6], [7], [12] that eliminate the cross-track
error and maintain the desired vehicle’s speed along the path.
Further development of this approach leads to algorithms for
the control of complex unmanned vehicles such as cars with
multiple trailers [13] and agricultural tractors [14]. Projectionbased sliding mode algorithms [7], [14] are capable to cope
with uncertainties, caused by the non-trivial geometry of the
path, lateral drift of the vehicle and actuator saturations.
The necessity to measure the distance to the track imposes
a number of limitations on the path-following algorithm. Even
if the nearest point is unique, the robot should either be
equipped with special sensors [7] or solve real-time optimization problems to find the cross-track error. In general,
the projection point cannot be uniquely determined when e.g.
the robot passes a self-intersection point of the path or its
position is far from the desired trajectory. A possible way
to avoid these difficulties has been suggested in [15] and is

2

referred to as the “virtual target” approach. The Serret-Frenet
frame, assigned to the projection point in the algorithm from
[6], [12], can be considered as the body frame of a virtual
target vehicle to be tracked by the real robot. A modification
of this approach, offered in [15], allows the virtual target
to have its own dynamics, taken as one of the controller’s
design parameter. The design from [15] is based on the path
following controller from [12], using the model representation
in the Serret-Frenet frame and taking the geometry of the
path into account. However, the controller from [15] implicitly
involves target tracking since the frame has its own dynamics.
Avoiding the projection problem, the virtual target approach
thus inherits disadvantages of the usual target tracking. In
presence of uncertainties the robot may slow down and turn
back in order to trace the target’s position [3].
A guidance strategy of a human helmsman inspired another
path following algorithm, referred to as the line-of-sight (LOS)
method [3], [5], [16] which is primarily used for air and marine
crafts. Maintaining the desired speed of the robot, the LOS
algorithm steers its heading along the LOS vector which starts
at the robot’s center of gravity and ends at the target point.
This target is located ahead of the robot either on the path
[3] or on the line, tangent to the path at the projection point
[5]. Unlike the virtual target approach, the target is always
chosen at a fixed prescribed distance from the robot, referred
to as the lookahead distance. The maneuvering characteristics
substantially depend on the lookahead distance: the shorter
distance is chosen, the more “aggressive” it steers. A thorough
mathematical examination of the LOS method has been carried
out in the recent paper [5], establishing the uniform semiglobal exponential stability (USGES) property.
Differential-geometric methods for invariant sets stabilization [8], [9], [17] have given rise to a broad class of “setbased” [18] path following algorithms. Treating the path as a
geometric set, the algorithm is designed to make it invariant
and attractive (globally or locally). Typically the path is
considered as a set where some (nonlinear) output of the
system vanishes, and thus the problem of its stabilization boils
down to the output regulation problem. To solve it, various
linearization techniques have been proposed [17], [19]–[25].
For stabilization of a closed strictly convex curve an elegant
passivity-based method has been established in [18].
In this paper we develop a path following strategy, based
on the idea of reference vector field. Vector field algorithms
are widely used in collision-free navigation and extremum
seeking problems [26]–[29]; their efficiency in path following
problems has been recently demonstrated in [3], [30]–[32].
A vector field is designed in a way that its integral curves
approach the path asymptotically. Steering the robot along the
integral curves, the control algorithm drives it to the desired
path. Unlike many path following algorithms, providing convergence only in a sufficiently small vicinity of the desired
path, the vector field algorithm guarantees convergence of any
trajectory, which does not encounter the “critical” points where
the vector field is degenerate. In particular, in any invariant
domain without critical points the convergence to the path can
be proved.
For holonomic robots described by a single or double

Fig. 1. The structure of motion control systems

integrator model a general vector-field algorithm for navigation
along a general smooth curve in n-dimensional space has been
discussed in [33]. However, for more realistic nonholonomic
vehicles the vector-field algorithms have been studied mainly
for straight lines and circular paths [30], [31], [34], where they
demonstrate better, in several aspects, performance compared
to other approaches [3]. Unlike [3], [30], [31], in this paper
we propose and analyze rigorously a vector-field algorithm
for guidance of a general nonholonomic robot along a general
smooth planar path, given in its implicit form.
The paper is organized as follows. The path following problem is set up in Section II. In Section III we design the guiding
vector field and discuss the properties of its integral curves.
Section IV offers the path following control algorithm and
establishes its main properties. This algorithm is practically
validated by experiments with wheeled robots, described in
Section V. In Section VI we give a detailed comparison of our
algorithm with several alternative approaches to path following
controllers design and also discuss its possible extensions.
II. P ROBLEM S ETUP
A widely used paradigm in path following control is to
decompose the controller into an “inner” and an “outer”
feedback loop [35]–[37] as illustrated in Fig. 1. The “inner”
dynamic controller is responsible for maintaining the vector
of generalized velocities (e.g. the longitudinal speed and turn
rate) by controlling the vector of forces and moments. The
design of the dynamic controller is based on the robot’s
mathematical model and may include rejection of external
disturbances, e.g. adaptive drift compensators [38]. Assuming
the dynamic controller to be sufficiently fast and precise, one
may consider the speed and the turn rate to be new control
inputs, describing thus the robot with a simpler kinematic
model. The path following algorithm is typically implemented
in the “outer” kinematic controller, steering the simplified
kinematic model to the prescribed path.
The kinematic model of a mobile robot can be holonomic
or nonholonomic. A holonomic robot is not restricted by
its angular orientation and able to move in any direction,
as exemplified by helicopters, wheeled robots with omnidirectional wheels [2] and fully actuated marine vessels at low
speeds [16]. For nonholonomic robots only some directions of
motion are possible. The simplest model of this type is the
unicycle, which can move only in the longitudinal direction
while the lateral motion is impossible. Examples of robots,
reducing under certain conditions to the unicycle-type kinematics, include differential drive wheeled mobile robots [13],

3

[39], fixed wing aircraft [30], [31], [36], marine vessels at
cruise speeds [35] and car-like vehicles, whose rear wheels are
not steerable [13], [20], [39]. The most interesting are unicycle
models where the speed is restricted to be sufficiently high,
making it impossible to reduce the unicycle to a holonomic
model via the feedback linearization [40]. The lift force of a
fixed-wing UAV, the rudder’s yaw moment of a marine craft
and the turn rate of a car-like robot depend on the longitudinal
speed; at a low speed their manoeuvrability is very limited.
A. The robot’s model
In this paper, we consider the unicycle-type model where
the longitudinal velocity ur > 0 is a predefined constant
 


ẋ
cos α
r̄˙ =
= ur m̄(α) ∈ R2 , m̄(α) =
,
ẏ
sin α
(1)
α̇ = ω.

Here r̄ is the position of the robot’s center of gravity C in the
inertial Cartesian frame of reference 0XY , α is the robot’s
orientation in this frame and m̄(α) is the unit orientation vector
(see Fig. 2). The only control input to the system (1) is the
angular velocity ω, and hence the system is underactuated.

m̄
n̄
α

ϕ(x, y) = ϕ(x∗ , y∗ )

τ̄

ϕ(x, y) = const > 0

y∗
δ

Y

m̄d

P : ϕ(x, y) = 0
v̄

∆

e(x, y) = ψ[ϕ(x, y)] ∈ R.

(4)

By definition, e = 0 if and only if (x, y) ∈ P. Our goal is
to design a path following algorithm, i.e. a causal feedback
law (x(·), y(·), α(·)) 7→ ω(·), which eliminates the tracking
error |e(t)| −−−→ 0, bringing thus the robot to the predefined
t→∞
path P. Upon reaching the desired path, the algorithm should
“stabilize” the robot on it, which means that the robot’s
heading is steered along the tangent vector to the path.
The mapping ψ(·) in (4) is a free parameter of the algorithm. Formally one can get rid of this parameter, replacing
ϕ by the composition ψ ◦ ϕ. However, it is convenient to
distinguish between the path-defining function ϕ(x, y) and the
tracking error, depending on the choice of ψ(s). The path
representation is usually chosen as simple as possible: for
instance, dealing with a straight line, it is natural to choose
linear ϕ(x, y), while the circular path is naturally described
by ϕ(x, y) = (x − x0 )2 + (y − y0 )2 . At the same time,
some mathematical properties of the algorithm, in particular,
the region where convergence of the algorithm is guaranteed,
depend on the way the tracking error is calculated. As will
be discussed, it may be convenient to choose ψ(·) bounded
with |ψ ′ (s)| −−−−→ 0, for instance, ψ(s) = arctan sp or
|s|→∞

r̄∗

ψ(s) = |s|p sign s/(1 + |s|p ) with p ≥ 1. The mentioned two
functions, as well as the simplest function ψ(s) = s satisfy
the condition, which is henceforth assumed to be valid

ϕ(x, y) = const < 0

−kn en̄

As illustrated in Fig. 2, the plane R2 is covered by the
disjoint level sets of the function ϕ, namely the sets where
ϕ(x, y) = c = const. If (3) holds, then the vector n̄(x, y) is
the normal vector to the corresponding level set at the point
(x, y). The path P is one of the level sets, corresponding to
c = 0; the value ϕ(x(t), y(t)) can be considered as a (signed)
“distance” from the robot to the path (differing, as usual, from
the Euclidean distance), or the tracking error [17], [22]. More
generally, choosing an arbitrary strictly increasing function
ψ ∈ C 1 (R → R) with ψ(0) = 0 (and thus ψ(s)s > 0 for
any s 6= 0), one may define the tracking error as follows

sup ψ ′ (ψ −1 (u)) < ∞.

(5)

u∈R

x∗

0

C. Technical assumptions

X

Henceforth the following three technical assumptions are
adopted, excluding some “pathological” situations. Our first
assumption enables one to use of the tracking error e(x, y) as a
“signed distance” to the path P and implies that the asymptotic
vanishing of the error e(r̄(t)) −−−→ 0 entails the convergence

Fig. 2. The robot orientation and level sets of the function ϕ(x, y)

B. The description of the desired path
The desired curvilinear path P is the zero set of a function
ϕ ∈ C 2 (R2 → R), i.e. it is described by the implicit equation
∆

P = {(x, y) : ϕ(x, y) = 0} ⊂ R2 .

(2)

The same curve P may be represented in the implicit form
(2) in many ways. The principal restrictions imposed by our
approach is regularity [17]: in some vicinity of P one has
i⊤
h
∆
∂ϕ(x,y)
6= 0. (3)
n̄(x, y) = ∇ϕ(x, y) = ∂ϕ(x,y)
;
∂x
∂y

t→∞

to the path in the usual Euclidean metric1 dist(r̄(t), P) → 0.
Assumption 1. For an arbitrary constant κ > 0 one has
inf{|e(r̄)| : dist(r̄, P) ≥ κ} > 0.

(6)

Our second assumption provides the regularity condition (3)
in a sufficiently small vicinity of P.
1 As usual, the distance from a point r̄ to the path P A is dist(r̄ , P) =
0
0
inf{|r̄0 − r̄| : r̄ ∈ P}. More generally, the distance between sets A, B is
defined as dist(A, B) = inf{|r̄1 − r̄2 | : r̄1 ∈ A, r̄2 ∈ B}.

4

Assumption 2. The set of critical points, where n̄ vanishes,
∆

2

C0 = {(x, y) ∈ R : ∇ϕ(x, y) = 0},
is separated from P by a positive distance dist(C0 , P) > 0.
In the most typical cases where either P is a closed curve
or C0 is compact, Assumption 2 boils down to the condition
C0 ∩ P = ∅. Our final assumption is similar in spirit to
Assumption 1 and guarantees that the asymptotic vanishing
of the normal vector n̄(r(t)) −−−→ 0 is possible only along a
t→∞
trajectory, converging to C0 , i.e. dist(r̄(t), C0 ) → 0.
Assumption 3. For an arbitrary constant κ > 0 one has
inf{|n̄(r̄)| : dist(r̄, C0 ) ≥ κ} > 0.

(7)

Assumption 3 implies the following useful technical lemma.
Lemma 1. Consider a Lipschitz vector-function r̄ : [0; ∞) →
R2 such that r̄(t) does not converge to C0 as t → ∞,
that is, lim sup d(t) > 0 where d(t) = dist(r̄(t), C0 ). Then
t→∞
R∞
p
|n̄(r̄(t))|
dt = ∞ for any p > 0.
0

Proof: It can be easily noticed that d(·) is a Lipschitz
function, since as can be easily shown, |d(t1 ) − d(t2 )| ≤
|r̄(t1 ) − r̄(t2 )| ≤ M |t1 − t2 | ∀t1 , t2 ≥ 0, where M > 0 is the
Lipschitz constant for r̄(·). Since lim sup d(t) > 0, a number
t→∞

ε > 0 and a sequence tk → ∞ exist such that d(tk ) ≥ 2ε and,
therefore, d(t) ≥ ε for t ∈ ∆k = [tk ; tk + M −1 ε]. Passing
to subsequences, one may assume without loss of generality
that tk + M −1 ε < tk+1 and thus the intervals ∆k are disjoint.
Assumption 3 implies that for sufficiently small
c > 0 one has
R
|n̄(r̄(t))| ≥ c whenever t ∈ ∆k . Hence ∆k |n̄(r̄(t))|p dt ≥
R∞
εM −1 cp and therefore 0 |n̄(r̄(t))|p dt = ∞.
III. T HE G UIDING V ECTOR F IELD AND I TS P ROPERTIES
In this section, we construct the Guiding Vector Field (GVF)
to be used in the path following control algorithm. We show
that the integral curves of this field lead either to the desired
path P or to the critical set C0 . Furthermore, we give efficient
criteria, ensuring that the integral curves of the second type
are either absent or cover a set of zero measure on R2 .
Besides the normal vector n̄(x, y) = ∇ϕ(x, y), at each point
we consider the tangent vector to the level set


0 1
τ̄ (x, y) = E n̄(x, y), E =
.
(8)
−1 0

If the point is regular (3), the basis (τ̄ , n̄) is right-handed
oriented2 (Fig. 2). Our goal is to find such a vector field
v̄(x, y), where the absolute tracking error |e| is decreasing
along each of its integral curves (unless e = 0), and the curves
starting on P do not leave it. We define this vector field by
∆

v̄(x, y) = τ̄ (x, y) − kn e(x, y)n̄(x, y).

(9)

The integral curves of the vector field (9) correspond to the
trajectories of the autonomous differential equation
d
¯
ξ̄(t) = v̄(ξ(t))
∈ R2 , t ≥ 0.
(10)
dt
2 In the subsequent constructions, one may replace τ̄ with −τ̄ and E to
−E, which causes the change of the path following direction.

A. Properties of the integral curves
We notice first that the vector field v̄ is C 1 -smooth, which
implies the local existence and uniqueness of solutions of the
system (10). Obviously, any point (x0 , y0 ) ∈ C0 corresponds
to the equilibrium of (10), being a trivial single-point integral
curve. The uniqueness property implies that non-constant
integral curves are free of the critical points; the corresponding
solutions, however, may converge to C0 asymptotically. In
general, a solution may also escape to infinity in finite time.
The following lemma establishes the principal dichotomy
property of the integral curves, stating that any curve leads
either to the desired path P or to the critical set C0 .
Lemma 2. Let ξ̄(t), t ∈ [0; t∗ ) be a maximally prolonged
solution to (10). Then two situations are possible
¯ P) −−−→ 0, that is, the solution con1) either dist(ξ(t),
t→t∗
verges to the desired path,
¯ C0 ) −−−→ 0.
2) or t∗ = ∞ and dist(ξ(t),
t→∞

Proof: The proof is based on the Lyapunov function
V (x, y) =

1
e(x, y)2 ≥ 0.
2

(11)

¯
A straightforward computation shows that ∇e(ξ)
=
¯ ∇ϕ(ξ)
¯ = ψ ′ (ψ −1 (e(ξ)))|n̄(
¯
¯ and hence V is nonψ ′ (ϕ(ξ))
ξ)|
¯ since its derivative is
increasing along the trajectory ξ(t)
⊤ ¯
¯ = e(ξ)ψ
¯ ′ (ψ −1 (e(ξ)))n̄
¯
¯ (9)
V̇ (ξ)
(ξ)v̄(ξ)
=
(9)

¯ 2 ψ ′ (ψ −1 (e(ξ)))|n̄(
¯
¯ 2 ≤ 0.
= −kn e(ξ)
ξ)|

(12)

¯
≥ 0.
In particular, there exists the limit e∗ = lim |e(ξ(t))|
t→t∗
If e∗ = 0 then, due to Assumption 1, statement 1) holds:
the solution converges to P. Suppose now that e∗ > 0. We
are going to prove that statement 2) is valid. Notice first
¯
that ψ ′ (ψ −1 (e(ξ(t)))
is uniformly bounded and positive; the
¯
same
holds
for
V
(
ξ(t)).
The equality (12) thus implies that
R t∗
2
¯
|n̄(
ξ(t))|
dt
<
∞.
Since
|v|2 = (1 + kn2 e2 )|n̄|2 , (10)
0
R t∗ ˙
2
¯
dt < ∞. If one had t∗ < ∞, the
implies that 0 |ξ(t))|
Cauchy-Schwartz inequality would imply that
!2
Z T
Z t∗
2
˙
¯
¯
¯
¯˙ 2 dt ∀T < t ,
|ξ(T )−ξ(0)| =
|ξ(t)|dt
|ξ(t)|
≤t
∗

∗

0

0

which contradicts the assumption that ξ̄(t) is a maximally
prolonged solution, escaping to infinity as t → t∗ . Therefore,
t∗ = ∞; statement 2 now follows from Lemma 1.
A natural question arises on how “large” the set of trajectories is, converging to the critical set C0 . In many practical examples, this set is finite. This holds, in particular if ϕ(x, y) 6= 0
when |x| + |y| is sufficiently large and its Hessian matrix
#
" 2
∂2
∂
⊤
∂x2 ϕ(x, y)
∂x∂y ϕ(x, y)
(13)
H(x, y) = H(x, y) = ∂ 2
∂2
∂x∂y ϕ(x, y)
∂y 2 ϕ(x, y)
is sign-definite at
H(x0 , y0 ) > 0 or
bounded (and thus
which implies that

any point (x0 , y0 ) ∈ C0 , i.e. either
H(x0 , y0 ) < 0. In this situation C0 is
compact) and all its points are isolated,
C0 is finite. If the set C0 is finite, the

5

¯ C0 ) −−−→ 0, obviously, means that the
convergence dist(ξ(t),
t→∞

∆

solution converges to a critical point: ξ¯∗ = lim ξ̄(t) ∈ C0 .
t→∞

Definition 1. Let ξ̄∗ ∈ C0 be an equilibrium point of (10).
The stable manifold of ξ¯∗ , denoted by W (ξ¯∗ ), is the set of all
points ξ̄0 such that the solution of (10), starting at ξ(0) = ξ0 ,
∆
exists for all t ≥ 0 and ξ¯∗ = lim ξ̄(t).
t→∞

We will use the following corollary of the central manifold
theorem (see e.g. Theorem 4.1 and Proposition 4.1 in [41]).
Lemma 3. If both eigenvalues of the Jacobian matrix J(ξ¯∗ ) =
∂
v̄(ξ¯∗ ) are strict unstable Re λ1,2 J(ξ¯∗ ) > 0, then W (ξ¯∗ ) =
∂ ξ̄
¯
{ξ∗ }. If J(ξ¯∗ ) has at least one strictly unstable eigenvalue,
¯ is a set of zero measure.
then W (ξ)
Lemma 3 in turn has the following important corollary.
Corollary 1. Let the set C0 be finite and for any of its points
ξ¯∗ ∈ C0 the matrix e(ξ¯∗ )H(ξ¯∗ ) has a negative3 eigenvalue.
Then the maximally prolonged solution of (10) (possibly,
existing on finite interval only) converges to P for almost all
initial conditions ξ̄(0). If e(ξ¯∗ )H(ξ¯∗ ) < 0 for any ξ∗ ∈ C0 ,
¯ 6∈ C0 .
this convergence takes place whenever ξ(0)
Proof: A straightforward computation shows that if
∇ϕ(ξ¯∗ ) = 0 then J(ξ¯∗ ) = (E − kn e(ξ¯∗ ))H(ξ¯∗ ) and
Tr J(ξ¯∗ ) = −kn e(ξ¯∗ ) Tr H(ξ¯∗ ),
det J(ξ¯∗ ) = (1 + k 2 e(ξ¯∗ )2 ) det H(ξ¯∗ ).
n

It can be easily shown that a symmetric 2 × 2 matrix M
is non-negatively definite M ≥ 0 if and only det M =
λ1 (M )λ2 (M ) ≥ 0 and Tr M = λ1 (M ) + λ2 (M ) ≥ 0.
Here λi (M ) stand for the eigenvalues of M . Since eH(e)
is not non-negatively definite, either det J < 0 or Tr J >
0, det J ≥ 0. In both situations, the matrix J(ξ¯∗ ) has at least
one strictly unstable eigenvalue. Furthermore, if eH(e) < 0
then Tr J > 0, det J > 0, i.e. both eigenvalues of J are strictly
unstable. The statement now follows from Lemma 3.
B. The GVF and the “ideal” motion of the robot
The main idea of the path following controller, designed in
the next section, is to steer the robot to the integral curve of
the field (9). In other words, when the robot is passing a point
r̄ = (x, y) ∈ R2 , its desired orientation is
m̄d (x, y) =

1
v̄(x, y).
|v̄(x, y)|

(14)

The field of unit vectors m̄d (x, y), henceforth referred to as
the guiding vector field (GVF), is defined
at any regular point
p
(x, y), where n̄ 6= 0 and thus |v̄| = 1 + kn2 e2 |n̄| 6= 0. Fig. 2
illustrates the relation between the vectors r̄, m̄, τ̄ , n̄, v̄, m̄d .
Consider the desired motion of the robot, “ideally” oriented
along the integral curves of the GVF at any point. Its position
vector r̄(t) obeys the differential equation
˙ = ur m̄d (r̄(t)),
r̄(t)
3 Recall

t ≥ 0.

that H is a symmetric matrix, so its eigenvalues are real

(15)

The following lemma is dual to Lemma 2.
Lemma 4. Let r̄(t), t ∈ [0; t∗ ) be a maximally prolonged
solution to (15). Then two situations are possible
1) either dist(r̄(t), C0 ) −−−→ 0,
t→t∗

2) or t∗ = ∞ and dist(r̄(t), P) −−−→ 0 that is, the solution
t→∞
converges to the desired path;
R t∗ Proof: If t∗ < ∞ then the limit r̄∗ = r̄(t∗ − 0) = r̄(0) +
m̄d (r̄(t))dt exists and, since the solution is maximally
0
prolonged, one obviously has |n̄(r̄∗ )| = 0, i.e. r̄∗ ∈ C0
and statement 1 holds. The case of t∗ = ∞ is considered
similar to the proof of Lemma 2. Introducing the Lyapunov
function (11), its derivative is shown to be
(9)

V̇ (r̄) = ur e(r̄)ψ ′ (ψ −1 (e(r̄)))n̄⊤ (r̄)m̄d (r̄) =
ur kn e(r̄)2
(9)
=−p
ψ ′ (ψ −1 (e(r̄))) |n̄(r̄)| ≤ 0.
1 + kn2 e(r̄)2
{z
}
|

(16)

Ψ(e(r̄))

There exists the limit e∗ = limt→∞ |e(r̄(t))|. If e∗ = 0,
statement 2 holds thanks to Assumption 1. Otherwise, e∗ > 0
Ψ(e(r̄(t))) is uniformly positive. In view of (16), one has
Rand
∞
|n̄(r̄(t))|dt
< ∞. Lemma 1 entails now statement 1.
0

Remark 1. The equation (16) implies that V̇ ≈ −kn V |n̄|θ(e),
where θ(e) → θ0 > 0 as e → 0. Assumptions 2 and 3 imply
that |n̄(r)| ≥ κ as e ≈ 0. Therefore, the desired path P
is locally exponentially attractive: if |e(0)| ≤ ε, where ε is
sufficiently small, then |e(t)| ≤ e−kn ur βt |e(0)|, where β > 0
is a constant, depending on ψ(·) and the normal vector n̄(x, y)
in the vicinity of P. The coefficient kn is responsible for the
attraction to P. In the limit case kn = 0 the path P is not
asymptotically stable; the higher kn > 0 one chooses, the
stronger is the attraction to the path P in its small vicinity.
Obviously, the integral curves of (15) are in one-to-one
correspondence with non-equilibrium integral curves of (10).
Corollary 1 can now be reformulated as follows.
Corollary 2. Let the set C0 be finite and for any of its points
ξ¯∗ ∈ C0 the matrix e(ξ¯∗ )H(ξ¯∗ ) has a negative eigenvalue.
Then for almost all initial conditions r̄(0) the solutions of (15)
can be prolonged up to ∞ and converge to P. If e(ξ¯∗ )H(ξ¯∗ ) <
0 for any ξ∗ ∈ C0 , this holds for any r̄(0) 6∈ C0 .
Note that the conditions of Corollary 2 always hold for
strictly convex function ϕ(x, y) (that is, H(x, y) > 0 at
any point) such that since the critical point (if it exists) is
the unique point of global minimum. Therefore, at this point
(x∗ , y∗ ) one has ϕ(x∗ , y∗ ) < 0 and thus e < 0, which implies
that eH < 0. This is exemplified by the function
(y − y0 )2
(x − x0 )2
+
− 1,
a2
b2
defining the elliptic path. Fig. 3 demonstrates the corresponding GVF with the unique critical point at the center of ellipse.
This critical point is “repulsive” and no trajectory of (15)
converges to it. Fig. 4 illustrates the GVF for the Cassini oval,
being the zero level set of the non-convex function
ϕ(x, y) =

(x − x0 )4 + (y − y0 )4 − 2a2 ((x − x0 )2 − (y − y0 )2 ) + a4 − b4 .

6

˙ d (t)⊤ E m̄d (t).
multiplier ωd (t) can be found from ωd (t) = −m̄
It remains to prove that ωd (t) in fact depends only on the trajectory (x(t), y(t), α(t)), and this dependence is continuous.
Introducing the vector field (9) v̄(t) = v̄(x(t), y(t)) and the
tracking error e(t) = e(x(t), y(t)) along the robot’s trajectory,
a straightforward computation shows that


I2
v̄v̄ ⊤
d v̄
˙
˙d=
v̄,
=
−
m̄
dt kv̄k
kv̄k kv̄k3
(18)
v̄˙ = ur [E − kn eI2 ]H(x, y)m̄(α) − kn ė n̄(x, y),

700

600

500

400

300

200

100

ė = ur ψ ′ (ϕ(x, y)) n̄(x, y)⊤ m̄(α).

0
0

200

400

600

800

1000

1200

Here I2 is 2 × 2 identity matrix and H is the Hessian (13).
˙ d and m̄d continuously depend on the triple
Since ϕ ∈ C 2 , m̄
˙⊤
(x(t), y(t), α(t)), the same holds for ωd = −m̄
d E m̄d .
To clarify the meaning of the function ωd , suppose for the
moment that the robot’s speed is ur = 1. If the robot is moving
strictly along the integral curves of the GVF, then ωd is the
signed curvature of the robot’s trajectory at its current position.
In general, ωd can be treated as the “desired” curvature of the
robot’s trajectory, which may differ from its real curvature.

Fig. 3. The GVF for an elliptic path

700

600

500

400

300

200

B. The algorithm of GVF steering

100

As was discussed in the foregoing, the idea of the path
following algorithm is to steer the robot’s orientation along the
guiding vector field m̄d = m̄d (x, y). We introduce the directed
angle δ = δ(x, y, α) ∈ (−π; π] between m̄d and m̄ (see
Fig. 2). The function δ is thus defined at any point (x, y, α)
and C 1 -smooth at the points where m̄(α) 6= −m̄d (x, y).
The orientation vector’s derivative along the trajectory is


d
− sin α(t)
˙
= −ω(t)E m̄(t).
m̄(t)
= m̄(α(t)) = ω(t)
cos α(t)
dt
(19)
On the other hand, m̄ may be decomposed (Fig. 2) as

0
0

200

400

600

800

1000

1200

Fig. 4. The GVF for a Cassini oval.

The corresponding set C0 consists of two “locus points”
(x0 ± a, y0 ) and the “center” (x0 , y0 ). In all these points
e < 0. At the locus points one has H > 0, whereas at the
center H has one positive and one negative eigenvalue, so the
robot potentially can be “trapped” at the center, but the set of
corresponding initial conditions has zero measure.
IV. T HE PATH F OLLOWING C ONTROLLER D ESIGN
In this section, we design the controller, steering the robot
to the GVF (14). This controller uses the GVF m̄d (x, y) at
the current robot’s position and the function ωd (x, y, α), measuring the GVF’s “rotation rate” along the robot’s trajectory.
A. Preliminaries
We start with the following technical lemma.

(14)
˙ = −δ̇ [(sin δ)I2 + (cos δ)E]m̄d +
−ωE m̄ = m̄
|
{z
}
=E m̄

˙ d (17)
+ [(cos δ)I2 − (sin δ)E]m̄
= −δ̇E m̄−

− ωd [(cos δ)I2 − (sin δ)E]E m̄d = −(δ̇ + ωd )E m̄,
|
{z
}
=E m̄

Lemma 5. Let m̄d (t) = m̄d (x(t), y(t)) stand for the GVF
˙ d (t) is
along a trajectory of the robot. Then its derivative m̄
ṁd (t) = −ωd (x(t), y(t), α(t))E m̄d (t),

m̄ = (cos δ)m̄d − (sin δ)E m̄d = [(cos δ)I2 − (sin δ)E]m̄d ,
(20)
At any point where δ(t) < π and thus δ̇(t) exists, one obtains

entailing the following principal relation between δ, ω and ωd
δ̇ = ω − ωd .

(21)

(17)

where ωd : R3 → R is continuous and uniquely determined
by the functions ϕ(·), ψ(·) and the constant kn from (9).
Proof: Differentiating the equality |m̄d (t)|2 = 1, one
˙ d (t)⊤ m̄d (t) = 0, that is, m̄
˙ d (t) ⊥ m̄d (t).
shows that m̄
˙ d (t) is proportional to the unit vector
Therefore, the vector m̄
˙ d (t) = −ωd (t)E m̄d (t) and the scalar
Emd (t), so that m̄

Furthermore, at any time the following equality is valid
d
d
sin δ = − m̄⊤ E m̄d = (ω − ωd ) cos δ.
dt
dt

(22)

Supposing that δ(t0 ) = π and ω(t) − ωd (t) ≤ 0 for t ≈
t0 , equation (22) entails that sin δ ≥ 0 and hence δ = π2 +
arcsin(sin δ) for t ∈ [t0 ; t0 + ε), where ε > 0 is sufficiently

7

small. In this situation, the function δ(x(t), y(t), α(t)) has the
right derivative4 δ̇ = D+ δ, satisfying (21) as t ∈ [t0 ; t0 + ε).
We are now in a position to describe our path-following
algorithm, employing the GVF “rotation rate” ωd and the angle
of discrepancy between the GVF and the robot’s orientation δ
ω(t) = ωd (x(t), y(t), α(t)) − kδ δ(x(t), y(t), α(t)).

(23)

Here kδ > 0 is a constant, determining the convergence rate.
When δ(x(t), y(t), α(t)) < π, equality (21) holds and thus
δ̇ = ω − ωd = −kδ δ.

(24)

Furthermore, even for δ(x(0), y(0), α(0)) = π one has
ω − ωd = −kδ δ < 0 as t ≈ t0 and hence (24) retains its
validity at t = 0, treating δ̇ as the right derivative D+ f .
Thus, considering δ̇ as a new control input, the algorithm (21)
is equivalent to a very simple proportional controller (24),
providing, in particular, that δ(x(t), y(t), α(t)) < π ∀t > 0.
C. Local existence and convergence of the solutions
In this subsection we examine the properties of the solutions
of the closed-loop system (1), (23), rewritten as follows
ẋ(t) = ur cos α(t),

ẏ(t) = ur sin α(t),

(25)

α̇(t) = ωd (x(t), y(t), α(t)) − kδ δ(x(t), y(t), α(t)).

The right-hand side of (25) is continuous at any point
(x0 , y0 , α0 ), where n̄(x0 , y0 ) 6= 0 and δ0 = δ(x0 , y0 , α0 ) < π.
However, the discontinuity at the points where δ0 = π makes
the usual existence theorem [42] inapplicable. To avoid this
problem, we consider the equivalent “augmented” system
ẋ(t) = ur cos α(t),

ẏ(t) = ur sin α(t),

α̇(t) = ωd (x(t), y(t), α(t)) − kδ δ∗ (t),

(26)

δ̇∗ (t) = −kδ δ∗ (t).
As was discussed in the previous subsection, any solution
of (25) satisfies (26) with δ∗ (t) = δ(x(t), y(t), α(t)), and
vice versa: choosing a solution (x(t), y(t), α(t), δ∗ (t)), where
δ∗ (0) = δ(x(0), y(0), α(0)) ∈ (−π; π], one has δ∗ (t) =
δ(x(t), y(t), z(t)) for any t ≥ 0 due to (24). Unlike (25), the
right-hand side of (26) is a C 1 -smooth function of (x, y, α, δ∗ )
at any point where n̄(x, y) 6= 0. The standard existence and
uniqueness theorem [42] implies the following lemma.
Lemma 6. For any point ζ0 = (x0 , y0 , α0 ), such that
n̄(x0 , y0 ) 6= 0, there exists the unique solution ζ(t) =
(x(t), y(t), α(t)) with ζ(0) = ζ0 . Extending this solution to
the maximal existence interval [0; t∗ ), one either has t∗ = ∞
or t∗ < ∞ and (x(t), y(t)) −−−→ (x∗ , y∗ ) ∈ C0 .
t→t∗

Proof: Reducing the Cauchy problem for the closed-loop
system (25) to the Cauchy problem for (26), one shows that the
solution exists locally and is unique [42]. Let its maximally
∆
prolongable solution (x(t), y(t), α(t)) be defined on ∆∗ =
˙
[0; t∗ ) with t∗ < ∞. Since |r̄(t)|
= ur , the limit exists
Z t∗
˙
r̄(t)dt.
r̄∗ = lim r̄(t) = r̄(0) +
t→t∗

0

4 By

definition, the right derivative of a function f (t) at t = t0 (written
(t0 )
.
f ′ (t0 + 0) or D + f (t0 )) is defined by D + f (t0 ) = lim f (t)−f
t−t
t→t0 +0

0

We are going to show that r̄∗ ∈ C0 , i.e. n̄(r̄∗ ) = 0. Suppose,
on the contrary, that |n̄(r̄∗ )| > 0; therefore, |v̄(x(t), y(t))|
is uniformly positive on ∆∗ . Using (18), this implies that
˙ d (t), and hence ωd (x(t), y(t), α(t)) and ω(t) are uniformly
m̄
bounded on ∆∗ . Thus there exists the finite limit
Z t∗
ω(t)dt,
α∗ = lim α(t) = α(0) +
t→t∗

0

enabling one to define the solution at t = t∗ and then to
prolong it to [t∗ , t∗ + ε), i.e. the solution is not maximally
prolonged. This contradiction proves that n̄(r̄∗ ) = 0.
Note that if δ(0) = 0, that is, the robot was perfectly
oriented along the GVF at the starting moment, (21) implies
that δ(t) ≡ 0 so that the robot follows the integral curve of
the GVF. As was shown in Section III, in this “ideal” situation
the robot either approaches the desired path P or is driven to
one of the critical points. The latter situation is practically
impossible if the conditions of Corollary 2 are valid since the
set of integral curves, leading to the set C0 , has zero measure.
One may consider (25) as a system with slow-fast dynamics. Informally, the controller (24) provides the exponential
convergence of the robot to an integral curve of the GVF;
after this “fast” transient process, the robot “slowly” follows
this integral curve and approaches the desired trajectory,
unless it is “trapped” in a critical point. Ignoring the “fast
dynamics”, one may suppose that the statement of Lemma 4
remains valid for a general solution of the system (25). This
argument, however, is not mathematically rigorous. Recalling
the proofs of Lemmas 2 and 4, one may notice that the central
argument was the non-increasing property of the Lyapunov
function (11). Although the deviation of the robot from the
integral curve exponentially decreases due to (24), the nonincreasing property, in general, fails; as will be shown below,
even if the robot is positioned very close to the desired path,
the tracking error may increase. However, this effect does not
destroy the dichotomy property (any solution converges either
to P or to C0 ) under the following assumption, which usually
holds in practice, being valid e.g. if ϕ(·) is a polynomial.
Assumption 4. There exist θ ∈ (0; kδ ) and C > 0 such that
|n̄(r̄)| ≤ Ceθ|r̄|

as |r̄| → ∞.

(27)

The latter condition can be relaxed in the case of a closed
path; however, we adopt it to consider both bounded and
unbounded paths in a unified way. Note that if ϕ(x, y) is a
polynomial function, then kδ > 0 can be arbitrarily small.
Henceforth all Assumptions 1-4 are supposed to be valid.
The following theorem, extending Lemma 4 to the case
where δ(0) 6= 0, is our first main result.
Theorem 1. For any maximally prolonged solution
(x(t), y(t), α(t)) of (25), defined for t ∈ [0; t∗ ), one of
the following statements holds:
1) either t∗ = ∞ and dist(r̄(t), P) −−−→ 0, or
t→∞
2) dist(r̄(t), C0 ) → 0 as t → t∗ .
In other words, for any initial condition the algorithm drives
the robot to either the desired path P or C0 .

8

Proof: In the case where t∗ = ∞ statement 2 holds
due to Lemma 6. Suppose that t∗ = ∞. Differentiating the
function (11) along the trajectories, it can be shown that
(20)

V̇ = ur eψ ′ (ψ −1 (e))n̄⊤ m̄ =
(20)

(14)

= ur eψ ′ (ψ −1 (e))n̄⊤ [m̄d cos δ − E m̄d sin δ] =
′
−1
(e)) ⊤
(9)
(14) ur eψ (ψ
n̄ [v̄ cos δ − Ev̄ sin δ] =
=
|v|
′
−1
(e))
(9) ur eψ (ψ
= p
|n̄|(−kn e cos δ + sin δ) =
2
1 + kn e2
= Φ(e)|n̄|(−kn e cos δ + sin δ).

(28)

ec = inf{|e(x, y)| : (x, y) ∈ C0 } > 0.

∆

ur eψ ′ (ψ −1 (e))
p
.
1 + kn2 e2

(30)

Consider the following set

Here Φ(e) denotes the bounded, in view of (5), function
Φ(e) =

criterion requires the initial condition (x(0), y(0), α(0)) to
belong to some invariant set, free of critical points. Similar
restrictions arise in most of the path following algorithms; for
example, in the projection-based algorithms the convergence
can be rigorously proved only in some region of attraction
where the projection to the desired curve is well defined [14].
Assumptions 1 and 2 imply the uniform positivity of the
error on C0 , that is, the following inequality holds

(29)

Since
| sin δ|
≤
|δ|, Assumption 4 entails that
R∞
|Φ(e(t))n̄(t)
sin
δ(t)|
dt
< ∞. Notice now that
0
(−Φ(e)e) ≤ 0 and cos δ(t) >
R ∞0 as t becomes sufficiently
large. Thus the integral I = 0 (−eΦ(e)|n̄|) cos δ dt exists,
being either finite or equal
R ∞to −∞. This implies, thanks
to (28), the existence of 0 V̇ dt = lim V (t) − V (0).
t→+∞
Since V ≥ 0, one has I > −∞ and therefore there exists the
limit e∗ = limt→∞ |e(x(t), y(t))|. If e∗ = 0, then statement 1
holds due to Assumption 1. Otherwise, eΦ(e) is uniformly
positive and thus,R recalling that cos δ(t) → 1 as t → ∞,
∞
one obtains that 0 |n̄(x(t), y(t))|dt < ∞, which implies
statement 2 thanks to Lemma 1.
Corollary 3. If C0 = ∅, then for any initial condition the
solution of (25) is infinitely prolongable and the algorithm
solves the path following problem dist(r̄(t), P) −−−→ 0.
t→∞

Corollary 3 is applicable to linear mappings ϕ(x, y) =
ax + by + c (with |a| + |b| 6= 0) and many other functions,
e.g. ϕ(x, y) = y + f (x). These functions, however, usually
correspond to unbounded desired curves, whereas for closed
paths the GVF m̄d is usually not globally defined.
The experiments show that under the assumptions of Corollary 2 the robot always “evades” the finite set of critical
points and converges to the desired trajectory. This looks
very natural since after very fast transient dynamics the robot
“almost precisely” follows some integral curve, which leads
to P “almost surely”. We formulate the following hypothesis.
Hypothesis. Under the assumptions of Corollary 2, for
almost all initial conditions (x(0), y(0), α(0)) the robot’s
trajectory (x(t), y(t)) converges to the desired path P.
Whereas the proof of this hypothesis remains a challenging
problem, it is possible to guarantee the global existence of
the solutions and their convergence to the desired path in
some broad invariant set, free of the critical points. The
corresponding result, which does not rely on the assumptions
of Corollary 2, is established in the next subsection.
D. An invariant set, free of critical points
In this subsection we give a sufficient condition, guaranteeing that a solution of (25) does not converge to C0 . This

M = {(x, y, α) : n̄ 6= 0, |δ| < arctan(kn ec ), |e| < ec } .
(31)
Recall that kn is the constant parameter from (9) and δ =
δ(x, y, α) is the angle between the robot’s heading and the
vector field direction (see Fig. 2). By definition, M ∩ C0 = ∅.
The following lemma states that in fact M is an invariant set,
i.e. any solution starting in M remains there.
Theorem 2. Any solution of (25), starting at
(x(0), y(0), α(0)) ∈ M, does not leave M, and is infinitely
prolongable satisfying the inequality
|e(t)| ≤ max{|e(0)|, kn−1 tan δ(0)} < ec .

(32)

For such a solution, one has dist(r̄(t), P) −−−→ 0, i.e. the
t→∞
algorithm (23) solves the path following problem in M.
Proof: Consider a solution (x(t), y(t), α(t)), starting at
(x(0), y(0), α(0)) ∈ M. Due to (24), one has |δ(t)| ≤
|δ(0)| < π2 ∀t ≥ 0, and hence cos δ(t) > 0. By noticing that
Φ(e)e ≥ 0 and thus |Φ(e)| = Φ(e) sign e, one has
(28)

−V̇ = −ur |Φ(e)| sign e |n̄|(e cos δ + sin δ) =
= ur |Φ(e)||n̄|(|e| + sign e tan δ) cos δ ≥
≥ ur |Φ(e)| |n̄|(|e| − | tan δ|) cos δ.
In particular, V̇ ≤ 0 whenever |e| ≥ | tan δ|.
Notice that if |e(t)| ≥ | tan δ(t)| for any t (where the solution exists), then, obviously |e(t)| ≤ |e(0)| and thus (32) holds.
Suppose now that |e(t0 )| < | tan δ(t0 )| at some t0 ≥ 0. We are
going to show that |e(t)| ≤ | tan δ(t0 )| ≤ | tan δ(0)| for any
t ≥ t0 . Indeed, had we |e(t1 )| > | tan δ(t0 )| at some point t1 ,
there would exist t∗ < t1 such that |e(t∗ )| = | tan δ(t0 )| and
|e(t)| > | tan δ(t0 )| as t ∈ (t∗ ; t1 ]. Using (24), it can be easily
shown that | tan δ(t0 )| ≥ | tan δ(t)| for t ≥ t0 , and hence
|e(t)| > | tan δ(t)| is non-increasing when t ∈ (t∗ ; t1 ], which
contradicts the assumption that |e(t1 )| > | tan δ(t0 )| = |e(t∗ )|.
We have proved that in both cases 1) and 2) the inequality (32)
holds at any point where the solution exists; thus the solution
stays in M. By definition (30) of ec , the vector r̄(t) cannot
converge to C0 in finite or infinite time, i.e. for the considered
solution the statement in Theorem 1 holds.
Remark 2. The condition (x(0), y(0), α(0)) ∈ M restricts
the robot to be “properly” headed in the sense that
|δ(x(0), y(0), α(0))| < arctan(kn ec ) < π/2.

(33)

Since ec > 0, (33) is valid for sufficiently large kn whenever
|δ| < π/2. In other words, Theorem 2 guarantees convergence

9

to the path from any starting position with |e(0)| < ec and
|δ(0)| < π/2, choosing large kn . Furthermore, if the desired
path P is closed and the direction of circulation along it is
unimportant, the condition |δ(0)| < π/2 can be provided by
reverting the vector field m̄d 7→ −m̄d (which corresponds
to the replacement of ϕ 7→ −ϕ and |δ| 7→ π − |δ|) unless
δ(0) = ±π/2 (in practice it is never possible).
Remark 3. Even if (33) is violated at the starting time
t = 0, it obviously holds when t > t0 = kδ−1 [ln |δ(0)| −
ln arctan(kn ec )] thanks to (24). Thus, if one is able
to prove that the solution is prolongable up to t0 and
|e(x(t0 ), y(t0 ))| < ec , Theorem 2 provides the convergence
of the robot’s position to the desired path.
Remark 3 suggests the way to relax the restriction on the
initial robot’s orientation (33). Since the robot moves at the
constant speed ur > 0, it covers the distance ur t0 until its
orientation satisfies (33). The path following is guaranteed if
ur t0 is less than the viability distance of the initial position,
i.e. the distance from it to the set where |e| ≥ ec .
Definition 2. Given a point r̄0 = (x0 , y0 ) with |e(r̄0 )| < ec ,
the number d0 = inf{dist(r̄0 , r̄) : e(r̄) ≥ ec } > 0 is said to
be its viability distance.
Theorem 2 and Remark 3 yield in the following.
Corollary 4. Let the initial position of the robot r̄(0) =
(x(0), y(0)) with e(r̄(0)) < ec have the viability distance
d0 > 0. If this viability distance satisfies the condition
ur |δ(x(0), y(0), α(0))|
ln
,
(34)
kδ
arctan(kn ec )
then the solution of (25) gets into the set M in finite time,
and is prolongable up to ∞ satisfying dist(r̄(t), P) −−−→ 0.
d0 >

t→∞

Note that (34) holds for any α(0) under the assumption
ur
π
d0 >
ln
.
(35)
kδ
arctan(kn ec )
The condition (35) gives an estimate for the region, starting in
which the robot necessarily converges to the desired path P.
Taking the original orientation relative to the field into account,
this estimate can be tightened by using (34). Typically d0 is
uniformly positive in the vicinity of P (this holds, for instance,
when P is a closed curve). Thus (35) guarantees the algorithm
to converge in the vicinity of P, choosing kδ /ur sufficiently
large. As was mentioned, practical experiments with natural
trajectories, satisfying the assumptions of Corollary 2, show
that the robot is always attracted to the desired path, although
the mathematical proof of this remains a non-trivial task.
Remark 4. Although explicit calculation of the viability distance is complicated, conditions (34) and (35) in fact require
only its lower bound. Such a bound can be explicitly obtained,
for instance, if the error is Lipschitz
|∇e(r̄)| = ψ ′ (ϕ(r̄))|n̄(r̄)| ≤ c = const

∀x, y.

(36)

If the condition (36) holds, then |e(r̄) − e(r̄0 )| ≤ c|r̄ − r̄0 | and
hence the viability distance of r̄0 is estimated as follows
d0 ≥ (ec − e(x0 , y0 ))/c.

Fig. 5. The E-puck robot with marker on top, used in experiments

Condition (36) can often be provided under an appropriate
choice of ψ(·). For instance, when ϕ(x, y) = (x − x0 )2 + (y −
y0 )2 − R2 (circular path), one can choose ψ(s) = arctan s so
that ψ ′ (s) = 1/(1 + s2 ). More generally, if |ϕ(r̄)| ≥ C1 |r̄|β
and |∇ϕ(r)| ≤ C2 |r̄|γ as |r̄| → ∞, then ∇e is globally
bounded; we choose ψ(s) = arctan(sp ), where p > 0 is
sufficiently large so that p − 1 + γ ≤ 2βp.
V. E XPERIMENTAL VALIDATION
In this section the experiments with real robots are reported.
In one of these experiments, the desired path is an ellipse,
and the other experiments deals with the more sophisticated
Cassini oval (see Section III).
We test the results using the E-puck mobile robotic platform [43]. The experimental setup consists of the differential
wheeled E-puck robot, a server computer, an overhead camera
and a communication module. The robot is identified by a
data matrix marker on its top (see Fig. 5). The position
of the central point and the orientation of the marker are
recognized by a vision algorithm running at a control computer
connected to the overhead camera. The available workspace
is a planar area of 2.6x2 meters covered by the image of
1280x720 Px (henceforth the acronym Px is used for pixels).
For convenience we operate with pixels as coordinates. The
PC runs real-time calculation of the control actions based on
the pose information of the robot and computes the control
inputs for the robot. The results of computation as the desired
angular and linear velocities of the robot are translated into
the commands for its left and right wheels. The commands
from the control computer to the robot are sent via Bluetooth
at the fixed frequency of 20 Hz.
We verify the guidance algorithm for two closed trajectories,
passed by the robot in the clockwise direction. In both cases
the forward velocity was set to ur = 50 Px/s, and we set the
parameters of the algorithm kn = 3 and kδ = 2. We choose
ψ(s) = s, so that e(x, y) = ϕ(x, y). Below we describe the
trajectories and show the numerical data from the experiments.
A. Circulation along the ellipse
For our first experiment, we choose the elliptic trajectory,
defined by the function


(y − y0 )2
(x − x0 )2
2
.
(37)
+
−
R
ϕ(x, y) = ks
p2
q2

10

Obviously, the level curves of ϕ are ellipses. The path P
corresponds to the ellipse with semiaxes pR and qR, centered
at the point (x0 , y0 ).
For the experiment we choose x0 = 600Px, y0 = 350Px,
R = 400Px, and the semiaxes scale factors p = 1, q = 0.5 (see
Fig. 6). To provide ϕ(x, y) ∈ [−5; 5] in the working area, we
choose the scaling factor ks = 10−5 . The robots’s trajectories
(labeled a, b, c and d) under four different initial conditions
are shown in Fig. 6. The respective starting conditions are
a : (x = 472, y = 311, α = 0.0768),
b : (x = 30, y = 555, α = 0.0278),
c : (x = 408, y = 369, α = 2.1515),

(38)

d : (x = 78, y = 133, α = 4.0419)
(the coordinates x, y are in pixels, and α is in radians). Fig. 7
illustrates the dynamics of the tracking error ϕ(x(t), y(t)).

700

a
b
c
d
P

600

Y, Px

500

the ellipse and sufficiently far from it, then it enters into the
domain (39) with a “proper” orientation, satisfying (33), and
thus the problem of path following is solved. The practical
experiments, however, show that the robot converges to the
path from any point, except for the ellipse’ center, independent
of the initial robot’s orientation, whereas the mathematically
rigorous proof remains an open problem.
B. Circulation along the Cassini oval
For our second experiment the path chosen is referred to as
the Cassini oval (Fig. 8) and determined by the function
i
h

2
ϕ(x, y) = ks ∆x2 + ∆y 2 − 2q 2 ∆x2 − ∆y 2 − p4 + q 4 ,
∆x = (x − x0 ), ∆y = (y − y0 ).

We choose here ks = 10−10 , p = 330Px and q = 300Px. The
oval is centered at (x0 , y0 ) = (600, 350)Px.
Fig. 8 illustrates four trajectories (labeled a, b, c and d),
corresponding to the initial conditions
a : (x = 233, y = 184, α = 2.9287),
b : (x = 106, y = 202, α = 4.2487),

400

c : (x = 355, y = 343, α = 5.4071),
d : (x = 503, y = 619, α = 0.1022).

300

(40)

200

(the coordinates x, y are in pixels, and α is in radians). The
corresponding tracking errors are displayed in Fig. 9.

100

0
0

200

400

600

800

1000

1200

X, Px
700

a
b
c
d
P

Fig. 6. Elliptical path: the vector field and robot’s trajectories.
600

Y, Px

500

400

300

200

100

0
0

200

400

600

800

1000

1200

X, Px

Fig. 8. Cassini oval: the vector field and robot’s trajectories.

Fig. 7. Elliptical path: the dynamics of the tracking error e = ϕ(x, y).

In this example the set of critical points is C0 = {(x0 , y0 )},
which corresponds to e(x0 , y0 ) = −ks R2 and hence ec = R2 .
The set {(x, y) : |e(x, y)| < ec } consists of all points, where
0<

(y − y0 )2
(x − x0 )2
+
< 2R2 .
2
p
q2

(39)

Theorem 2 guarantees that starting at any of these points
and (33) holds, the robot converges to the desired path. Using
the geometry of the specific set (39), Remark 3 allows to prove
convergence for many other initial conditions: it is clear, for
instance, that if the robot starts from some point, lying outside

Fig. 9. Cassini oval: the dynamics of the tracking error e = ϕ(x, y).

11

As discussed in Section III, the set of critical points is C0 =
{(x0 ± q, y0 ), (x0 , y0 )}. If q 2 < p2 < 2q 2 , such inequalities
being valid in our example, then it can be checked that ec =
|e(x0 , y0 )| = p4 − q 4 , and the set {(x, y) : |e(x, y)| < ec } is a
domain, shown in Fig. 10. Theorem 2 guarantees that starting

800

700

600

The second path following algorithm to be compared with
the GVF method is the “nonlinear guidance law” (NGL) [3],
[44]. At the current robot position r̄, draw a circle of radius
R. The circle intersects the path P at two points q, q ′ . It is
supposed that the robot moves sufficiently close to the desired
path and the direction of the path following is fixed, so it is
always possible to choose the intersection point, lying ahead
of the robot on the path P (in Fig. 11, this is point q). The
absolute bearing αR to this point is the desired direction of
motion. The angular velocity controller is designed as

500

ωR = −kR (α − αR ) ,

400

300

200

100

0

-100
0

200

400

600

800

1000

1200

The constants R > 0, kR > 0 are the algorithm’s parameters.
For the sake of comparison, we will use the same elliptic
path, determined by the function (37) and the same parameters
of the GVF algorithm as in Subsect. V-A. The initial position
and heading of the robot are
x = 200; y = 450; α = 0.0278.

Fig. 10. Cassini oval: the set {(x, y) : |e(x, y)| < ec }.

at any of these points with the orientation satisfying (33), the
robot converges to the desired path. Similar to the case of the
elliptic path, Remark 3 allows to prove convergence for many
other initial conditions; experiments show that in fact the robot
converges to the desired path from any non-critical point.
VI. D ISCUSSION
In this section, we compare the proposed method with the
existing path following control algorithms. Potential extensions of the proposed method are also discussed.
A. The GVF method vs. other path following algorithms
It should be noticed that many approaches to path following, discussed in Introduction, are not applicable to general
smooth paths of non-constant curvature and nonholonomic
robots; for example, various virtual target approaches [15]
typically require to control the robot’s velocity. We compare
the GVF method with two commonly used algorithms [3]: the
projection-based line-of-sight (LOS) method [5], [16] and the
“nonlinear guidance law” (NGL) [3], [44].
The LOS approach assumes the existence of the unique
projection of the robot location r̄ onto the path P, that is, the
point P closest to r̄ in the Euclidean metric. The distance to
the path d = d(r̄, P ) serves as the tracking (or “cross-track”)
error. The desired robot’s orientation αLOS is the absolute
bearing to the unique point P∆ (Fig. 11), such that the vector
−−−→
P P∆ has length ∆ > 0 and is tangent to P at the projection
point. The direction of this vector is the desired direction of
path following (in our example, the robot circulates along the
path clockwise). The controller (23) is replaced by

It is clear that any path following controller is very sensitive
to the parameter choice, and varying the parameters one can
always make the convergence rate faster or slower. To compare
the behavior of different algorithms, one thus needs to choose
parameters, providing approximately the same convergence
rate; in our situation, the dynamics of the Euclidean distance
d(t) = dist(r̄(t), P) is similar to the GVF method, choosing
kLOS = 2 and ∆ = 70 in the LOS method and kR = 2
and R = 40 in the NGL algorithm (Fig. 12). As one can
see, the GVF method gives much better transient behavior
than the other algorithms, which have visible overshoots. To
eliminate the overshoot of the LOS method, one has to increase
the lookahead distance, but in this case the convergence to
the desired path becomes slower. The NGL method has the
largest overshoot; using this method, it is also impossible to
eliminate completely the tracking error since the controller
has no information about the changing curvature of the path.
To decrease the overshoot, one has to increase the radius R;
this however amplifies the oscillations in the tracking error d.
Thus the GVF method demonstrates better performance than
the considered alternative methods.

ωLOS = c(P )ur − kLOS (α − αLOS ) ,
where kLOS > 0 is a constant gain and c(P ) is the curvature
of the curve P at the projection point P . The parameters of the
algorithm are the lookahead distance ∆ and the gain kLOS .

Fig. 11. Geometric objects, employed by the LOS (red) and the NGL (blue)
path following algorithms

12

40

Guiding Vector Field
LOS Path Following
Nonlinear Guidance Law

35
30
25
20
15
10
5
0
-5
0

5

10

15

20

25

Fig. 12. The dynamic of the euclidean distance d

B. The previous works on GVF algorithms
The vector-field path following algorithms have been mainly
considered for straight lines and circular paths [3], [31]. The
GVF designed for these paths is the same as constructed in
Section III; the only difference is the path following controller:
unlike [31], in this paper we do not consider sliding-mode
algorithms and confine ourselves to simple linear controllers.
However, a general approach to the Lyapunov-based GVF design has been suggested in [30]. The desired path is supposed
to be the set of global minima of some predefined Lyapunov
function; e.g. for a path (2), our usual Lyapunov function (11)
can be chosen. Departing from this Lyapunov function, a
broad class of GVF has been suggested, including (14) as
a special case. The paper [30] also offers a methodology of
proving the convergence of integral curves to the desired path
(considered in our Section III) and design of path following
controllers. For the case of unicycle-like robot the control
algorithm, suggested in [30, Eq. (35)] is equivalent to our
controller (24). Describing a general framework, the paper [30]
focuses however on some specific problems of flight control,
namely, following a loiter circle (planar circle at the fixed
altitude) or some other planar path, which can be obtained
from a loiter circle by a smooth “warping” transformation,
preserving the vector field’s properties. Many assumptions
adopted in [30] are restrictive, e.g. the Lyapunov function
should be radially unbounded, excluding straight lines and
other unbounded paths from the consideration. The analysis
of the closed-loop system in [30] is not fully rigorous as it
completely ignores the problem of the solution existence up
to ∞, which becomes non-trivial when C0 6= ∅.
Unlike the previous papers [3], [30], [31], we consider a
general smooth path (2), which can be closed or unbounded;
for closed trajectories, we do not restrict the Lyapunov function to be radially unbounded. We provide a mathematically
rigorous analysis of our algorithm and specify an explicit set
of initial conditions, starting at which the robot evades the
critical points and converges to the desired trajectory.
C. Further extensions of our GVF method
The results presented in the paper can be extended in several
directions. First, the idea of GVF steering can be applied
to more general models than the simplest unicycle robot,
considered in this paper.
Second, the algorithm can be extended to the threedimensional case. Some special results in this direction have

Fig. 13. The desired path as the intersection of two surfaces

been reported in [24], [30]. The desired curvilinear path
P in the three dimensional space can be describing as an
intersection of two surfaces (Fig. 13)
∆

P = {(x, y, z) : ϕ1 (x, y, z) = 0 ∧ ϕ2 (x, y, z) = 0} ⊂ R3 ,
where ϕi ∈ C 2 (R3 → R), i = 1, 2. The main design idea
is that if the robot simultaneously approaches both surfaces,
eventually it will be brought to the prescribed path P. Thus,
we may define the guiding vector field by
v̄ = τ̄ − kn1 e1 n̄1 − kn2 e2 n̄2 ,

τ̄ = n̄1 × n̄2

where n̄i = ∇ϕi and ei = ϕi (x, y, z) are the “tracking errors”
(i = 1, 2). Subsequent design of the controller, steering the
robot to the GVF, is similar to the planar case (two angles are
controlled instead of a single angle). The technical analysis
of this algorithm requires however some extra assumptions of
non-degeneracy and is beyond the scope of this paper; the
convergence rate of the path following algorithm appears to
be very sensitive to the choice of ϕ1 and ϕ2 .
Third, the results can be extended to time-varying vector
fields [30], allowing to work with some types of transformations (such as translation, rotation and scaling) of the
predefined trajectory and enabling to use the GVF approach
for stand-off tracking of slowly moving targets [45].
Fourth, it is important to consider external disturbances,
unavoidable in practice and leading, in particular, to the lateral
drift of the robot. Some preliminary results on using draft
compensators in parallel with the path following controller
are reported in the conference papers [34], [46].
VII. C ONCLUSIONS
In this paper, we offer a new algorithm for path following
control of nonholonomic robots exploiting the idea of a
guiding vector field. Unlike the existing results, the desired
path can be an arbitrary smooth curve in its implicit form, i.e.
the zero set of a smooth function. We examine mathematical
properties of the algorithm and give global conditions for following asymptotically the desired path. The results have been
experimentally validated, using the E-puck wheeled robots.

13

R EFERENCES
[1] T. Kurfess, Ed., Robotics and Automation Handbook. CRC Press, 2005.
[2] B. Siciliano and O. Khatib, Eds., Springer Handbook of Robotics.
Springer, 2008.
[3] P. Sujit, S. Saripalli, and J. Borges Sousa, “Unmanned aerial vehicle path
following: A survey and analysis of algorithms for fixed-wing unmanned
aerial vehicless,” IEEE Control Systems, vol. 34, no. 1, pp. 42–59, 2014.
[4] E. Peymania and T. I. Fossen, “Path following of underwater robots using
Lagrange multipliers,” Robotics and Autonomous Systems, vol. 67, pp.
44–52, 2015.
[5] T. I. Fossen and K. Y. Pettersen, “On uniform semiglobal exponential
stability (USGES) of proportional line-of-sight guidance laws,” Automatica, vol. 50, no. 11, pp. 2912 – 2917, 2014.
[6] C. Samson, “Path following and time-varying feedback stabilization of
a wheeled mobile robot,” in Proc. of the ICARCV’92, Singapore, 1992.
[7] J. Ackermann, J. Guldner, W. Sienel, R. Steinhauser, and V. I. Utkin,
“Linear and nonlinear controller design for robust automatic steering,”
IEEE Trans. Control Syst. Tech., vol. 3, no. 1, pp. 132–143, 1995.
[8] A. Isidori, Nonlinear Control Systems, 3rd ed., M. Thoma, E. D. Sontag,
B. W. Dickinson, A. Fettweis, J. L. Massey, and J. W. Modestino, Eds.
Secaucus, NJ, USA: Springer-Verlag New York, Inc., 1995.
[9] L. Marconi and A. Isidori, “Mixed internal model-based and feedforward
control for robust tracking in nonlinear systems,” Automatica, vol. 36,
no. 7, pp. 993 – 1000, 2000.
[10] M. M. Seron, J. H. Braslavsky, P. V. Kokotovic, and D. Q. Mayne,
“Feedback limitations in nonlinear systems: From Bode integrals to
cheap control,” IEEE Trans. Autom. Control, vol. 44, no. 4, pp. 829–833,
1999.
[11] A. Aguiar, J. Hespanha, and P. Kokotovic, “Path-following for nonminimum phase systems removes performance limitations,” IEEE Trans.
Autom. Control, vol. 50, no. 2, pp. 234–239, 2005.
[12] A. Micaelli, C. Samson et al., “Trajectory tracking for unicycle-type and
two-steering-wheels mobile robots,” Tech. Rep. No. 2097, INRIA, 1993.
[13] C. Samson, “Control of chained systems application to path following
and time-varying point-stabilization of mobile robots,” IEEE Trans.
Autom. Control, vol. 40, no. 1, pp. 64–77, 1995.
[14] A. S. Matveev, M. Hoy, J. Katupitiya, and A. V. Savkin, “Nonlinear
sliding mode control of an unmanned agricultural tractor in the presence
of sliding and control saturation,” Robotics and Autonomous Systems,
vol. 61, no. 9, pp. 973 – 987, 2013.
[15] D. Soetanto, L. Lapierre, and A. Pascoal, “Adaptive, non-singular pathfollowing control of dynamic wheeled robots,” in Proc. of 42nd IEEE
Conf. on Decision and Control, vol. 2, Dec 2003, pp. 1765–1770 Vol.2.
[16] M. Breivik and T. Fossen, “Principles of guidance-based path following
in 2D and 3D,” in Proc. of 44th IEEE Conf. on Decision and Control
and European Control Conf. CDC-ECC’05., Dec 2005, pp. 627–634.
[17] A. Fradkov, I. Miroshnik, and V. Nikiforov, Nonlinear and Adaptive
Control of Complex Systems, ser. Mathematics and Its Applications.
Springer, 1999.
[18] M. I. El-Hawwary and M. Maggiore, “Case studies on passivity-based
stabilisation of closed sets,” International Journal of Control, vol. 84,
no. 2, pp. 336–350, 2011.
[19] A. Hladio, C. Nielsen, and D. Wang, “Path following for a class of
mechanical systems,” IEEE Trans. Control Syst. Tech., vol. 21, no. 6,
pp. 2380–2390, Nov 2013.
[20] A. Akhtar, C. Nielsen, and S. Waslander, “Path following using dynamic transverse feedback linearization for car-like robots,” IEEE Trans.
Robotics, vol. 31, no. 2, pp. 269–279, April 2015.
[21] S. Kolyubin and A. Shiriaev, “Dynamics-consistent motion planning for
underactuated ships using virtual holonomic constraints,” in Proc. of
IEEE OCEANS 2014, Sept 2014, pp. 1–7.
[22] Y. Kapitanyuk and S. Chepinsky, “Control of mobile robot following a
piecewise-smooth path,” Gyroscopy and Navigation, vol. 4, no. 4, pp.
198–203, 2013.
[23] Y. Kapitanyuk, S. Chepinskiy, and A. Kapitonov, “Geometric path
following control of a rigid body based on the stabilization of sets,”
in Proc. of the 19th IFAC World Congress, Cape Town, South Africa,
2014, pp. 7342–7347.
[24] J. Wang, I. A. Kapitaniuk, S. A. Chepinskiy, D. Liu, and A. J.
Krasnov, “Geometric path following control in a moving frame,” IFACPapersOnLine, vol. 48, no. 11, pp. 150 – 155, 2015.
[25] R. Gill, D. Kulić, and C. Nielsen, “Spline path following for redundant
mechanical systems,” IEEE Trans. Robotics, vol. 31, no. 6, pp. 1378–
1392, 2015.

[26] D. Panagou and V. Kumar, “Cooperative visibility maintenance for
leader-follower formations in obstacle environments,” IEEE Trans.
Robotics, vol. 30, no. 4, pp. 831–844, 2014.
[27] M. Hoy, A. Matveev, and A. V. Savkin, “Algorithms for collisionfree navigation of mobile robots in complex cluttered environments:
a survey,” Robotica, vol. 33, no. 3, pp. 463–497, 2015.
[28] A. S. Matveev, K. S. Ovchinnikov, and A. V. Savkin, “A method of
reactive 3d navigation for a tight surface scan by a nonholonomic mobile
robot,” Automatica, vol. 75, pp. 119–126, 2017.
[29] A. S. Matveev, M. C. Hoy, and A. V. Savkin, “Extremum seeking
navigation without derivative estimation of a mobile robot in a dynamic
environmental field,” IEEE Transactions on Control Systems Technology,
vol. 24, no. 3, pp. 1084–1091, 2016.
[30] D. A. Lawrence, E. W. Frew, and W. J. Pisano, “Lyapunov vector fields
for autonomous unmanned aircraft flight control,” Journal of Guidance,
Control, and Dynamics, vol. 31, no. 5, pp. 1220–1229, 2008.
[31] D. Nelson, D. Barber, T. McLain, and R. Beard, “Vector field path
following for miniature air vehicles,” IEEE Trans. Robotics, vol. 23,
no. 3, pp. 519–529, June 2007.
[32] L. Valbuena Reyes and H. Tanner, “Flocking, formation control, and
path following for a group of mobile robots,” IEEE Trans. Control Syst.
Tech., vol. 23, no. 4, pp. 1268–1282, July 2015.
[33] V. Goncalves, L. Pimenta, C. Maia, B. Dutra, and G. Pereira, “Vector
fields for robot navigation along time-varying curves in n-dimensions,”
IEEE Trans. Robotics, vol. 26, no. 4, pp. 647–659, Aug 2010.
[34] R. Beard, J. Ferrin, and J. Humpherys, “Fixed wing UAV path following
in wind with input constraints,” IEEE Trans. Control Syst. Tech., vol. 22,
no. 6, pp. 2103–2117, Nov 2014.
[35] T. I. Fossen, Guidance and control of ocean vehicles. John Wiley &
Sons, Ltd, 1994.
[36] R. W. Beard and T. W. McLain, Small unmanned aircraft: Theory and
practice. Princeton university press, 2012.
[37] A. Tsourdos, B. White, and M. Shanmugavel, Cooperative Path Planning of Unmanned Aerial Vehicles, ser. Aerospace Series. Wiley, 2010.
[38] T. I. Fossen, K. Y. Pettersen, and R. Galeazzi, “Line-of-sight path
following for Dubins paths with adaptive sideslip compensation of drift
forces,” IEEE Trans. Control Syst. Tech., vol. 23, no. 2, pp. 820–827,
2015.
[39] P. Morin and C. Samson, “Trajectory tracking for non-holonomic
vehicles: overview and case study,” in Proc. of 4th Int. Workshop on
Robot Motion and Control RoMoCo’04., June 2004, pp. 139–153.
[40] G. Oriolo, A. D. Luca, and M. Vendittelli, “WMR control via dynamic
feedback linearization: Design, implementation, and experimental validation,” IEEE Trans. Control Syst. Tech., vol. 10, no. 6, pp. 835–852,
2002.
[41] R. Potrie and P. Monzo’n, “Local implications of almost global stability,”
Dynam. Syst., vol. 24, no. 1, pp. 109–115, 2009.
[42] H. Khalil, Nonlinear systems. Englewood Cliffs, NJ: Prentice-Hall,
1996.
[43] F. Mondada, M. Bonani, X. Raemy, J. Pugh, C. Cianci, A. Klaptocz,
S. Magnenat, J.-C. Zufferey, D. Floreano, and A. Martinoli, “The epuck, a robot designed for education in engineering,” in Proc. of the 9th
Conf. on Auton. Robot Syst. and Competitions, 2009, pp. 59–65.
[44] S. Park, J. Deyst, and J. P. How, “Performance and lyapunov stability
of a nonlinear path following guidance method,” Journal of Guidance,
Control, and Dynamics, vol. 30, no. 6, pp. 1718–1728, 2007.
[45] H. Oh, S. Kim, H. Shin, and A. Tsourdos, “Coordinated standoff tracking
of moving target groups using multiple UAVs,” IEEE Transactions on
Aerospace and Electronic Systems, vol. 51, no. 2, pp. 1501–1514, April
2015.
[46] H. G. de Marina, Y. A. Kapitanyuk, M. Bronz, G. Hattenberger, and
M. Cao, “Guidance algorithm for smooth trajectory tracking of a fixed
wing UAV flying in wind flows,” accepted to IEEE Int. Conf. Robotics
and Automation (ICRA), 2017.

