Event-driven Trajectory Optimization for Data Harvesting in
Multi-Agent Systems
Yasaman Khazaeni and Christos G. Cassandras
Division of Systems Engineering
and Center for Information and Systems Engineering
Boston University, MA 02446

arXiv:1608.06336v1 [math.OC] 22 Aug 2016

yas@bu.edu,cgc@bu.edu

Abstract— We propose a new event-driven method for on-line
trajectory optimization to solve the data harvesting problem: in
a two-dimensional mission space, N mobile agents are tasked
with the collection of data generated at M stationary sources
and delivery to a base with the goal of minimizing expected
collection and delivery delays. We define a new performance
measure that addresses the event excitation problem in eventdriven controllers and formulate an optimal control problem.
The solution of this problem provides some insights on its
structure, but it is computationally intractable, especially in
the case where the data generating processes are stochastic.
We propose an agent trajectory parameterization in terms of
general function families which can be subsequently optimized
on line through the use of Infinitesimal Perturbation Analysis
(IPA). Properties of the solutions are identified, including
robustness with respect to the stochastic data generation process
and scalability in the size of the event set characterizing
the underlying hybrid dynamical system. Explicit results are
provided for the case of elliptical and Fourier series trajectories
and comparisons with a state-of-the-art graph-based algorithm
are given.

I. I NTRODUCTION
Systems consisting of cooperating mobile agents have
been extensively studied and used in a broad spectrum of
applications such as environmental sampling [1],[2], surveillance [3], coverage [4],[5],[6], persistent monitoring [7],[8],
task assignment [9], and data harvesting and information
collection [10],[11],[12].
The data harvesting problem in particular (and its variant,
the “minimum latency” problem [13]) arises in many settings where wireless sensor networks (WSNs) are deployed
for purposes of monitoring the environment, road traffic,
infrastructure for transportation and for energy distribution,
surveillance, and a variety of other specialized purposes [14],
[15]. Although many efforts focus on the analysis of the vast
amount of data gathered, we must first ensure the existence of
robust means to collect all data in a timely fashion when the
size of the network and the level of node interference do not
allow for a fully connected wireless system. In such cases,
sensors can locally gather and buffer data, while mobile
elements (e.g., vehicles, aerial drones) retrieve the data from
each part of the network. Similarly, mobile elements may
themselves be equipped with sensors and visit specific points
The authors’ work is supported in part by NSF under grants CNS1239021, ECCS-1509084, and IIP-1430145, by AFOSR under grant
FA9550-12-1-0113, by ONR under grant N00014-09-1-1051.

of interest, called “targets”, where a direct communication
path does not exist between them and the central sink or base
node where the data must be delivered. In a delay-tolerant
system, the control scheme commonly used is to deploy
mobile agents referred to as “data mules”, “message ferries”
or simply “ferries” [16],[17],[18],[19]. The mobile agents
visit the data generation nodes and collect data which are
then delivered to the base. Moreover, since agents generally
have limited buffer sizes, visits to the base are also needed
once a buffer is full; the same is true due to limited energy
that requires them to be periodically recharged. In general,
the paths followed by the mobile agents need to be optimized
(in some sense to be defined) so as to ensure timely delivery
of data through sufficiently frequent visits at each data
source and the base and within the constraints of a given
environment (e.g., an urban setting).
Interestingly, there are analogs to the data harvesting
problem outside the sensor network realm. For instance,
in disaster planning, evacuation and rescue operations,
pickup/delivery and transportation systems, and UAV surveillance operations, the general theme involves a network
of cooperating mobile agents that need to frequently visit
points of interest and transfer data/goods/people to a base.
Base visits may also be needed to recharge/renew power
supply/fuel. For example, the flying time span of a drone on
a single battery charge is limited, so that flight trajectories
need to be optimized and returns to base may be scheduled
for recharge or loading/unloading. Thus, we may view data
harvesting in the broader context of a multi-agent system
where mobile agents must cooperatively design trajectories
to visit a set of targets and a base so as to optimize one or
more performance ctiteria.
Having its root in wireless sensor networks, the data
harvesting problem is normally studied on a directed or
undirected graph where minimum length tours or sub-tours
are to be found. This graph topology view of the problem
utilizes a multitude of routing and scheduling algorithms
developed for wireless sensor networks (e.g., [20],[21],[22]
and references therein). One of its main advantages is the
ability to accommodate environment constraints (e.g., obstacles) by properly selecting graph edges; thus, movements
inside a building or within a road network are examples
where a graph topology is suitable. On the other hand, these
methods also have several drawbacks: they are generally

combinatorially complex; they do not account for limitations
in motion dynamics which should not, for instance, allow
an agent to form a trajectory consisting of a sequence of
straight lines; they become computationally infeasible as online methods in the presence of stochastic effects such as
random target rewards or failing agents, since the graph
topology has to be re-evaluated as new information becomes
available.
As an alternative to the graph-based approach, we view
data harvesting as a trajectory optimization problem in a
two-dimensional space, where the mobile agents are freely
(or with some specified constraints) moving and “visit”
targets whenever they reach their vicinity. For example, each
target can be assumed to have a finite range within which a
mobile agent can initiate wireless communication with it and
exchange data. These trajectories do not necessarily consist
of straight lines, i.e., edges in an underlying graph topology;
therefore, an advantage of this approach is that an agent trajectory can be designed to conform to physical limitations in
the agent’s mobility. In addition, trajectories can be adjusted
on line when target locations are uncertain. Constraining such
trajectories to obstacles is also still possible. Such continuous
topologies have been used in [11], where the problem is
viewed as a polling system with a mobile server visiting
data queues at fixed targets and trajectories are designed
to stabilize the system, keeping queue contents (modeled
as fluid queues) uniformly bounded; and in [23] where
parameterized trajectories are optimized to solve a multiagent persistent monitoring problem.
A key benefit of a trajectory optimization view of the
data harvesting problem is the ability to parameterize the
trajectories with different types of functional representations
and then optimize them over the given parameter space.
This reduces a dynamic optimization problem into a much
simpler parametric optimization one. If the parametric trajectory family is broad enough, we can recover the true
optimal trajectories; otherwise, we can approximate them
within a desired accuracy. Moreover, adopting a parametric
family of trajectories has several additional benefits. First, it
allows trajectories to be periodic, often a desirable property
in practice. Second, it allows one to restrict solutions to
trajectories with desired features that the true optimal cannot
have, e.g., smoothness properties required for physically
realizable agent motion.
In this paper, we cast data harvesting as an optimal control
problem. Defining an appropriate optimization criterion is
nontrivial in this problem (as we will explain) and introducing appropriate performance metrics is the first contribution
of this work. Obtaining optimal agent trajectories ultimately
requires the solution of a two point boundary value problem
(TPBVP). Although a complete solution of such a TPBVP is
computationally infeasible in general, we identify structural
properties of the optimal control policy which allow us to
reduce the agent-target/base interaction process to a hybrid
system with a well-defined set of events that cause discrete
state transitions. The second contribution is to formulate
and solve an optimal parametric agent trajectory problem.

In particular, similar to the idea introduced in [23], we
represent an agent trajectory in terms of general function
families characterized by a set of parameters that we seek to
optimize, given an objective function. We consider elliptical
trajectories as well as the much richer set of Fourier series
trajectory representations. We then show that we can make
use of Infinitesimal Perturbation Analysis (IPA) for hybrid
systems [24] to determine gradients of the objective function
with respect to these parameters and subsequently obtain
(at least locally) optimal trajectories. This approach also
allows us to exploit (i) robustness properties of IPA to
allow stochastic data generation processes, (ii) the eventdriven nature of the IPA gradient estimation process which
is scalable in the event set of the underlying hybrid dynamic
system, and (iii) the on-line computation which implies that
trajectories adjust as operating conditions change (e.g., new
targets); in contrast, the solution of a TPBVP is computationally challenging even for strictly off line methods. Finally, we
provide comparisons of our approach to algorithms based on
a graph topology of the mission space. These comparisons
show that while the latter generate a spatial partitioning
of the target set among agents, our approach results in a
temporal partitioning which adds robustness with respect to
agent failures or other environmental changes. The graphbased approaches are mostly offline and normally assume
the agents would be able to travel straight lines and meet
targets in exact locations. On the other hand the trajectory
optimization approach, allows us to accommodate limitations
in agent mobility and to adjust trajectories on line.
In Section II we formulate the data harvesting problem using a queueing model and present the underlying
hybrid system. In Section III we provide a Hamiltonian
analysis leading to a TPBVP. In Section III we formulate
the alternative problem of determining optimal trajectories
based on general function representations and provide solutions through a gradient-based algorithm using IPA for two
particular function families. Section IV presents numerical
results and comparisons with state of the art data harvesting
algorithms and Section VI contains conclusions.
II. P ROBLEM F ORMULATION
We consider a data harvesting problem where N mobile
agents collect data from M stationary targets in a twodimensional mission space S. Each agent may visit one or
more of the M targets, collect data from them, and deliver
them to a base. It then continues visiting targets, possibly the
same as before or new ones, and repeats this process. The
objective of the agent team is to minimize data collection and
delivery delays over all targets within a fixed time interval
T . This minimization problem is formalized in the sequel.
The data harvesting problem described above can be
viewed as a polling system where mobile agents are serving
the targets by collecting data and delivering it to the base.
As seen in Fig. 1, there are three sets of queues. The
first set includes the data contents Xi (t) ∈ R+ at each
target i = 1, ..., M where we use σi (t) as the instantaneous
inflow rate. In general, we treat {σi (t)} as a random process

X1
p11

.X. .i
...

pB1

...

Y1

...
Yi

XM
pM N
Zij

Fig. 1.

pBN
YM

Data harvesting queueing model for M targets and N agents

assumed only to be piecewise continuous; we will treat it as
a deterministic constant only for the Hamiltonian analysis in
the next section. Thus, at time t, Xi (t) is a random variable
resulting from the random process {σi (t)}.
The second set of queues consists of data contents
Zij (t) ∈ R+ onboard agent j collected from targets i =
1, ..., M . The last set consists of queues Yi (t) ∈ R+ containing data at the base, one queue for each target, delivered
by some agent j. Note that {Zij (t)} and {Yi (t)} are also
random processes.
In Fig. 1 collection and delivery switches are shown by
pij and pBj (formally defined in the sequel). These switches
are “on” when agent j is connected to target i or the
base respectively. All queues are modeled as flow systems
whose dynamics are given next (however, as we will see, the
agent trajectory optimization is driven by events observed in
the underlying system where queues contain discrete data
packets so that this modeling device has minimal effect on
our analysis).
Let sj (t) = [sxj (t), syj (t)] ∈ S be the position of agent j
at time t, Then, the state of the system can be defined as
X(t) = [X1 (t), . . . , XM (t), Y1 (t), . . . , YM (t),
(1)
Z11 (t), . . . , ZM N (t), sx1 (t), sy1 (t), . . . , sxN (t), syN (t)]
The position of the agent follows single integrator dynamics
at all times:
ṡxj (t) = uj (t) cos θj (t),
ṡyj (t) = uj (t) sin θj (t) (2)
y
x
sj (0) = XB
sj (0) = YB , ∀j
where uj (t) is the scalar speed of the agent (normalized so
that 0 ≤ uj (t) ≤ 1), 0 ≤ θj (t) < 2π is the angle relative
to the positive direction and [XB , YB ] is the location of the
base. Thus, we assume that the agent controls its orientation
and speed. The agent states {sj (t)}, j = 1, . . . , N , are also
random processes since the controls are generally dependent
on the random queue states. Thus, we ensure that all random
processes are defined on a common probability space.
An agent is represented as a particle, so that we will
omit the need for any collision avoidance control. The agent
dynamics above could be more complicated without affecting
the essence of our analysis, but we will limit ourselves here
to (2).
We consider a set of data sources as points wi ∈ S, i =
1, . . . , M, with associated ranges rij , so that agent j can
collect data from wi only if the Euclidean distance dij (t) =
kwi − sj (t)k satisfies dij (t) ≤ rij . Similarly, the base is at
wB = [XB , YB ] ∈ S which receives all data collected by
the agents and an agent can only deliver data to the base
if the Euclidean distance dBj (t) = kwBj − sj (t)k satisfies

dBj (t) ≤ rBj . Using p : S × S → [0, 1], we define a function
pij (t) representing the collection switches in Fig. 1 as:
pij (t) = p(wi , sj (t))
(3)
pij (t) is viewed as the normalized data collection rate from
target i when the agent is at sj (t) and we assume that:
(A1) it is monotonically non-increasing in the value of
dij (t) = kwi − sj (t)k, and (A2) it satisfies pij (t) = 0 if
dij (t) > rij . Thus, pij (t) can model communication power
constraints which depend on the distance between a data
source and an agent equipped with a receiver (similar to
the model used in [11]) or sensing range constraints if an
agent collects data using on-board sensors. For simplicity, we
will also assume that: (A3) pij (t) is continuous in dij (t).
Similarly, we define:
pBj (t) = p(wB , sj (t))
(4)
The maximum rate of data collection from target i by agent
j is µij , so that the instantaneous rate is µij pij (t) if j is
connected to i. We will assume that: (A4) only one agent at
a time is connected to a target i even if there are other agents
l with pil (t) > 0; this is not the only possible model, but we
adopt it based on the premise that simultaneous downloading
of packets from a common source creates problems of proper
data reconstruction at the base.
We can now define the dynamics of the queue-related
components of the state vector in (1). The dynamics of Xi (t),
assuming 
that agent j is connected to it, are
Ẋi (t) =

0
if Xi (t) = 0 and σi (t) ≤ µij pij (t)
σi (t) − µij pij (t)
otherwise

(5)

Obviously, Ẋi (t) = σi (t) if pij (t) = 0, j = 1, . . . , N .
In order
the dynamics
of Zij (t), let
( to express


µ̃ij (t) =

min
µij

σi (t)
pij (t) , µij

if Xi (t) = 0 and pij (t) > 0
otherwise

(6)
This gives
 us the dynamics:

0
if Zij (t) = 0 and µ̃ij (t)pij (t) − βij pBj (t) ≤ 0
µ̃ij (t)pij (t) − βij pBj (t)
otherwise

(7)
where βij is the maximum rate of data from target i delivered
to B by agent j. For simplicity, we assume that: (A5)
kwi − wB k > rij + rBj for all i = 1, . . . , M and j =
1, . . . , N , i.e., the agent cannot collect and deliver data at
the same time. Therefore, in (7) it is always the case that
for all i and j, pij (t)pBj (t) = 0. Finally, the dynamics of
Yi (t) depend on Zij (t), the content of the on-board queue
of each agent j from target i as long as pBj (t) > 0. We
PN
define βi (t) =
j=1 βij pBj (t)1[Zij (t) > 0] as the total
instantaneous delivery rate for target i data, so that the
dynamics of Yi (t) are:
Ẏi (t) = βi (t)
(8)
Żij (t) =

Hybrid System model: Taking into account the state vector
in (1) and the dynamics in (2), (5), (7) and (8), the data
harvesting process is a stochastic hybrid system. Discrete
modes of the system are defined by intervals over which (i)
agents are visiting a target, (ii) agents are visiting the base,
and (iii) agents are moving when not connected to any target
or base. The events that trigger mode transitions are defined
in Table I (the superscript 0 denotes events causing a variable
to reach a value of zero from above and the superscript +

where the weight αi represents the relative importance factor
of target i. Similarly, we define a weighted sum of expected
base queue content:
M
X

J2 (t) = E
αi Yi (t)
(11)
i=1

Fig. 2.

One target i and one agent j hybrid automaton

denotes events causing a variable to become strictly positive
from a zero value). We also use the following definitions:
+
d+
ij (t) = max(0, dij (t) − rij ), dBj (t) = max(0, dBj (t) − rBj ) (9)
The variables above are zero if agent j is within range of
target i or the base respectively.
Event Name
1. ξi0
2. ξi+
0
3. ζij
+
4. δij
0
5. δij
6. ∆+
j
7. ∆0j

TABLE I
H YBRID S YSTEM E VENTS
Description
Xi (t) hits 0, for i = 1, . . . , M
Xi (t) leaves 0, for i = 1, . . . , M .
Zij (t) hits 0, for i = 1, . . . , M , j = 1, . . . , N
d+
ij (t) leaves 0, for i = 1, . . . , M , j = 1, . . . , N
d+
ij (t) hits 0, for i = 1, . . . , M , j = 1, . . . , N
+ (t) leaves 0, for j = 1, . . . , N
dBj
+ (t) hits 0, for j = 1, . . . , N
dBj

Observe that each of the events in Table I causes a change
in at least one of the state variables in (5), (7), (8). For
example, ξi0 (i.e., the queue at target i is emptied) causes a
switch in (5) from Ẋi (t) = σi (t) − µij pij (t) to Ẋi (t) = 0.
+
Also note that we have omitted an event ζij
for Zij (t)
becoming strictly positive since this event is immediately
0
induced by δij
when agent j comes within range of target i
and starts collecting data causing Zij (t) > 0 if Zij (t) = 0
and Xi (t) > 0. Finally, note that all events are directly
observable during the execution of any agent trajectory and
they do not depend on our flow-based queueing model.
For example, if Xi (t) becomes zero, this defines event
ξi0 regardless of whether the corresponding queue is based
on a flow or on discrete data packets; this observation is
very useful in the sequel. A high-level hybrid automaton is
presented in Fig. 2 for a single target i and one agent j
system. This automaton becomes much more complicated
once more targets and agents are included.
A. Performance Measures
Our objective is to maintain minimal data content at all
target queues while also maximizing the contents of the
delivered data at the base queues. Thus, we define J1 (t) to be
the weighted sum of expected target queue content (recalling
that {σi (t)} are random processes):
M
X

J1 (t) = E
αi Xi (t)
(10)
i=1

Therefore, a tentative optimization objective is the convex
combination of (10) and (11) leading to the minimization
problem:
Z

1 T
min J(T ) =
qJ1 (t) − (1 − q)J2 (t) dt (12)
T 0
u(t),θ(t)
where u and θ are the vectors formed by the agent speed
and headings and q ∈ [0, 1] is a weight capturing the relative
importance of collected data as opposed to delivered data.
This performance measure captures the collection and
delivery of data which are processes taking place while an
agent is connected to any of the targets or the base. However,
it lacks any information regarding the interaction of an agent
with the environment when this agent is not connected to
any target or base and is due to the fact that the environment
has only a finite number of points of interest (targets). This
motivates two new performance measures we introduce next.
Agent Utilization: In accessing the targets, we must ensure
that the agents maximize their utilization, i.e., the fraction of
time spent performing a useful task by being within range
of a target or the base. Equivalently, we aim to minimize
the non-productive idling time of each agent during which
it is not visiting any target or the base. Using (9), agent j is
+
idling when d+
ij (t) > 0 for all i and dBj (t) > 0. We define
the idling function Ij (t) as follows:
!
M
Y
+
+
Ij (t) = log 1 + dBj (t)
dij (t)
(13)
i=1

This function has the following properties. First, Ij (t) = 0
if and only if the product term inside the bracket is zero,
i.e., agent j is visiting a target or the base; otherwise,
Ij (t) > 0. Second, Ij (t) is monotonically nondecreasing in
the number of targets M . The logarithmic function is selected
to prevent the value of Ij (t) from dominating those of J1 (·)
and J2 (·) when included in a single objective function. Thus,
we define:
N
X

J3 (t) = E
Ij (t)
(14)
j=1

Note that Ij (t) is also a random variable since it is a function
of the agent states sj (t), j = 1, . . . , N .
Event Excitation: As mentioned in the Introduction, our
goal is to develop an event-driven approach for on-line
trajectory optimization. In other words, we seek a controller
whose actions are based on events observed during the
operation of the hybrid system described earlier. Clearly,
the premise of this approach is that the events involved
are observable so as to “excite” the underlying event-driven
controller. However, it is not always obvious that these events
actually take place under every feasible control, in which
case the controller may be useless. This is illustrated in Fig.
3 where two different trajectories are shown for the agent.
The blue and red trajectories pass through none of the targets.
Consequently, there is an infinite number of trajectories for

Fig. 3.

Two trajectories with same objective function value

which the value of the objective
Z function in (12) is given by
q T
J(T ) =
tσi (t)dt
(15)
T 0
which is simply the total amount of data generated at all
targets through (5). This cannot be affected by any eventdriven control action, since none of the events in Table I is
excited. Clearly, the same is true for J3 (t) in (14).
To address this issue, our goal is to “spread” each target
cost (accumulated data) Xi (t) over all w ∈ S. This will create a potential field throughout the mission space. Following
[25], we begin by determining the convex hull produced by
the targets, since the trajectories need not go outside this
polygon. Let T = {w1 , w2 , · · · , wM } be the set of all target
points. Then, their convex hull is

X
M
X
βi w i |
βi = 1, ∀i, βi ≥ 0
(16)
C=
i=1

i

Given that C ⊂ S, we seek a function R(w, t) that satisfies
the following property for some constants ci > 0:
Z
M
X
R(w, t)dw =
ci Xi (t)
(17)
C

i=1

Thus, R(w, t) can be viewed as a time-varying density
function defined for all points w ∈ C which generates a total
cost equivalent to a weighted sum of the target data Xi (t),
i = 1, . . . , M . Letting d+
i (w) = max(kw − wi k, ri ), where
ri = minj rij , we then define:
M
X
αi Xi (t)
R(w, t) =
(18)
d+
i (w)
i=1
Intuitively, a target’s cost (numerator above) is spread over
all w ∈ S so as to obtain the “total weighted cost density”
at w. Note that d+
i (w) is defined to ensure that the target
cost remains positive and fixed for all points w ∈ C(wi ). In
order to illustrate this construction, Fig. 4(a) shows a sample
mission space with 9 target locations and Fig 4(b) shows the
value of R(w, t) at a specific time t.
Proposition 1: There exist ci > 0, i = 1, . . . , M , such
that:
Z
M
X
R(w, t)dw =
ci xi (t)
(19)
C

i=1

Proof: See Appendix I.
Using the same idea for theP
base, we define
M
αi Zij
RBj (w, t) = i=1
(20)
dB+ (w)
+
where dB (w) = max(kwB − wk, rB ) is a constant and rB =
minj rBj .
Proposition 1 asserts that the total cost due to data accumulated at targets may indeed be spread over all points

(a) Mission Space with dots as tar- (b) R Function at a sample time t
get locations
Fig. 4. R function illustration

in the mission space allowing an agent to “interact” with
these points through the resulting potential field. In order to
capture this interaction (see also [25]), we define the travel
cost for an agent j to reach point w as the quadratic of the
distance between them ksj (t) − wk2 and the total travel cost
as
N
X
ksj (t) − wk2
(21)
P (w, s(t)) =
j=1

Using these definitions we can now introduce a new performance metric:


 PN R

J4 (t) = E
j=1 S R(w, t) + RBj (w, t) Pj (w, t)dw (22)
Terminal Cost: Since we address the data harvesting
problem over a finite interval T , we define a terminal cost
at T capturing the expected value of the amount of data left
on board the agents:
M N

1 XX
αi Zij (T )
(23)
Jf (T ) = E
T
i=1 j=1
Clearly, the effect of this term vanishes as T → ∞ as long as
all E[Zij (T )] remains bounded. Moreover, if we constrain
trajectories to be periodic, this terminal cost may be omitted.
Finally, for simplicity, we will assume that αi = 1 for all i.
B. Optimization Problem
We can now formulate a stochastic optimization problem
P1 where the control variables are the agent speeds and
headings denoted by the vectors u(t) = [u1 (t), . . . , uN (t)]
and θ(t) = [θ1 (t), . . . , θN (t)] respectively (omitting their
dependence on the full system state at t). Combining the
components in (14),Z (22) and (23) we obtain:
1 T q
(1 − q)
min J(T ) =
J1 (t) −
J2 (t)
(24)
T 0
MX
MY
u(t),θ(t)

1
1
1
+
J3 (t) +
J4 (t) dt +
Jf (T )
MI
MR
MZ
where we introduce the normalizing factors MX , MY , MI ,
MR and MZ . This normalization ensures that all different
components of J(T ) are in the same range so that none of
them may dominate any other. We use an upper bound for
the value of each component as follows, where we assume
that σi (0) > 0 w.p. 1:
X
MX = M Y = MZ = T
σi (0)
(25)
i



MI = log 1 +

q

L21 + L22

M +1 

(26)

P
ri
T L1 L2 (L21 + L22 ) X
MR =
σi (0), r = i
(27)
r
M
i
where L1 and L2 define the size of the rectangular mission
space (the normalization factors can easily be adapted to
different mission space shapes). Observe that an unattainable
lower bound of the total objective function is −(1−q) which
occurs if J1 = J3 = J4 = 0 and J2 is at its maximum of 1.
If q = 0, then the lower bound is at its minimum of -1.
III. O PTIMIZATION M ETHODOLOGY
In this section, we consider problem P1 in a setting
where all data arrival processes are deterministic, so that all
expectations in (10)-(23) degenerate to their arguments. We
proceed with a standard Hamiltonian analysis leading to a
Two Point Boundary Value Problem (TPBVP) [26] where
the states and costates are known at t = 0 and t = T
respectively. We define the costate vector associated to (1):
λ(t) = [λ1 (t), . . . , λM (t), γ1 (t), . . . , γM (t),
y
x
φ11 (t), . . . , φM N (t), η1x (t), η1y (t), . . . , ηN
(t), ηN
(t)] (28)
The Hamiltonian is h
i
1
qJ1 (t) − (1 − q)J2 (t) + J3 (t) + J4 (t)
H(X, λ, u, θ) =
T X
X
XX
+
λi (t)Ẋi (t) +
γi (t)Ẏi (t) +
φij (t)Żij (t)
i

+

i

X

i

ηjx (t)uj (t) cos θj (t)

+

j


ηjy (t)uj (t) sin θj (t)

(29)

j

where the costate equations are
λ̇i (t) = −

Z
αi Pj (w, t) 
∂H
1 q
1 X
dw λi (T ) = 0
=−
+
∂Xi
T MX
MR j S d+
i (w)

∂H
1−q
=
γi (T ) = 0
∂Yi
T MY
Z
αi Pj (w, t)
∂H
1
φ̇ij (t) = −
=−
dw
∂Zij
MR S dB+ (w)

γ̇i (t) = −

φij (T ) =

∂Jf
∂Zij

T

∂H
∂sxj
"
Z
 ∂Pj (w, t)
1 ∂Ij (t)
1 X
=−
+
R(w, t) + RBj (w, t)
dw
x
T MI ∂sj
T MR j S
∂sxj
#
X ∂
X ∂
X ∂
+
λ (t)Ẋi (t) +
γ (t)Ẏi (t) +
φij (t)Żij (t)
x i
x i
∂s
∂s
∂sxj
j
j
i
i
i
∂H
η̇jy (t) = − y
∂sj
"
Z
 ∂Pj (w, t)
1 X
1 ∂Ij (t)
=−
R(w, t) + RBj (w, t)
dw
y +
T MI ∂sj
T MR j S
∂syj
#
X ∂
X ∂
X ∂
+
φij (t)Żij (t)
y λi (t)Ẋi (t) +
y γi (t)Ẏi (t) +
∂s
∂s
∂syj
j
j
i
i
i

η̇jx (t) = −

ηjx (T ) = ηjy (T ) = 0

From (29), after some
trigonometric manipulations, we
i get
1h
H(X, λ, u, θ) =
qJ1 (t) − (1 − q)J2 (t) + J3 (t) + J4 (t)
T X
X
XX
+
λi (t)Ẋi (t) +
γi (t)Ẏi (t) +
φij (t)Żij (t)
i

+

X
j

i

q

uj (t)sgnηjy (t)

i
2

ηjx (t) +

2
ηjy (t)

η x (t)

j

(30)

sin(θj (t) + ψj (t))

where tan ψj (t) = ηjy (t) for ηjy (t) 6= 0 and ψj (t) =
j
sgnηjx (t) π2 if ηjy (t) = 0. Applying the Pontryagin principle
to (29) with (u∗ , θ ∗ ) being the optimal control, we have:
H(X∗ , λ∗ , u∗ , θ ∗ ) = min H(X, λ, u, θ)
(31)
u(t),θ(t)

From (30) we see that we can always set the control θj (t) to
ensure that sgnηjy (t) sin(θj (t) + ψj (t)) < 0. Hence, recalling
that 0 ≤ uj (t) ≤ 1,
u∗j (t) = 1
(32)
and
sin(θj∗ (t) + ψj (t)) = 1 if µyj (t) < 0
sin(θj∗ (t) + ψj (t)) = −1 if µyj (t) > 0
(33)
Following the Hamiltonian definition in (29) we have:
∂H
= −ηjx (t)uj (t) sin θj (t) + ηjy (t)uj (t) cos θj (t) (34)
∂θj
∂H
and setting ∂θ
= 0 the optimal heading θj∗ (t) should satisfy:
j
ηjy (t)
(35)
ηjx (t)
Since u∗j (t) = 1, we only need to evaluate θj∗ (t) for all
t ∈ [0, T ]. This is accomplished by discretizing the problem
in time and numerically solving a TPBVP with a forward
integration of the state and a backward integration of the
costate. Solving this problem becomes intractable as the
number of agents and targets grows. The fact that we are
dealing with a hybrid dynamic system further complicates
the solution of a TPBVP. On the other hand, it enables us to
make use of Infinitesimal Perturbation Analysis (IPA) [24]
to carry out the parametric trajectory optimization process
discussed in the next section. In particular, we propose a
parameterization of agent trajectories allowing us to utilize
IPA to obtain an unbiased estimate for the objective function
gradient with respect to the trajectory parameters.
tan θj∗ (t) =

A. Agent Trajectory Parameterization and Optimization
The key idea is to represent each agent’s trajectory through
general parametric equations
sxj (t) = f (Θj , ρj (t)),
syj (t) = g(Θj , ρj (t))
(36)
where the function ρj (t) controls the position of the agent
on its trajectory at time t and Θj is a vector of parameters
controlling the shape and location of the agent j trajectory.
Let Θ = [Θ1 , . . . , ΘN ]. We now replace problem P1 in (25)
by problem
Z P2:
1 Th
qJ1 (Θ, t) − (1 − q)J2 (Θ, t) + J3 (Θ, t)
min
Θ∈FΘ T 0
(37)
i
+ J4 (Θ, t) dt + Jf (Θ, T )
where we return to allowing arbitrary stochastic data arrival
processes {σi (t)} so that P2 is a parametric stochastic
optimization problem with the feasible parameter set FΘ
appropriately defined depending on (36). The cost function
in (37) is written as
J(Θ, T ; X(Θ, 0)) = E[L(Θ, T ; X(Θ, 0))]
where L(Θ, T ; X(Θ, 0)) is a sample function defined over
[0, T ] and X(Θ, 0) is the initial value of the state vector. For
convenience, in the sequel we will use Li , i = 1, . . . , 4, and
Lf to denote sample functions of Ji , i = 1, . . . , 4, and Jf
respectively. Note that in (37) we suppress the dependence
of the four objective function components on the controls
u(t) and θ(t) and stress instead their dependence on the
parameter vector Θ.
In the rest of the paper, we will consider two families
of trajectories motivated by a similar approach used in the

multi-agent persistent monitoring problem in [23]: elliptical
trajectories and a more general Fourier series trajectory representation better suited for non-uniform target topologies.
The hybrid dynamics of the data harvesting system allow us
to apply the theory of IPA [24] to obtain on line the gradient
of the sample function L(Θ, T ; X(Θ, 0)) with respect to Θ.
The value of the IPA approach is twofold: (i) The sample
gradient ∇L(Θ, T ) can be obtained on line based on observable sample path data only, and (ii) ∇L(Θ, T ) is an unbiased
estimate of ∇J(Θ, T ) under mild technical conditions as
shown in [24]. Therefore, we can use ∇L(Θ, T ) in a standard
gradient-based stochastic optimization algorithm
Θl+1 = Θl − ν l ∇L(Θl , T ), l = 0, 1, . . .
(38)
to converge (at least locally) to an optimal parameter vector
Θ∗ with a proper selection of a step-size sequence {ν l } [27].
We emphasize that this process is carried out on line, i.e.,
the gradient is evaluated by observing a trajectory with given
Θ over [0, T ] and is iteratively adjusted until convergence is
attained.
1) IPA Calculus Review and Implementation: Based on
the events defined earlier, we will specify event time derivative and state derivative dynamics for each mode of the
hybrid system. In this process, we will use the IPA notation
from [24] so that τk is the kth event time in an observed
dX
0
k
sample path of the hybrid system and τk0 = dτ
dΘ , X (t) = dΘ
are the Jacobian matrices of partial derivatives with respect
to all components of the controllable parameter vector Θ.
Throughout the analysis we will be using (·)0 to show such
derivatives. We will also use fk (t) = dX
dt to denote the state
dynamics in effect over an interevent time interval [τk , τk+1 ).
We review next the three fundamental IPA equations from
[24] based on which we will proceed.
First, events may be classified as exogenous or endogenous. An event is exogenous if its occurrence time is
independent of the parameter Θ, hence τk0 = 0. Otherwise, an
endogenous event takes place when a condition gk (Θ, X ) =
0 is satisfied, i.e., the state X (t) reaches a switching surface
described by gk (Θ, X ). Inthiscase, it is shown in
[24] that
−1
dgk
dgk 0 − 
−
0
0
τk = −
fk (τk )
gk +
X (τk )
(39)
dX
dX
∂gk
−
as long as ∂X fk (τk ) 6= 0. It is also shown in [24] that the
state derivative X 0 (t) satisfies
d 0
dfk 0
X (t) =
X (t) + fk0 (t), t ∈ [τk , τk+1 )
(40)
dt0 +
dX
−
+
0 −
0
X (τk ) = X (τk ) + [fk−1 (τk ) − fk (τk )]τk
(41)
Then, X 0 (t) for t ∈ [τk , τk+1 ) isZcalculated through
t
d 0
X 0 (t) = X 0 (τk+ ) +
X (t)dt
(42)
τk dt
Table I contains all possible endogenous event types for
our hybrid system. To these, we add exogenous events κi ,
i = 1, ..., M , to allow for possible discontinuities (jumps)
in the random processes {σi (t)} which affect the sign of
σi (t) − µij pij (t) in (5). We will use the notation e(τk ) to
denote the event type occurring at t = τk with e(τk ) ∈ E,
the event set consisting of all endogenous and exogenous
events. Finally, we make the following assumption which is
needed in guaranteeing the unbiasedness of the IPA gradient

estimates: (A6) Two events occur at the same time w.p. 0
unless one is directly caused by the other.
2) Objective Function Gradient: The sample function
gradient ∇L(Θ, T ) needed in (38) is obtained from (37)
assuming a total number of K events over [0 T ] with
τK+1 = T and τ0 = 0:
∇L(Θ, T ; X(Θ; 0)) =
Z
1 h T
∇
qL1 (Θ, t) − (1 − q)L2 (Θ, t) + L3 (Θ, t)
T
0
 i
+ L4 (Θ, t) dt + ∇Lf (Θ, T )
K Z
1 h X τk+1 
= ∇
qL1 (Θ, t) − (1 − q)L2 (Θ, t) + L3 (Θ, t)
T
k=0 τk
 i
+ L4 (Θ, t) dt + ∇Lf (Θ, T )
Z
K

1 h X  τk+1
0
=
q
∇L1 (Θ, t)dt + L1 (Θ, τk+1 )τk+1
− L1 (Θ, τk )τk0
T
τk
k=0
 Z τk+1

0
− (1 − q)
∇L2 (Θ, t)dt + L2 (Θ, τk+1 )τk+1
− L2 (Θ, τk )τk0
τk
 Z τk+1

0
+
∇L3 (Θ, t)dt + L3 (Θ, τk+1 )τk+1
− L3 (Θ, τk )τk0
τk
 Z τk+1
i
0
+
∇L4 (Θ, t)dt + L4 (Θ, τk+1 )τk+1
− L4 (Θ, τk )τk0

(43)

τk

+ ∇Lf (Θ, T )
K Z
1 h X τk+1 
q∇L1 (Θ, t) − (1 − q)∇L2 (Θ, t) + ∇L3 (Θ, t)
=
T
k=0 τk
 i
+ ∇L4 (Θ, t) dt + ∇Lf (Θ, T )

The last step follows from the continuity of the state variables
which causes adjacent limit terms in the sum to cancel out.
Therefore, ∇L(Θ, T ) does not have any direct dependence
on any τk0 ; this dependence is indirect through the state
derivatives involved in the four individual gradient terms.
Referring to (10), the first term in (43) involves ∇L1 (Θ, t)
which is as a sum of Xi0 (t) derivatives. Similarly, ∇L2 (Θ, t)
is a sum of Yi0 (t) derivatives and ∇Lf (Θ, T ) requires only
0
Zij
(T ). The third term, ∇L3 (Θ, t), requires derivatives of
Ij (t) in (13) which depend on the derivatives of the max
function in (9) and the agent state derivatives s0j (t) with
respect to Θ. The term ∇L4 (Θ, t) needs the values of Xi0 (t)
0
and Zij
(t). The gradients of the last two terms are derived
in the appendix. Possible discontinuities in these derivatives
occur when any of the last four events in Table I takes place.
In summary, the evaluation of (43) requires the state
0
derivatives Xi0 (t), Zij
(t), Yi0 (t), and s0j (t). The latter are
easily obtained for any specific choice of f and g in (36)
and are shown in Appendix II. The former require a rather
laborious use of (39)-(41) which, however, reduces to a
simple set of state derivative dynamics as shown next.
Proposition 2: : After an event occurrence at t = τk , the
0
state derivatives Xi0 (τk+ ), Yi0 (τk+ ), Zij
(τk+ ), with respect to
the controllable
 parameter Θ satisfy the following: 0
if e(τk ) = ξi
 0
+
Xi0 (τk− ) − µil pil (τk )τk0 if e(τk ) = δij
Xi0 (τk+ ) =
 0 −
Xi (τk )
otherwise
where l 6= j with pil (τk ) > 0 if such l exists and
τk0 =

∂dij (sj ) 0
∂sj sj

Yi0 (τk+ ) =

−1
∂dij (sj )
.
∂sj ṡj (τk )
 0 −
0
Yi (τk ) + Zij
(τk− )
Yi0 (τk− )

0
if e(τk ) = ζij
otherwise


0
if e(τk ) = ζij
 0
−
+
0
0
0 −
0
Z (τk ) + Xi (τk ) if e(τk ) = ξi
Zij (τk ) =
 ij
0
Zij
(τk− )
otherwise
0
where e(τk ) = ξi occurs when j is connected to target i.
Proof: : See (76), (87), (95), (93), (79), (88), (90), (82)
in Appendix IV.
This result shows that only three of the events in E can
actually cause discontinuous changes to the state derivatives.
Further, note that Xi0 (t) is reset to zero after a ξi0 event.
0
Moreover, when such an event occurs, note that Zij
(t)
0
0
0
is coupled to Xi (t). Similarly for Zij (t) and Yi (t) when
0
event ζij
occurs, showing that perturbations in Θ can only
propagate to an adjacent queue when that queue is emptied.
−
−
Proposition 3: : The state derivatives Xi0 (τk+1
), Yi0 (τk+1
)
with respect to the controllable parameter Θ satisfy the
following after
 an event occurrence at t = τk :
0
if e(τk ) = ξi0
−
R τk+1
Xi0 (τk+1
)=
0 +
0
Xi (τk ) − τk µij pij (u)du otherwise
Z τk+1
−
Yi0 (τk+1
) = Yi0 (τk+ ) +
βi0 (u)du
τk

where j is such that pij (t) > 0, t ∈ [τk , τk+1 ).
Proof: : See (75), (78) and (80) in Appendix IV.
+
0
Proposition 4: : The state derivatives Zij
(τk+1
) with respect to the controllable parameter Θ satisfy the following
after an event occurrence at t = τk :
i- If j is connected
 0 to+target i,
+
0
Zij (τk ) R
if e(τk ) = ξi0 , ζij
or δij
−
0
Zij (τk+1 ) =
τ
k+1
0
Zij
(τk+ ) + τk µij p0ij (u)du otherwise
ii- If j is connected to B with ZZij (τk ) > 0,
−
0
0
Zij
(τk+1
) = Zij
(τk+ ) −

τk+1

βij p0Bj (u)du

τk
0
Zij
(τk+ ).

−
0
iii- Otherwise, Zij
(τk+1
)=
Proof: : See (83), (84), (91) and (98) in Appendix IV.
0
Corollary 5: The state derivatives Xi0 (t), Zij
(t), Yi0 (t)
with respect to the controllable parameter Θ are independent
of the random data arrival processes {σi (t)}, i = 1, . . . , M .
Proof: : Follows directly from the three Propositions.

There are a few important consequences of these results.
First, as the Corollary asserts, one can apply IPA regardless
of the characteristics of the random processes {σi (t)}. This
robustness property does not mean that these processes do
0
not affect the values of the Xi0 (t), Zij
(t), Yi0 (t); this happens
through the values of the event times τk , k = 1, 2, . . ., which
are observable and enter the computation of these derivatives
as seen above.
Second, the IPA estimation process is event-driven:
0
Xi0 (τk+ ), Yi0 (τk+ ), Zij
(τk+ ) are evaluated at event times
and then used as initial conditions for the evaluations of
−
−
−
0
Xi0 (τk+1
), Yi0 (τk+1
), Zij
(τk+1
) along with the integrals
appearing in Propositions 2,3 which can also be evaluated
at t = τk+1 . Consequently, this approach is scalable in the
number of events in the system as the number of agents and
targets increases.
Third, despite the elaborate derivations in the Appendix,
the actual implementation reflected by the three Propositions

is simple. Finally, returning to (43), note that the integrals
involving ∇L1 (Θ, t), ∇L2 (Θ, t) are directly obtained from
Xi0 (t), Yi0 (t), the integral involving ∇L3 (Θ, t) is obtained
from straightforward differentiation of (13), and the final
0
term is obtained from Zij
(T ).
3) Objective Function Optimization: This is carried out
using (38) with an appropriate diminishing step size sequence.
Elliptical Trajectories: Elliptical trajectories are described by their center coordinates, minor and major axes and
orientation. Agent j’s position sj (t) = [sxj (t), syj (t)] follows
the general parametric equation of the ellipse:
sxj (t) = Aj + aj cos ρj (t) cos φj − bj sin ρj (t) sin φj
syj (t) = Bj + aj cos ρj (t) sin φj + bj sin ρj (t) cos φj

(44)

Here, Θj = [Aj , Bj , aj , bj , φj ] where Aj , Bj are the coordinates of the center, aj and bj are the major and minor axis
respectively while φj ∈ [0, π) is the ellipse orientation which
is defined as the angle between the x axis and the major axis
of the ellipse. The time dependent parameter ρj (t) is the
eccentric anomaly of the ellipse. Since the agent is moving
with constant speed of 1 on this trajectory from (32), we
have ṡxj (t)2 + ṡyj (t)2 = 1 which gives


ρ̇j (t) = 



2

− 21

a sin ρj (t) cos φj + bj cos ρj (t) sin φj


2 
+ a sin ρj (t) sin φj − bj cos ρj (t) cos φj

In the data harvesting problem, trajectories that do not pass
through the base are inadmissible since there is no delivery
of data. Therefore, we add a constraint to force the ellipse
to pass through wB = [wBx , wBy ] where:
wBx =Aj + aj cos ρj (t) cos φj − bj sin ρj (t) sin φj
(45)
wBy =Bj + aj cos ρj (t) sin φj + bj sin ρj (t) cos φj
Using the fact that sin2 ρ(t) + cos2 ρ(t) = 1 we define a
quadratic constraint term added to J(Θ, T ; X(Θ, 0)) with a
sufficiently large multiplier. This can ensure the optimal path
passes through the base location wB . We define Cj (Θj ):
2
Cj (Θj ) = 1 − fj1 cos2 φj − fj2 sin2 φj − fj3 sin 2φj (46)
wy −Bj 2
wx −Aj 2
wx −Aj 2
B
B
+ Bbj
, fj2 =
+
where fj1 =
aj
bj
y
2
2
x
y
w −Bj 2
(b
−a
)(w
−A
)(w
−B
)
j
j
B
, fj3 = j j Ba2 b2 B
.
aj
j j
Multiple visits to the base may be needed during the
mission time [0, T ]. We can capture this by allowing an agent
trajectory to consist of a sequence of admissible ellipses.
For each agent, we define Ej as the number of ellipses in
its trajectory. The parameter vector Θκj with κ = 1, . . . , Ej ,
defines the κth ellipse in agent j’s trajectory and Tjκ is the
time that agent j completes ellipse κ. Therefore, the location
of each agent is described through κ during [Tjκ−1 , Tjκ ]
where Tj0 = 0. Since we cannot optimize over all possible
Ej for all agents, an iterative process needs to be performed
in order to find the optimal number of segments in each
agent’s trajectory. At each step, we fix Ej and find the optimal
trajectory with that many segments. The process is stopped
once the optimal trajectory with Ej segments is no better
than the optimal one with Ej − 1 segments (obviously, this
is generally not a globally optimal solution). We can now
formulate the parametric optimization problem P2e where

(a) TPBVP Trajectories for Case I
Fig. 5.

(b) Elliptical Trajectories for case I

Simulation results for the two-target and two-agent case

E

Θj = [Θ1j , . . . , Θj j ] and Θ = [Θ1 , . . . , ΘN ]:
1
min Je =
Θ∈FΘ
T

Z

T

constraints in Zthis
case:
T h

h
qJ1 (Θ, t) − (1 − q)J2 (Θ, t) + J3 (Θ, t)

0

min Jf =

Θ∈FΘ

(47)

N
i
X
+ J4 (Θ, t) dt + MC
Cj (Θj ) + Jf (Θ, T )
j=1

where MC is a large multiplier. The evaluation of ∇Cj is
straightforward and does not depend on any event (details
are shown in Appendix II).
Fourier Series Trajectories: The elliptical trajectories are
limited in shape and may not be able to cover many targets in
a mission space. Thus, we next parameterize the trajectories
using a Fourier series representation of closed curves [28].
Using a Fourier series function for f and g in (36), agent j’s
trajectory can be described as follows with base frequencies
fjx and fjy :
PΓxj
sxj (t) = a0,j + n=1
an,j sin(2πnfjx ρj (t) + φxn,j )
y
(48)
P
Γj
syj (t) = b0,j + n=1
bn,j sin(2πnfjy ρj (t) + φyn,j )
The parameter ρ(t) ∈ [0, 2π], similar to elliptical trajectories,
represents the position of the agent along the trajectory. In
this case, forcing a Fourier series curve to pass through the
base is easier. For simplicity, we assume a trajectory to start
at the base and set sxj (0) = wBx , syj (0) = wBy . Assuming
ρ(0) = 0, with no loss of generality, we can calculate the
zero frequency terms by means of the remaining parameters:
Γy
Γx
j
j
X
X
a0,j = wBx −
an,j sin(φxn,j ), b0,j = wBy −
bn,j sin(φyn,j )
n=1

n=1

The parameter vector for agent j is Θj
=
[fjx , a0,j , . . . , aΓxj , b0,j , . . . , bΓyj , φ1,j , . . . , φΓxj , ξ1,j , . . . , ξΓyj ]
and Θ = [Θ1 , . . . , ΘN ]. Note that the shape of the curve is
fully captured by the ratio fjx /fjy , so that one of these two
parameters can be kept constant. For the Fourier trajectories,
the fact that u∗j = 1 allows us to calculate ρ̇j (t) as follows:
Γx


ρ̇j (t) =

1
2π








fjx

j
X

!2
an,j n cos(2πfjx ρj (t)

+

−1/2

φxn,j )




!2 
Γx
j

X

bn,j n cos(2πfjy ρj (t) + φyn,j )
+ fjy
n=1

n=1

(c) Fourier Trajectories for case I

Problem P2f is the same as P2 but there are no additional

1
T

qJ1 (Θ, t) − (1 − q)J2 (Θ, t) + J3 (Θ, t)
(49)
i
+ J4 (Θ, t) dt + Jf (T )
0

IV. N UMERICAL R ESULTS
In this section numerical results are presented to illustrate
our approach. The mission space S is considered to be
[0, 10] × [0, 10] in all cases. The first case we consider is
a small mission to obtain the TPBVP results and confirm
the fact that it is not scalable to bigger problems.
In Case I we consider a two-target, two-agent setting. We
assume deterministic arrival processes with σi = 0.5 for
all i. For (3) and (4) we have used p(w, v) = max(0, 1 −
d(w,v)
) where r is the corresponding value of rij or rBj .
r
We set µij = 100 and βij = 500 for all i and j. Other
parameters used are q = 0.5, rij = rBj = 0.5 and T =
20. The trajectory comparison from TPBVP, Elliptical and
Fourier parametric solutions is shown in Figs. 5(a), 5(b), 5(c).
In each figure, the trajectories are shown in the top part, while
the actual objective function convergence behavior is hown
in the middle graph. The lower graph shows the total amount
of data at targets at any time (in blue) and the total amount
of data at the base (in green).
In the TPBVP results, the main limitation is the number
of the time steps in the discretization of the interval [0, T ],
since the number of control values grows with it. To bring
this into prospective, for this sample problem with T = 20
we considered 300 time steps, i.e., 300 values for the heading
of each agent need to be calculated, which brings the total
number of controls to 600. In contrast, for the same problem
the total number of controls for the elliptical trajectories are
10 parameters and for the Fourier trajectories it is 28. This
explains why the TPBVP cannot be a viable solution for
larger values of T . Note that, in this scenario the total time is
only 20 time steps so we can obtain a TPBVP solution. This,
however, causes a poor representation for the parameterized
trajectories which are spending significant time outside the
convex hull since they have not converged after only 20 time
steps.

the same performance by the optimization algorithm in the
stochastic setting.

V. C OMPARISON WITH A G RAPH BASED A LGORITHM

(a) Elliptical Trajectories for Case (b) Fourier Trajectories for Case II
II
Fig. 6. Simulation results for 9-target and two-agent case
TABLE II
R ESULT C OMPARISON FOR C ASE I

Method
TPBVP
Elliptical
Fourier

J∗
0.272
0.255
0.202

J1∗
0.098
0.092
0.089

−J2∗
-0.038
-0.095
-0.095

In Table II, the actual values for J ∗ , J1∗ , J2∗ are shown for
the three different trajectories of Fig. 5(a).,5(b),5(c). Note
that the objective is to minimize J by minimizing J1 − J2 .
Next, in Case II we consider 9 targets and 2 agents. The
base is located at the center of the mission space. We have
σi (t) = 0.5, µij = 50 and βij = 500 for all i and j. Other
parameters used are q = 0.5, rij = 0.55, rBj = 0.65 and
T = 50. In Fig. 6(a) the solution with two ellipses in each
agent’s trajectory is shown. As can be seen the trajectory
correctly finds all the target locations and empties the target
queues periodically. Fig. 6(b) shows the Fourier trajectories.
The two graphs on the bottom show the objective function
value and instantaneous total content at targets and base. The
results for Case II are summarized in Table III.
Case III has 12 targets that are uniformly distributed in
the mission space. Here, we try to examine the robustness
of our approach with respect to the arrival rate process at
targets. We use the same parameters as in case II and solve
the problem for deterministic σi (t) = 0.5 using the elliptical
trajectories and Fourier trajectories. The same mission is
simulated assuming that σi (t) is a stochastic process with
piecewise linear arrival rate. The average arrival rate is
kept at 0.5. The results in Figs. 7(a),7(b),7(c) and Table IV
show that the Fourier parametric trajectories achieve almost
TABLE III
R ESULTS C OMPARISON FOR C ASE II

Method
Elliptical
Fourier

J∗
0.19
0.18

J1∗
0.090
0.069

−J2∗
-0.124
-0.117

TABLE IV
R ESULTS C OMPARISON FOR C ASE III

Method
Elliptical
Fourier
Fourier (Stochastic Arrival)

J∗
0.35
0.23
0.23

J1∗
0.12
0.09
0.13

−J2∗
-0.09
-0.1
-0.13

We begin with the observation that the final parametric
trajectories provide a sequence of targets visits, similar to
the functionality of a tour selection algorithm that uses the
underlying graph topology of the mission space to determine
such sequences.
We have compared the results of our approach with
a graph topology algorithm called Path Splitter Heuristic
(PSH) developed in [12]. This algorithm starts with the
best Hamiltonian sequence and then uses a heuristic method
to divide the Hamiltonian tour into several sub-tours that
go through a few targets and then return to the base. The
algorithm then provides a sequence of these sub-tours for
each agent. We compare the sequences from Case I and Case
II in both elliptical and Fourier trajectories and results are
shown in Tables V and VI. For a fair comparison, we adopt
each sequence and apply it with the system dynamics in our
model, i.e., an agent can collect the data once within range
of a target and the data collection does not happen instantaneously. This, however, is not the basic modeling assumption
used in PSH, where agents pick up all the data at the target
instantaneously once at the target location. We compare the
sum total of data at targets and the base for T = 200. A
larger value of T is used for this comparison in order to
approximate infinite time results. These sequences are shown
in Figs. 8(a), 8(b) and 8(c) where each color represents one
agent trajectory. We can see from these comparisons that
in the graph-based approach targets are completely divided
between agents. This generates a spatial partitioning, giving
each agent full responsibility for a set of targets. However,
in the trajectory planning results, in most cases, we see a
temporal partitioning where agents can visit the same targets
but at different times of the mission. This clearly allows
for more robustness with respect to potential agent failures
or changes in an agent’s operation. Moreover, even though
the computational complexity of the PSH algorithm and the
parametric trajectory optimization approach are comparable,
methods such as PSH need to re-solve the complete problem
each time a new target may appear in the mission space.
In contrast, the on-line event-driven parametric optimization
process is a methodology designed to adapt to targets which
may randomly appear in the mission space. These results are
not necessarily aiming to prove performance enhancement
but to put the two approaches into contrast in terms of
suitability for online and offline applications. Also, PSH
agent trajectories consist of straight line segments. These are
not physically realizable, given limitations on the motion of
agents which must smoothly turn direction from one target to
the next. Whereas the parametrical trajectories are easier to
realize by most of the agents given these motion limitations.

(a) Elliptical Trajectories for case III
Fig. 7.

(a) PSH Sequences for case III

(b) Fourier Trajectories for case III

(c) Fourier Trajectories for case III
with Stochastic Arrival
Simulation results for the 12-target and two-agent case

(b) Elliptical Sequences for case III
Fig. 8. Sequence Comparison

TABLE V
R ESULT C OMPARISON WITH PSH FOR C ASE II

Method
PSH Sequence
Elliptical Sequence
Fourier Sequence

J1∗

−J2∗

0.023
0.027
0.024

-0.22
-0.21
-0.21

TABLE VI
R ESULT C OMPARISON WITH PSH FOR C ASE III

Method
PSH Sequence
Elliptical Sequence
Fourier Sequence

J1∗
0.0257
0.0304
0.0212

−J2∗
-0.21
-0.199
-0.21

VI. C ONCLUSIONS
We have demonstrated a new event-driven methodology
for on-line trajectory optimization with application in the
data harvesting problem. We proposed a new performance
measure that addresses the event excitation problem in eventdriven controllers. The proposed optimal control problem is
then formulated as a parametric trajectory optimization utilizing general function families which can be subsequently optimized on line through the use of Infinitesimal Perturbation
Analysis (IPA). Several numerical results are provided for
the case of elliptical and Fourier series trajectories and some
properties of the solution are identified, including robustness
with respect to the stochastic data generation processes and
scalability in the size of the event set characterizing the

(c) Fourier Sequences for case III

underlying hybrid dynamic system.
Although the proposed methodology is focused on applying the event-driven optimization approach to the data
harvesting problem, it should be noted that the new metric
which was introduced in [25] to ensure event excitation
allows us to generalize the methodology to other applications
as well. The new metric introduces a potential field or density
map over the entire mission space. This can viewed as a
probability distribution of potential targets in problems where
the exact locations of targets are unknown. Used as a prior
distribution, it can be improved while the agents move within
the mission space and gather more information. In addition,
this density can be dynamically changing if the targets are
moving and their location changes with time assuming some
prior information about the target paths. This provides a tool
to apply agent trajectory optimization in a much broader
range of problems tracking moving points of interest.
A PPENDIX I
P ROOFS
Proposition 1:
Proof: Fixing Xi (t) = xi (t), from (18) we have
Z
Z
Z X
M
M
X
αi xi (t)
xi (t)
dw
=
α
dw
R(w, t) =
i
+
+
d
C
C
C i=1 di (w)
i (w)
i=1
R
(t)
To evaluate C dx+i(w)
for each target i, we first look at the
i
case of a single target in a 2D space and temporarily replace

3. Target i is in the interior of C but C(wi ) is not completely
in the interior of C: This is the case of the yellow target in
Fig. 9. Carrying out the integration for the appropriate limits,
since C(wi ) is not fully contained in C:
Z 2π Z Λ(θ)
xi (t)
xi (t)
dw
=
drdθ
+
+
d
(w)
d
0
0
C i
i (r, θ)
Z 2π Z r(θ)
Z θ2 Z Λ(θ)
xi (t)
xi (t)
=
drdθ +
drdθ
r
r
i
0
ri
0
θ1
Z
Z θ2
 2π r(θ)

Λ(θ)
log(
= xi (t)
dθ +
)dθ
(54)
ri
ri
0
θ1
Here, Λ(θ) is only defined for θ ∈ [θ1 θ2 ] as shown in Fig.
9 for the yellow target. Again, since Λ(θ) ≥ ri , it follows
that log( Λ(θ)
ri ) > 0 and with r(θ) > 0 the term in brackets
above is positive.ZWe can define Zci in (19) as
θ2
 2π r(θ)

Λ(θ)
ci = αi
dθ +
log(
)dθ
(55)
ri
ri
0
θ1
Z

Fig. 9.

One Target R(w, t) Calculation

C by a disk with radius Λ around the target (black circle with
radius Λ in Fig. 9). We can now calculate the integral above
using polar
Z coordinates:Z 2π Z Λ
xi (t)
xi (t)
drdθ
dw =
+
max(r
d
(w)
i , r)
0
0
C i
Z 2π Z ri
Z 2π Z Λ
xi (t)
xi (t)
=
drdθ
drdθ +
ri
r
0
ri
0
0

Λ 
= xi (t) 2π 1 + log( )
ri
Since C is actually the convex Rhull of all targets, we will
(t)
dw where C is the
use the same idea to calculate C dx+i(w)
i
convex hull by considering three separate cases depending
on the target location.
1. Target i and C(wi ) are completely in the interior of C:
This is shown in Fig. 9 for the red target. Using the same
polar coordinate for each θ, we define Λ(θ) to be the distance
of theZtarget to the edge
of ZC in the direction of θ and get:
Z 2π
Λ(θ)
xi (t)
xi (t)
dw =
drdθ
+
+
d
(w)
d
C i
0
0
i (r, θ)
Z 2π Z ri
Z 2π Z Λ(θ)
xi (t)
xi (t)
=
drdθ
drdθ +
r
r
i
0
0
0
ri
Z 2π


Λ(θ)
= xi (t) 2π +
)dθ
(50)
log(
ri
0
The second part in (50) is calculated knowing that Λ(θ) ≥
ri . This means log( Λ(θ)
ri ) > 0 and the term in brackets is
positive. We can then defineZ ci in (19) as
2π


Λ(θ)
ci = αi 2π +
log(
)dθ
(51)
r
i
0
2. Target i is on an edge of C: This is shown in Fig. 9 for the
green target. Since C(wi ) is not entirely contained within C,
we have
Z
Z 2π Z Λ(θ)
xi (t)
xi (t)
dw =
drdθ
+
+
d
(w)
d
0
0
C i
i (r, θ)
Z θ2 Z r i
Z θ2 Z Λ(θ)
xi (t)
xi (t)
=
drdθ +
drdθ
ri
r
θ1
0
θ1
ri
Z θ2


Λ(θ)
= xi (t) θ2 − θ1 +
log(
)dθ
(52)
ri
θ1
The second part in (52) is calculated knowing that for θ ∈
[θ1 θ2 ] we have Λ(θ) ≥ ri . Therefore, log( Λ(θ)
ri ) > 0 and by
definition θ2 −θ1 > 0 so that the term in brackets is positive.
We can then define ci in (19) Zas
θ2


Λ(θ)
ci = αi θ2 − θ1 +
log(
)dθ
(53)
ri
θ1

A PPENDIX II
E LLIPTICAL T RAJECTORIES
In order to calculate the IPA derivatives we need the
derivatives of state variable associated with agents with
respect to all entries in the parameter vector Θj =
[Aj , Bj , aj , bj , φj ] for all agents j. These derivatives do
not depend on the events in the system, since the agent
trajectories are fixed at each iteration. For now we assume
Ej = 1 for all j = 1, . . . , N hence, we drop the superscript.
We have:
∂sxj
∂sxj
= 1,
=0
(56)
∂Aj
∂Bj
x
x
∂sj
∂sj
= cos ρj (t) cos φj ,
= − sin ρj (t) sin φj
∂aj
∂bj
(57)
∂sxj
= −aj cos ρj (t) sin φj − bj sin ρj (t) cos φj
(58)
∂φj
y
y
∂sj
∂sj
= 0,
=1
(59)
∂Aj
∂Bj
y
y
∂sj
∂sj
= cos ρj (t) sin φj ,
= sin ρj (t) cos φj (60)
∂aj
∂bj
y
∂sj
= aj cos ρj (t) cos φj − bj sin ρj (t) sin φj
(61)
∂φj
The time derivative of the position state variables are calculated as follows:
ṡxj (t) = −aj ρ̇j (t) sin ρj (t) cos φj + bj ρ̇j (t) cos ρj (t) sin φj (62)
ṡyj (t) = −aj ρ̇j (t) sin ρj (t) sin φj + bj ρ̇j (t) cos ρj (t) cos φj (63)
The gradient of the last term in the Je in (47) needs to be
∂C
calculated separately. We have for j 6= l, ∂Θjl = 0 and for
j = l:
∂fj1
∂fj2
∂fj3 
∂Cj
= 2Cj − cos2 φj
− sin2 φj
− sin 2φj
∂Aj
∂Aj
∂Aj
∂Aj
1
2
∂f
∂f
∂fj3 
∂Cj
j
j
2
= 2Cj − cos2 φj
− sin φj
− sin 2φj
∂Bj
∂Bj
∂Bj
∂Bj
1
2
3
∂f
∂f
∂f
∂Cj
j
j
j 
= 2Cj − cos2 φj
− sin2 φj
− sin 2φj
∂aj
∂aj
∂aj
∂aj

∂fj2
∂fj3 
∂fj1
∂Cj
= 2Cj − cos2 φj
− sin2 φj
− sin 2φj
∂bj
∂bj
∂bj
∂bj

∂Cj
1
2
3
= 2Cj (fj − fj ) sin 2φj − 2fj cos 2φj
∂φj
where1
 wx − A 
 wy − B 
∂fj1
∂fj
j
j
,
= −2 B 2
= −2 B 2
∂Aj
aj
∂Bj
bj
1
 (wx − A )2 
 (wy − B )2 
∂f
∂fj1
j
j
j
B
B
,
= −2
= −2
∂aj
a3j
∂bj
b3j
 wx − A  ∂f 2
 wy − B 
∂fj2
j
j
j
,
= −2 B 2
= −2 B 2
∂Aj
bj
∂Bj
aj
 (wy − B )2  ∂f 2
 (wx − A )2 
∂fj2
j
j
j
B
B
,
= −2
= −2
∂aj
a3j
∂bj
b3j
 (b2 − a2 )(wy − Bj ) 
∂fj3
j
j
B
=−
∂Aj
a2j b2j
 (b2 − a2 )(wx − Aj ) 
∂fj3
j
j
B
=−
∂Bj
a2j b2j
 (wx − A )(wy − B ) 
∂fj3
j
j
B
B
= −2
∂aj
a3j
 (wx − A )(wy − B ) 
∂fj3
j
j
B
B
=2
∂bj
b3j

A PPENDIX III
F OURIER S ERIES T RAJECTORIES
In the Fourier parametric trajectories the agent state derivative is calculated as follows. The parameter vector is Θj =
[fjx , a0,j , . . . , aΓxj , b0,j , . . . , bΓyj , φ1,j , . . . , φΓxj , ξ1,j , . . . , ξΓyj ].
Thus, we have:
∂sxj
∂sxj
= 1,
=0
(64)
∂a0,j
∂b0,j
x
x
∂sj
∂sj
= sin(2πnfjx ρj (t) + φxn,j ),
= 0 (65)
∂an,j
∂bn,j
x
x
∂sj
∂sj
= an,j cos(2πnfjx ρj (t) + φxn,j )
= 0 (66)
x
∂φn,j
∂φyn,j
Γx

j
X
∂sxj
=
2πρ
(t)
an,j n cos(2πnfjx ρj (t) + φxn,j ), (67)
j
∂fjx
n=1
∂syj
∂syj
= 1,
=0
(68)
∂b0,j
∂a0,j
y
y
∂sj
∂sj
= sin(2πnfjy ρj (t) + φyn,j ),
= 0 (69)
∂bn,j
∂an,j
y
y
∂sj
∂sj
= bn,j cos(2πnfjy ρj (t) + φyn,j )
= 0 (70)
y
∂φn,j
∂φxn,j
y
∂sj
=0
(71)
∂fjx
The time derivative of the position state variables are calculated as follows: x
Γj
X
ṡxj (t) = ρ̇j (t)
2πnfjx an,j cos(2πnfjx ρj (t)+φxn,j ), (72)

n=1

Γy

ṡyj (t)

= ρ̇j (t)

j
X

2πnfjy bn,j cos(2πnfjy ρj (t) + φyn,j ), (73)

n=1

A PPENDIX IV
IPA EVENTS AND DERIVATIVES
In this section, we derive all event time derivatives and
state derivatives with respect to the controllable parameter
Θ for each event by applying the IPA equations.
1. Event ξi0 : This event causes a transition from Xi (t) > 0,
t < τk to Xi (t) = 0, t ≥ τk . The switching function is
∂gk
= 1. From (39) and (5):
gk (Θ, X) = Xi so ∂X
i
 ∂g
−1 
∂gk 0 − 
k
τk0 = −
fk (τk− )
X (τ )
gk0 +
∂Xi
∂Xi i k
(74)
Xi0 (τk− )
=−
σi (τk ) − µij pij (τk )
where agent j is the one connected to i at t = τk and we
have used the assumption that two events occur at the same
time w.p. 0, hence σi (τk− ) = σi (τk ). From (40)-(41), since
Ẋi (t) = 0, for τk ≤ t < τk+1 :
d 0
∂ Ẋi (t) 0
Xi (t) =
Xi (t) + Ẋi0 (t) = 0
(75)
dt
∂X
(t)
i
h

i
Xi0 (τk+ ) = Xi0 (τk− ) + σi (τk ) − µij pij (τk ) − 0 τk 0


Xi0 (τk− ) σi (τk ) − µij pij (τk )
= Xi0 (τk− ) −
=0
σi (τk ) − µij pij (τk )
(76)
For Xr (t), r 6= i, the dynamics of Xr (t) in (5) are unaffected
and we have:
Xr0 (τk+ ) = Xr0 (τk− )
(77)
If Xr (τk ) > 0 and agent l is connected to it, then
∂ Ẋr (t) 0
d 0
Xr (t) =
X (t) + Ẋr0 (t)
(78)
dt
∂Xr (t) r
= σr0 (t) − µrl p0rl (τk ) = −µrl p0rl (t)
and if Xr (t) = 0 in [τk , τk+1 ] or if no agents are connected
d
to i, then dt
Xr0 (t) = 0.
For Yr (t), r = 1, . . . , M , the dynamics of Yr (t) in (8) are
not affected by the event ξi0 at τk , hence
Yr0 (τk+ ) = Yr0 (τk− )
(79)
and since Ẏr (t) = βr (t), for τk ≤ t < τk+1 :
d 0
∂ Ẏr (t) 0
Y (t) =
Y (t) + Ẏr0 (t) = βr0 (t)
(80)
dt r
∂Yr (t) r
For Zij (t), we must have Zij (τk ) > 0 since Xi (τk− ) > 0,
hence µ̃ij (τk− ) > 0 and from
h (41):
i
0
0
Zij
(τk+ ) = Zij
(τk− ) + Żij (τk− ) − Żij (τk+ ) τk0
h
i
(81)
0
= Zij
(τk− ) + µ̃ij (τk− ) − µ̃ij (τk+ ) pij (τk )τk0
Since Xi (τk− ) > 0, from (6) we have µ̃ij (τk− ) = µij .
At τk+ , j remains connected to target i with µ̃ij (τk+ ) =
σi (τk+ )/pij (τk ) = σi (τk )/pij (τk )hand we get
i
−Xi0 (τk− ) µij pij (τk ) − σi (τk )
0
0
Zij
(τk+ ) = Zij
(τk− ) +
(82)
σi (τk ) − µij pij (τk )
−
0
0 −
= Zij (τk ) + Xi (τk )

From (40) for τk ≤ t < τk+1 :
∂ Żij (t) 0
d 0
0
Z (t) =
Z (t) + Żij
(t)
dt ij
∂Zij (t) ij
(83)


0
0
0
= Żij (t) = µ̃ij (t)Pij (t) − βij PBj (t)
Since µ̃ij (t) = σi (t)/pij (t) for the agent which remains
connected to target i after this event, it follows that
∂
∂Θ [µ̃ij (t)pij (t)] = 0. Moreover, pBj (t) = 0 by our assumption that agents cannot be within range of the base and targets
at the same time and we get
d 0
Z (t) = 0
(84)
dt ij
Otherwise, for r 6= j, we have µ̃ir (t) = 0 and we get:
d 0
0
(t)
(85)
Z (t) = −βir pBr
dt ir
+
−
0
0
Finally, for Zrj (t), r 6= i we have Zrj (τk ) = Zrj (τk ). If
d 0
Zrj (t) = 0. Otherwise, we
Zrj (t) = 0 in [τk , τk+1 ), then dt
d 0
get dt Zrj (t) from (83) with i replaced by r.
2. Event ξi+ : This event causes a transition from Xi (t) = 0,
t ≤ τk to Xi (t) > 0, t > τk . Note that this transition can
occur as an exogenous event when an empty queue Xi (t)
gets a new arrival in which case we simply have τk0 = 0
since the exogenous event is independent of the controllable
parameters. In the endogenous case, however, we have the
switching function gk (Θ, X) = σi (t) − µij pij (t) in which
agent j is connected to target i at t = τk . Assuming s0j (t) =
 ∂sxj ∂syj >
and ṡj = [ṡxj ṡyj ]> , from (39):
∂Θ ∂Θ
 ∂g

−1
k 0
τk 0 = −
sj (τk ) gk0 ṡj (τk )
(86)
∂sj
At τk we have σi (τk ) = µij pij (τk ). Therefore from (41):
Xi0 (τk+ ) = Xi0 (τk− ) + [Ẋi (τk− ) − Ẋi (τk+ )]τk 0


= Xi0 (τk− ) + 0 − σi (τk ) + µij pij (τk ) τk 0 = Xi0 (τk− )
(87)
Having Xi (t) > 0 in [τk , τk+1 ) we know Ẋi (t) = σi (t) −
d
µij pij (t) therefor, we can get dt
Xi0 (t) from (78) with r and
l replaced by i and j. For Xr (t), r 6= i, if Xr (τk ) > 0 and
agent l is connected to r then Ẋr (τk ) = σr (τk )−µrl prl (τk ),
therefor, we get Xr0 (τk+ ) from (77) while in [τk , τk+1 ) we
d
have dt
Xr0 (t) from (78). If Xr (τk ) = 0 or if no agent is
connected to target r, Ẋr (τk ) = 0. Thus, Xr0 (τk+ ) = Xr0 (τk− )
d
and dt
Xr0 (t) = 0.
For Yr (t), r = 1, . . . , M the dynamics of Yr (t) in (8) are
not affected by the event at τk hence, we can get Yr0 (τk+ )
d
Yr0 (t) in [τk , τk+1 ) from (79) and (80) respectively.
and dt
For Zij (t) assuming agent j is the one connected to target
i, we have:
h
i
0
0
Zij
(τk+ ) = Zij
(τk− ) + Żij (τk− ) − Żij (τk+ ) τk0
h
i
0
0
= Zij
(τk− ) + µ̃ij (τk− ) − µ̃ij (τk+ ) pij (τk )τk0 = Zij
(τk− )
(88)
In the above equation, µ̃ij (τk+ ) = µij because Xi (τk+ ) > 0.
k)
Also, µij pij (τk ) = σi (τk ) and µ̃ij (τk− ) = pσiji (τ
(τk ) results
in µ̃ij (τk+ ) = µij . For Zil (t), l 6= j , agent l cannot be
connected to target i at τk so we have, Zil0 (τk+ ) = Zil0 (τk− )
d 0
and dt
Zil (t) = 0 in [τk , τk+1 ). For Zrl (t) ,r 6= i and l 6= j
using the assumption that two events occur at the same time

w.p. 0, the dynamics of Zrl (t) are not affected at τk , hence
d 0
Zrl (t) from (83) for i and j replaced by r and l.
we get dt
0
3. Event ζij
: This event causes a transition from Zij (t) > 0
for t < τk to Zij (t) = 0 for t ≥ τk . The switching function
∂gk
is gk (Θ, X) = Zij (t) so ∂Z
= 1. From (39):
 ∂g
−1ij
∂gk 0 − 
k
τk 0 = −
fk (τk− )
Z (τ )
gk0 +
∂Zij
∂Zij ij k
(89)
0
0
Zij
(τk− )
Zij
(τk− )
=−
=
βij pBj (τk )
µ̃ij (τk− )pij (τk− ) − βij pBj (τk− )
Since Zij (t) is being emptied at τk , by the assumption that
agents can not be in range with the base and targets at the
same time, we have pij (τk )h= 0. Then from (41):
i
0
0
Zij
(τk+ ) = Zij
(τk− ) + − βij pBj (τk ) − 0 τk 0
h
i Z 0 (τ − )
(90)
ij k
0
=0
= Zij
(τk− ) − βij pBj (τk )
βij pBj (τk )
Since Żij (t) = 0 in [τk , τk+1 ):
∂ Żij (t) 0
d 0
0
(t) = 0
(91)
Zij (t) =
Z (t) + Żij
dt
∂Zij (t) ij
For Zrl (t), r 6= i or l 6= j, the dynamics in (7) are not
affected at τk , hence:
0
0
Zrl
(τk+ ) = Zrl
(τk− )
(92)
d 0
if Zrl (τk ) > 0, the value for dt Zrl (t) is calculated by (83)
with r and l replacing i and j respectively. If Zrl (τk ) = 0
d 0
then dt
Zrl (t) = 0.
For Yi (t) we have βi (τk+ ) = 0 since the agent has emptied
its queue, hence:
h
i
Yi0 (τk+ ) = Yi0 (τk− ) + Ẏi (τk− ) − Ẏi (τk+ ) τk0
= Yi0 (τk− ) + [βij pBj (τk ) − 0]

0
Zij
(τk− )
βij pBj (τk )

(93)

0
= Yi0 (τk− ) + Zij
(τk− )
d
In [τk , τk+1 ) we can get dt
Yi0 (t) = 0. For Yr (t), r 6= i the
dynamics of Yr (t) in (8) are not affected by the event at
d
τk hence, Yr0 (τk+ ) and dt
Yr0 (t) in [τk , τk+1 ) are calculated
from (79) and (80) respectively. The dynamics of Xr (t),
r = 1, . . . , M is are not affected at τk since the event at
τk is happening at the base. We have Xr0 (τk+ ) = Xr0 (τk− ).
d
If Xr (τk ) > 0 then we have dt
Xr0 (t) from (78) and if
d
0
Xr (τk ) = 0 then dt Xr (t) = 0 in [τk , τk+1 ).
+
+
4. Event δij
: This event causes a transition from Dij
(t) = 0
+
for t ≤ τk to Dij (t) > 0 for to t > τk . It is the moment
that agent j leaves target i’s range. The switching function
is gk (Θ, X) = dij (t) − rij , from (39):
−1
∂dij 0  ∂dij
τk 0 = −
sj (t)
ṡj (τk )
(94)
∂sj
∂sj
If agent j was connected to target i at τk then by leaving the
target, it is possible that another agent l which is within range
with target i connects to that target. This means Ẋi (τk+ ) =
σi (τk ) − µil pil (τk ) and Ẋi (τk− ) = σi (τk ) − µij pij (τk ), with
pij (τk ) = 0, from (41) we have
Xi0 (τk+ ) = Xi0 (τk− ) − µil pil (τk )τk0
(95)
d
If Xi (τk ) > 0, dt
Xi0 (t) in [τk , τk+1 ) is as in (78) with r
d
replaced by i and if Xi (τk ) = 0 then dt
Xi0 (t) = 0. On the
other hand, if agent j was not connected to target i at τk , we
know that some l 6= j is already connected to target i. This

means agent j leaving target i cannot affect the dynamics
d
Xi0 (t) is
of Xi (t) so we have Xi0 (τk+ ) = Xi0 (τk− ) and dt
calculated from (78) with r replaced by i.
For Xr (t), r 6= i the dynamics in (5) are not affected by the
event at τk hence, we get Xr0 (τk+ ) from (77). If Xr (τk ) > 0
d
Xr0 (t) in [τk , τk+1 ) can be calculated
the time derivative dt
d
from (78) and if Xr (τk ) = 0 then dt
Xr0 (t) = 0.
For Yr (t), r = 1, . . . , , M , the dynamics in (8) are not also
affected by the event at τk hence, we get Yr (τk+ ) from (79)
d
Yr0 (t) is calculated from (80).
and in [τk , τk+1 ) the dt
For Zij (t), the dynamics in (7) are not affect at τk , regardless
of the fact that agent j is connected to target i or not. We have
Żij (τk− ) = µ̃ij (τk )pij (τk ) with pij (τk ) = 0 and Żij (τk+ ) =
0, hence from (41):
h
i
0
0
Zij
(τk+ ) = Zij
(τk− ) + Żij (τk− ) − Żij (τk+ ) τk0
(96)
0
0
= Zij
(τk− ) + µ̃ij (τk )pij (τk )τk0 = Zij
(τk− )
d 0
Zij (t) = 0 using (83) knowing
and in [τk , τk+1 ) , we have dt
pij (τk ) = pBj (τk ) = 0. For Zrl (t), r 6= i or l 6= j, the
dynamics of Zrl (t) are not affected at τk hence (92) holds
and in [τk , τk+1 ) again we can use (83) with i and j replaced
by r and l.
+
0
5. Event δij
: This event causes a transition from Dij
(t) > 0
+
for t < τk to Dij (t) = 0 for to t ≥ τk . The event is the
moment that agent j enters target i’s range. The switching
function is gk (Θ, X) = dij (t)−rij . From (39) we can get τk 0
from (94). If no other agent is already connected to target i,
agent j connects to it. Otherwise, if another agent is already
connected to target i, no connection is established. For Xi (t),
the dynamics in (5) are not affected in both cases, hence, (87)
d
Xi0 (t) using
holds. If Xi (t) > 0 in [τk , τk+1 ) we calculate dt
(78) with l being the appropriate connected agent to target i.
d
If Xi (τk− ) = 0, dt
Xi0 (t) = 0. For Xr (t), r 6= i the dynamics
in (5) are not affected by the event at τk . Hence, we get
d
Xr0 (τk+ ) from (77). If Xr (τk ) > 0 we calculate dt
Xr0 (t)
from (78) with i replaced by r and if Xr (τk ) = 0 then
d
0
dt Xr (t) = 0.
For Yr (t), r = 1, . . . , M again the dynamics in (8) are not
affected at τk so both (79) and (80) hold.
For Zij (t), with agent j being connected or not to target i at
τk the dynamics of Zij (t) are unaffected at τk , hence (92)
d 0
holds for i and j and in [τk , τk+1 ) the dt
Zij (t) is calculated
through (83). For Zrl (t), r 6= i or l 6= j the dynamics are
d 0
unaffected (92) holds again. In [τk , τk+1 ), dt
Zrl (t) is given
through (83) with i and j replaced by r and l.
+
6. Event ∆+
j : This event causes a transition from DBj (t) = 0
+
for t ≤ τk to DBj (t) ≥ 0 for t > τk . The switching function
is gk (Θ, X) = dBj (t) − rBj .
 ∂d
−1
∂DBj 0
Bj
τk 0 = −
sj (τk )
ṡj (τk )
(97)
∂sj
∂sj
Similar to the previous event, the dynamics of Xi (t) are
unaffected at τk hence, we have Xi0 (τk+ ) calculated from
d
Xi0 (t) through
(87). If Xi (t) > 0 in [τk , τk+1 ) we calculate dt
−
d
0
(78) and if Xi (τk ) = 0, dt Xi (t) = 0.
For Yr (t), r = 1, . . . , , M , the dynamics of Yr (t) in (8) are
not affected at τk , hence, we get Yr (τk+ ) from (79) and in
d
[τk , τk+1 ), dt
Yr0 (t) is calculated from (80).

For Zij (t), Using the fact that agent j can only be connected
to one target or the base, we have Żij (τk− ) = βij (τk )pBj (τk )
with pBj (τk ) = 0 and Żij (τk+ ) = 0, hence (92) holds with i
and j replacing r and l. In [τk , τk+1 ) from (40):
∂ Żij (t) 0
d 0
0
Z (t) =
Z (t) + Żij
(t)
dt ij
∂Zij (t) ij
(98)
0
0
(t)
= Żij
(t) = −βij PBj
As for Zrl (t), r 6= i or l 6= j the dynamics are unaffected so
d 0
(92) holds. In [τk , τk+1 ) we can calculate dt
Zrl (t) through
(83) with j replacing l.
+
7. Event ∆0j : This event causes a transition from DBj
(t) > 0
+
for t < τk to DBj (t) = 0 for t ≥ τk . The switching function
is gk (Θ, X) = dBj (t) − rBj . Using (39) we can get τk 0 from
(97). Similar with the previous event we have Xi0 (τk+ ) from
d
Xi0 (t) from (78) and if
(87). If Xi (t) > 0 we can get dt
−
d
0
Xi (τk ) = 0 then dt Xi (t) = 0.
For Yr (t), r = 1, . . . , , M , we again follow the previous
event analysis so (79) and (80) hold.
For Zij (t), the analysis is similar to event ∆+
j so we can
+
d 0
0
calculate Zij (τk ) and dt Zij (t) in [τk , τk+1 ) from (88) and
(83) respectively. Also for Zrl (t), r 6= i or l 6= j, (92)
holds with same reasoning as previous event. In [τk , τk+1 )
d 0
we calculate dt
Zrl (t) from (83).
A PPENDIX V
O BJECTIVE FUNCTION GRADIENT
From (43) we have:
K

∇L(Θ, T ; X(Θ; 0))) =

1 hX
T

k=0

Z

τk+1



q∇L1 (Θ, t)

τk

 i (99)
− (1 − q)∇L2 (Θ, t) + ∇L3 (Θ, t) + ∇L4 (Θ, t) dt
+ ∇Lf (Θ, T )
We calculate each term separately:
M
1 X
αi Xi0 (t)
∇L1 (Θ, t) =
MX i=1
M
1 X
∇L2 (Θ, t) =
αi Yi0 (t)
MY i=1

(100)
(101)

∇L3 (Θ, t) =
!
M
M
M
Y
X
Y
1
+
+0
+
+0
+
d (t)
dij (t) + dBj (t)
dlj (t)
dij (t)
MI Ij (t) Bj
i=1
l=1

i=1,i6=l

(102)
N hZ
X



R(w, t) + RBj (w, t) Pj0 (w, t)dw

1
∇L4 (Θ, t) =
MR j=1 S

Z 
i
0
+
R0 (w, t) + RBj
(w, t) Pj (w, t)dw
S


Z 
N
1 Xh
=
R(w, t) + RBj (w, t) 2hsj (t) − w, s0j (t)idw
MR j=1 S
PM
Z X
M
0 
i
αi Xi0 (t)
i=1 αi Zij
+
Pj (w, t)dw
+
+
+
dB (w)
di (w)
S
i=1

(103)

M
1 X
0
∇Lf (Θ, T ) =
αi Zij
(T )
MZ i=1

(104)

R EFERENCES
[1] P. Corke, T. Wark, R. Jurdak, W. Hu, P. Valencia, and D. Moore,
“Environmental wireless sensor networks,” in Proc. of the IEEE,
vol. 98, pp. 1903–1917, 2010.
[2] R. N. Smith, M. Schwager, S. L. Smith, D. Rus, and G. S. Sukhatme,
“Persistent ocean monitoring with underwater gliders: Towards accurate reconstruction of dynamic ocean processes,” in Proc. - IEEE Int.
Conf. on Robotics and Automation, pp. 1517–1524, 2011.
[3] Z. Tang and U. Özgüner, “Motion planning for multitarget surveillance
with mobile sensor agents,” IEEE Trans. on Robotics, vol. 21, pp. 898–
908, 2005.
[4] M. Zhong and C. G. Cassandras, “Distributed coverage contorol
and data collection with mobile sensor networks,” IEEE Trans. on
Automatic Cont.,, vol. 56, no. 10, pp. 2445–2455, 2011.
[5] K. Chakrabarty, S. S. Iyengar, H. Qi, and E. Cho, “Grid coverage for
surveillance and target location in distributed sensor networks,” IEEE
Trans. on Computers, vol. 51, no. 12, pp. 1448–1453, 2002.
[6] M. Cardei, M. T. Thai, Y. Li, and W. Wu, “Energy-efficient target
coverage in wireless sensor networks,” 24th Annual INFOCOM 2005.,
pp. 1976–1984, 2005.
[7] S. Alamdari, E. Fata, and S. L. Smith, “Persistent monitoring in
discrete environments: Minimizing the maximum weighted latency
between observations,” The Int. J. of Robotics Research, 2013.
[8] C. G. Cassandras, X. Lin, and X. Ding, “An optimal control approach
to the multi-agent persistent monitoring problem,” IEEE Trans. on Aut.
Cont., vol. 58, pp. 947–961, April 2013.
[9] D. Panagou, M. Turpin, and V. Kumar, “Decentralized goal assignment
and trajectory generation in multi-robot networks: A multiple lyapunov
functions approach,” in Robotics and Automation (ICRA), 2014 IEEE
Int. Conf. on, pp. 6757–6762, May 2014.
[10] A. T. Klesh, P. T. Kabamba, and A. R. Girard, “Path planning
for cooperative time-optimal information collection,” Proc. of the
American Cont. Conf., pp. 1991–1996, 2008.
[11] J. L. Ny, M. a. Dahleh, E. Feron, and E. Frazzoli, “Continuous path
planning for a data harvesting mobile server,” Proc. of the IEEE Conf.
on Decision and Cont., pp. 1489–1494, 2008.
[12] R. Moazzez-Estanjini and I. C. Paschalidis, “On delay-minimized data
harvesting with mobile elements in wireless sensor networks,” Ad Hoc
Networks, vol. 10, pp. 1191–1203, 2012.
[13] A. Blum, P. Chalasani, D. Coppersmith, B. Pulleyblank, P. Raghavan,
and M. Sudan, “The minimum latency problem,” in Proceedings of
the twenty-sixth annual ACM symposium on Theory of computing,
pp. 163–171, ACM, 1994.
[14] C. G. Cassandras and W. Li, “Sensor networks and cooperative
control,” European Journal of Control, vol. 11, no. 4-5, pp. 436–463,
2005.
[15] J. K. Hart and K. Martinez, “Environmental sensor networks: A revolution in the earth system science?,” Earth-Science Reviews, vol. 78,
no. 34, pp. 177–191, 2006.
[16] O. Tekdas, V. Isler, J. H. Lim, and A. Terzis, “Using mobile robots
to harvest data from sensor fields,” Wireless Communications, IEEE,
vol. 16, no. 1, pp. 22–28, 2009.
[17] W. Wei, V. Srinivasan, and C. Kee-Chaing, “Extending the lifetime
of wireless sensor networks through mobile relays,” Networking,
IEEE/ACM Transactions on, vol. 16, no. 5, pp. 1108–1120, 2008.
[18] W. Zhao, M. Ammar, and E. Zegura, “Controlling the mobility
of multiple data transport ferries in a delay-tolerant network,” in
INFOCOM 2005. 24th Annual Joint Conference of the IEEE Computer
and Communications Societies. Proceedings IEEE, vol. 2, pp. 1407 –
1418 vol. 2, March 2005.
[19] A. Pandya, A. Kansal, and G. Pottie, “Goodput and delay in networks
with controlled mobility,” in Aerospace Conference, 2008 IEEE, pp. 1–
8, 2008.
[20] K. Akkaya and M. Younis, “A survey on routing protocols for wireless
sensor networks,” Ad Hoc Networks, vol. 3, pp. 325–349, 2005.
[21] M. Liu, Y. Yang, and Z. Qin, “A survey of routing protocols and
simulations in delay-tolerant networks,” Lecture Notes in Computer
Science, vol. 6843 LNCS, pp. 243–253, 2011.
[22] C. Chang, G. Yu, T. Wang, and C. Lin, “Path Construction and Visit
Scheduling for Targets using Data Mules,” IEEE Trans. on Sys, Man,
and Cybernetics: Systems, vol. 44, no. 10, pp. 1289–1300, 2014.
[23] X. Lin and C. G. Cassandras, “An optimal control approach to the
multi-agent persistent monitoring problem in two-dimensional spaces,”
IEEE Transactions on Automatic Control, vol. 60, pp. 1659–1664, June
2015.

[24] C. G. Cassandras, Y. Wardi, C. G. Panayiotou, and C. Yao, “Perturbation analysis and optimization of stochastic hybrid systems,” European
Journal of Cont., vol. 16, no. 6, pp. 642 – 661, 2010.
[25] Y. Khazaeni and C. G. Cassandras, “Event excitation for event-driven
control and optimization of multi-agent systems,” in 13th International
Workshop on Discrete Event Systems, 2016.
[26] A. E. Bryson and Y. C. Ho, Applied optimal control: optimization,
estimation and control. CRC Press, 1975.
[27] H. Kushner and G. Yin, Stochastic Approximation and Recursive
Algorithms and Applications. Springer, 2003.
[28] C. T. Zahn and R. Z. Roskies, “Fourier Descriptors for Plane Closed
Curves,” IEEE Trans. on Computers, vol. C-21, no. 3, 1972.

