Event excitation for event-driven control and optimization of
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:1604.00691v1 [math.OC] 3 Apr 2016

yas@bu.edu,cgc@bu.edu

Abstract— We consider event-driven methods in a general
framework for the control and optimization of multi-agent
systems, viewing them as stochastic hybrid systems. Such
systems often have feasible realizations in which the events
needed to excite an on-line event-driven controller cannot occur,
rendering the use of such controllers ineffective. We show that
this commonly happens in environments which contain discrete
points of interest which the agents must visit. To address this
problem in event-driven gradient-based optimization problems,
we propose a new metric for the objective function which
creates a potential field guaranteeing that gradient values
are non-zero when no events are present and which results
in eventual event excitation. We apply this approach to the
class of cooperative multi-agent data collection problems using the event-driven Infinitesimal Perturbation Analysis (IPA)
methodology and include numerical examples illustrating its
effectiveness.

I. I NTRODUCTION
The modeling and analysis of dynamic systems has historically been founded on the time-driven paradigm provided
by a theoretical framework based on differential (or difference) equations: we postulate the existence of an underlying
“clock” and with every “clock tick” a state update is performed which synchronizes all components of the system. As
systems have become increasingly networked, wireless, and
distributed, the universal value of this paradigm has come
to question, since it may not be feasible to guarantee the
synchronization of all components of a distributed system,
nor is it efficient to trigger actions with every time step when
such actions may be unnecessary. The event-driven paradigm
offers an alternative to the modeling, control, communication, and optimization of dynamic systems. The main idea
in event-driven methods is that actions affecting the system
state need not be taken at each clock tick. Instead, one can
identify appropriate events that trigger control actions. This
approach includes the traditional time-driven view if a clocktick is considered a system “event”. Defining the right events
is a crucial modeling step and has to be carried out with a
good understanding of the system dynamics.
The importance of event-driven behavior in dynamic systems was recognized in the development of Discrete Event
The authors work is supported in part by NSF under grants CNS1239021, ECCS-1509084, and IIP-1430145, by AFOSR under grant
FA9550-15-1-0471, by ONR under grant N00014-09-1-1051, and by the
Cyprus Research Promotion Foundation under Grant New Infrastructure
Project/Strategic/0308/26.

Systems (DES) and later Hybrid Systems (HS) [1]. More
recently there have been significant advances in applying
event-driven methods (also referred to as “event-based” and
“event-triggered”) to classical feedback control systems; e.g.,
see [2], [3], [4], as well as [5] and [6] and references
therein. Event-driven approaches are also attractive in receding horizon control, where it is computationally inefficient
to re-evaluate an optimal control value over small time
increments as opposed to event occurrences defining appropriate planning horizons for the controller (e.g., see [7]).
In distributed networked systems, event-driven mechanisms
have the advantage of significantly reducing communication
among networked components which cooperate to optimize
a given objective. Maintaining such cooperation normally
requires frequent communication among them; it was shown
in [8] that we can limit ourselves to event-driven communication and still achieve optimization objectives while
drastically reducing communication costs (hence, prolonging
the lifetime of a wireless network), even when delays are
present (as long as they are bounded).
Clearly, the premise of these methods 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: it is possible that under some control no such events
are excited, in which case the controller may be useless.
In such cases, one can resort to artificial “timeout events”
so as to eventually take actions, but this is obviously inefficient. Moreover, in event-driven optimization mechanisms
this problem results in very slow convergence to an optimum
or in an algorithm failing to generate any improvement in the
decision variables being updated.
In this work, we address this issue of event excitation in
the context of multi-agent systems. In this case, the events
required are often defined by an agent “visiting” a region
or a single point in a mission space S ⊂ R2 . Clearly, it is
possible that such events never occur for a large number
of feasible agent trajectories. This is a serious problem
in trajectory planning and optimization tasks which are
common in multi-agent systems seeking to optimize different
objectives associated with these tasks, including coverage,
persistent monitoring or formation control [9], [10], [11],
[12], [13], [14], [15], [16]. At the heart of this problem is
the fact that objective functions for such tasks rely on a non-

zero reward (or cost) metric associated with a subset S + ⊂ S
of points, while all other points in S have a reward (or cost)
which is zero since they are not “points of interest” in the
mission space. We propose a novel metric which allows all
points in S to acquire generally non-zero reward (or cost),
thus ensuring that all events are ultimately excited. This leads
to a new method allowing us to apply event-based control
and optimization to a large class of multi-agent problems. We
will illustrate the use of this method by considering a general trajectory optimization problem in which Infinitesimal
Perturbation Analysis (IPA) [1] is used as an event-driven
gradient estimation method to seek optimal trajectories for
a class of multi-agent problems where the agents must
cooperatively visit a set of target points to collect associated
rewards (e.g., to collect data that are buffered at these points.)
This defines a family within the class of Traveling Salesman
Problems (TSPs) [17] for which most solutions are based on
techniques typically seeking a shortest path in the underlying
graph. These methods have several drawbacks: (i) they are
generally combinatorially complex, (ii) they treat agents as
particles (hence, not accounting for limitations in motion
dynamics which should not, for instance, allow an agent to
form a trajectory consisting of straight lines), and (iii) they
become computationally infeasible as on-line methods in the
presence of stochastic effects such as random target rewards
or failing agents. As an alternative we seek solutions in terms
of parameterized agent trajectories which can be adjusted on
line as a result of random effects and which are scalable,
hence computationally efficient, especially in problems with
large numbers of targets and/or agents. This approach was
successfully used in [18], [19].
In section II we present the general framework for multiagent problems and address the event excitation issue. In
section III we overview the event-driven IPA methodology
and how it is applied to a general hybrid system optimization
problem. In section IV we introduce a data collection problem as an application of the general framework introduced
in section II and will show simulation results of applying the
new methodology to this example in section V.
II. E VENT-D RIVEN O PTIMIZATION IN M ULTI -AGENT
S YSTEMS
Multi-agent systems are commonly modeled as hybrid
systems with time-driven dynamics describing the motion
of the agents or the evolution of physical processes in a
given environment, while event-driven behavior characterizes
events that may occur randomly (e.g., an agent failure) or
in accordance to control policies (e.g., an agent stopping
to sense the environment or to change directions). In some
cases, the solution of a multi-agent dynamic optimization
problem is reduced to a policy that is naturally parametric. As
such, a multi-agent system can be studied with parameterized
controllers aiming to meet certain specifications or to optimize a given performance metric. Moreover, in cases where
such a dynamic optimization problem cannot be shown to be
reduced to a parametric policy, using such a policy is still
near-optimal or at least offers an alternative.

Fig. 1.

Multi-agent system in a dynamic setting, blue areas are obstacles

In order to build a general framework for multi-agent
optimization problems, assuming S as the mission space,
we introduce the function R(w) : S → R as a “property” of
point w ∈ S. For instance, R(w) could be a weight that gives
relative importance to one point in S compared to another.
Setting R(w) > 0 for only a finite number of points implies
that we limit ourselves to a finite set of points of interest
while the rest of S has no significant value.
Assuming F to be the set of all feasible agent states, We
define P (w, s) : S × F → R to capture the cost/reward
resulting from how agents with state s ∈ F interact with
w ∈ S. For instance, in coverage problems if an “event”
occurs at w, then P (w, s) is the probability of agents jointly
detecting such events based on the relative distance of each
agent from w.
In general settings, the objective is to find the best state
vector s1 , · · · , sN so that N agents achieve a maximal
reward (minimal cost) from interacting with the mission
space S:
Z
min J =
P (w, s)R(w)dw
(1)
s∈F
S
This static problem can be extended to a dynamic version
where the agents determine optimal trajectories si (t), t ∈
[0, T ], rather than static states:
Z TZ
min J =
P (w, s(u(t)))R(w, t)dwdt
(2)
u(t)∈U

0

S

subject to motion dynamics:
ṡj (t) = fj (sj , uj , t), j = 1, · · · , N
(3)
In Fig. 1, such a dynamic multi agent system is illustrated.
As an example, consensus problems are just a special case
of (1). Suppose that we consider a finite set of points w ∈ S
which coincide with the agents states s1 , ..., sN (which are
not necessarily their locations). Then we can set P (w, s) =
ksi − sj k2 and, therefore, replace the integral in (1) by a
sum. In this case, R(w) = Ri is just the weight that an
agent carries in the consensus algorithm. An optimum occurs
when ksi − sj k2 = 0 for all i, j, i.e., all agents “agree”
and consensus is reached. This is a special case because
of the simplicity in P (w, s) making the problem convex so
that a global optimum can be achieved, in contrast to most
problems we are interested in.
As for the formulation in (2), consider a trajectory planning problem where N mobile agents are tasked to visit
M stationary targets in the mission space S. Target behavior is described through state variables xi (t) which may
model reward functions, the amount of data present at i, or
other problem-dependent target properties. More formally, let
(Ω, F, P) be an appropriately defined probability space and

ω ∈ Ω a realization of the system where target dynamics are
subject to random effects:
ẋi (t) = gi (xi (t), ω)
(4)
gi (·) is as such that xi (t) is monotonically increasing by t
and it resets to zero each time a target is completely emptied
by an agent. In the context of (2), we assume the M targets
are located at pointswi , i = 1, · · · , M and define
R(xi (t), w) if w ∈ C(wi )
R(w, t) =
(5)
0
otherwise
to be the value of point w, where C(wi ) is a compact 2manifold in R2 containing wi which can be considered to be
a region defined by the sensing range of that target relative
to agents (e.g., a disk centered at wi ). Note that R(w, t) is
also a random variable defined on the same probability space
above. Given that only points w ∈ C(wi ) have value for the
agents, there is an infinite number of points w ∈
/ C(wi ) such
that R(w, t) = 0 provided the following condition holds:
Condition 1: If ∃i such that w ∈ C(wi ) then w ∈
/ C(wj )
holds ∀j 6= i.
This condition is to ensure that two targets do not share
any point w in their respective sensing ranges. Also it ensures
that the set {C(wi ) | i = 1 : · · · , M } does not create a
compact partitioning of the mission space and there exist
points w which do not belong to any of the C(wi ).
Viewed as a stochastic hybrid system, we may define different modes depending on the states of agents or targets and
events that cause transitions between these modes. Relative
to a target i, any agent has at least two modes: being at a
point w ∈ C(wi ), i.e., visiting this target or not visiting it.
Within each mode, agent j’s dynamics, dictated by (3), and
target i’s dynamics in (4) may vary. Accordingly, there are
0
at least two types of events in such a system: (i) δij
events
+
occur when agent j initiates a visit at target i, and (ii) δij
events occur when agent j ends a visit at target i. Additional
event types may be included depending on the specifics of
a problem, e.g., mode switches in the target dynamics or
agents encountering obstacles.
An example is shown in Fig. 2, where target sensing ranges
are shown with green circles and agent trajectories are shown
in dashed lines starting at a base shown by a red triangle.
In the blue trajectory, agent 1 moves along the trajectory
that passes through points A → B → C → D. It is easy
to see that when passing through points A and C we have
0
δi1
and δi00 1 events, while passing through B and D we
+
have δi1
and δi+0 1 events. The red trajectory is an example
where none of the events is excited. Suppose we consider
an on-line trajectory adjustment process in which the agent
improves its trajectory based on its performance measured
through (5). In this case, R(w, t) = 0 over all t, as long
as the agent keeps using the red trajectory, i.e., no event
ever occurs. Therefore, if an event-driven approach is used
to control the trajectory adjustment process, no action is ever
triggered and the approach is ineffective. In contrast, in the
blue trajectory the controller can extract useful information
from every observed event; such information (e.g., a gradient
of J with respect to controllable parameters as described in
the next section) can be used to adjust the current trajectory

Fig. 2.

Sample trajectories

so as to improve the objective function J in (1) or (2).
Therefore, if we are to build an optimization framework
for this class of stochastic hybrid systems to allow the application of event-driven methods by calculating a performance
measure gradient, then a fundamental property required is the
occurrence of at least some events in a sample realization. In
particular, the IPA method [20] is based on a single sample
realization of the system over which events are observed
along with their occurrence times and associated system
states. Suppose that the trajectories can be controlled through
a set of parameters forming a vector θ. Then, IPA provides
an unbiased estimate of the gradient of a performance metric
J(θ) with respect to θ. This gradient is then used to improve
the trajectory and ultimately seek an optimal one when
appropriate conditions hold.
As in the example of Fig. 2, it is possible to encounter
trajectory realizations where no events occur in the system.
In the above example, this can easily happen if the trajectory
does not pass through any target. The existence of such
undesirable trajectories is the direct consequence of Condition 1. This lack of event excitation results in event-based
controllers being unsuitable.
New Metric: In order to overcome this issue we propose
a new definition for R(w, t) in (5) as follows:
M
X
R(w, t) =
hi (xi (t), di (w))
(6)
i=1

where w ∈ S, hi (·) is a function of the target’s state xi (t)
and di (w) = kwi −wk. Note that, if hi (·) is properly defined,
(6) yields R(w, t) > 0 at all points.
While the exact form of hi (·) depends on the problem, we
impose the condition that hi (·) is monotonically decreasing
in di (w). We can think of hi (·) as a value function associated
with point wi . Using the definition of R(w, t), this value
is spread out over all points w ∈ S rather than being
concentrated at the single point wi . This creates a continuous
potential field for the agents leading to a non-zero gradient of
the performance measure even when the trajectories do not
excite any events. This non-zero gradient will then induce
trajectory adjustments that naturally bring them toward ones
with observable events.
Finally, recalling the definition in (2), we also define:
N
X
P (w, s) =
ksj (t) − wk2
(7)
j=1

the total quadratic travel cost for agents to visit point w.
In Section IV, we will show how to apply R(w, t) and
P (w, s) defined as above in order to determine optimal agent
trajectories for a class of multi-agent problems of the form

(2). First, however, we review in the next section the eventdriven IPA calculus which allows us to estimate performance
gradients with respect to controllable parameters.
III. E VENT-D RIVEN IPA C ALCULUS
Let us fix a particular value of the parameter θ ∈ Θ and
study a resulting sample path of a general SHS. Over such a
sample path, let τk (θ), k = 1, 2, · · · denote the occurrence
times of the discrete events in increasing order, and define
τ0 (θ) = 0 for convenience. We will use the notation τk
instead of τk (θ) when no confusion arises. The continuous
state is also generally a function of θ, as well as of t, and is
thus denoted by x(θ, t). Over an interval [τk (θ), τk+1 (θ)),
the system is at some mode during which the time-driven
state satisfies ẋ = fk (x, θ, t), in which x is any of the
continuous state variables of the system and ẋ denotes ∂x
∂t .
Note that we suppress the dependence of fk on the inputs
u ∈ U and d ∈ D and stress instead its dependence on
the parameter θ which may generally affect either u or d
or both. The purpose of perturbation analysis is to study
how changes in θ influence the state x(θ, t) and the event
times τk (θ) and, ultimately, how they influence interesting
performance metrics that are generally expressed in terms of
these variables.
An event occurring at time τk+1 (θ) triggers a change
in the mode of the system, which may also result in new
dynamics represented by fk+1 . The event times τk (θ) play
an important role in defining the interactions between the
time-driven and event-driven dynamics of the system.
Following the framework in [20], consider a general
performance function J of the control parameter θ:
J(θ; x(θ, 0), T ) = E[L(θ; x(θ, 0), T )]
(8)
where L(θ; x(θ, 0), T ) is a sample function of interest evaluated in the interval [0, T ] with initial conditions x(θ, 0).
For simplicity, we write J(θ) and L(θ). Suppose that there
are K events, with occurrence times generally dependent on
θ, during the time interval [0, T ] and define τ0 = 0 and
τN +1 = T . Let Lk : Rn × Θ × R+ → R be a function and
define L(θ) by
K Z τk+1
X
L(θ) =
Lk (x, θ, t)dt
(9)
k=0

τk

where we reiterate that x = x(θ, t) is a function of θ and
t. We also point out that the restriction of the definition of
J(θ) to a finite horizon T which is independent of θ is made
merely for the sake of simplicity. Returning to the stochastic
setting, the ultimate goal of the iterative process shown is to
maximize Eω [L(θ, ω)], where we use ω to emphasize dependence on a sample path ω of a SHS (clearly, this is reduced to
L(θ) in the deterministic case). Achieving such optimality is
possible under standard ergodicity conditions imposed on the
underlying stochastic processes, as well as the assumption
that a single global optimum exists; otherwise, the gradientbased approach is simply continuously attempting to improve
the observed performance L(θ, ω). Thus, we are interested
in estimating the gradient
dEω [L(θ, ω)]
dJ(θ)
=
(10)
dθ
dθ

based on directly observed data. We
by evaluating dL(θ,ω)
dθ
obtain θ ∗ by optimizing J(θ) through an iterative scheme
of the form
θn+1 = θn − ηn Hn (θn ; x(θ, 0), T, ωn ), n = 0, 1, · · · (11)
where ηn is a step size sequence and Hn (θn ; x(θ, 0), T, ωn )
at θ = θn . In using IPA,
is the estimate of dJ(θ)
dθ
,
Hn (θn ; x(θ, 0), T, ωn ) is the sample derivative dL(θ,ω)
dθ
which is an unbiased estimate of dJ(θ)
if
the
condition
dθ
(dropping the symbol ω for simplicity)
 dL(θ)  dE[L(θ)]
dJ(θ)
E
=
=
(12)
dθ
dθ
dθ
is satisfied, which turns out to be the case under mild
technical conditions. The conditions under which algorithms
of the form (11) converge are well-known (e.g., see [21]).
Moreover, in addition to being unbiased, it can be shown that
such gradient estimates are independent of the probability
laws of the stochastic processes involved and require minimal
information from the observed sample path. The process
is based on analyzing
through which IPA evaluates dL(θ)
dθ
how changes in θ influence the state x(θ, t) and the event
times τk (θ). In turn, this provides information on how L(θ)
is affected, because it is generally expressed in terms of these
variables. Given θ = [θ1 , ..., θl ]T , we use the Jacobian matrix
notation:
∂x(θ, t)
∂τk (θ)
x0 (θ, t) =
, τk 0 =
, k = 1, · · · , K (13)
∂θ
∂θ
for all state and event time derivatives. For simplicity of
notation, we omit θ from the arguments of the functions
above unless it is essential to stress this dependence. It is
shown in [20] that x0 (t) satisfies:
dx0 (t)
∂fk (t) 0
∂fk (t)
=
x (t) +
(14)
dt
∂x
∂θ
for t ∈ [τk (θ), τk+1| (θ)) with boundary condition
x0 (τk+ ) = x0 (τk− ) + [fk−1 (τk− ) − fk (τk+ )]τk0
(15)
for k = 0, · · · , K. We note that whereas x(t) is often
continuous in t, x0 (t) may be discontinuous in t at the event
times τk ; hence, the left and right limits above are generally
different. If x(t) is not continuous in t at t = τk (θ), the value
of x(τk+ ) is determined by the reset function r(q, q 0 , x, ν, δ)
and
dr(q, q 0 , x, ν, δ)
x0 (τk+ ) =
(16)
dθ
0 +
Furthermore, once the initial condition x (τk ) is given,
the linearized state trajectory x0 (t) can be computed in the
interval t ∈ [τk (θ), τk+1 (θ)) by solving (14) to obtain
R t ∂fk (u)
h Z t ∂f (v) R t ∂fk (u)
i
du
du
−
k
0
∂x
τk
x (t) = e
e τk ∂x
dv + ξk
∂θ
τk
(17)
with the constant ξk determined from x0 (τk+ ). In order to
complete the evaluation of x0 (τk+ ) we need to also determine
τk0 . If the event at τk (θ) is exogenous τk0 = 0 and if the event
at τk (θ) is endogenous:
h ∂g
i ∂g
∂gk 0 − 
k
k
τk0 = −
fk (τk− )
+
x (τk )
(18)
∂x
∂θ
∂x
∂gk
+
where gk (x, θ) = 0 and it is defined as long as ∂x fk (τk ) 6=
0 (details may be found in [20].)
The derivative evaluation process involves using the IPA
calculus in order to evaluate the IPA derivative dL
dθ . This is

accomplished by taking derivatives in (9) with respect to θ:
Z τk+1
K
dL(θ) X d
Lk (x, θ, t)dt
(19)
=
dθ
dθ τk
k=0
Applying the Leibnitz rule, we obtain, for every k =
0, · · · , ZK,
τk+1
d
Lk (x, θ, t)dt
dθ τk
Z τk+1 h
∂Lk (x, θ, t) i
∂Lk (x, θ, t) 0
dt
x (t) +
=
∂x
∂θ
τk
0
+ Lk (x(τk+1 ), θ, τk+1 )τk+1
− Lk (x(τk ), θ, τk )τk0
(20)
In summary the three equations (15), (17) and (18) form
the basis of the IPA calculus and allow us to calculate the
final derivative in (20). In the next section IPA is applied to
a data collection problem in a multi-agent system.

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. This means that j in (22) is the index
of the agent that is connected to target i at time t.
The dynamics of xi (t) in (22) results in two new event
types added to what was defined earlier, (i) ξi0 events occur
when xi (t) reaches zero, and (ii) ξi+ events occur when xi (t)
leaves zero.
The performance measure is the total content of data left
at targets at the end of a finite mission time T . Thus, we
define J1 (t) to be the following (recalling that {σi (t)} are
random processes):
M
X
J1 (t) =
αi E[xi (t)]
(23)
i=1

IV. T HE DATA C OLLECTION P ROBLEM
We consider a class of multi-agent problems where the
agents must cooperatively visit a set of target points to collect
associated rewards (e.g., to collect data that are buffered at
these points.). The mission space is S ⊂ R2 . This class of
problems falls within the general formulation introduced in
(2). The state of the system is the position of agent j time
t, sj (t) = [sxj (t), syj (t)] and the state of the target i, xi (t).
The agent’s dynamics (3) follow a single integrator:
ṡxj (t) = uj (t) cos θj (t),
ṡyj (t) = uj (t) sin θj (t) (21)
where uj (t) is the scalar speed of the agent (normalized so
that 0 ≤ uj (t) ≤ 1) and θj (t) is the angle relative to the
positive direction, 0 ≤ θj (t) < 2π. Thus, we assume that
each agent controls its speed and heading.
We assume the state of the target xi (t) represents the
amount of data that is currently available at target i (this can
be modified to different state interpretations). The dynamics
of xi (t) in (4) for this problem are:
0
if xi (t) = 0 and σi (t) ≤ µij p(sj (t), wi )
ẋi (t) =
σi (t) − µij p(sj (t), wi )
otherwise
(22)
i.e., we model the data at the target as satisfying simple flow
dynamics with an exogenous (generally stochastic) inflow
σi (t) and a controllable rate with which an agent empties
the data queue given by µij p(sj (t), wi ). For brevity we set
p(sj (t), wi ) = pij (t) which is the normalized data collection
rate from target i by agent j and µij is a nominal rate
corresponding to target i and agent j.
Assuming M targets are located at wi ∈ S, i = 1, . . . , M,
and have a finite range of ri , then agent j can collect data
from wi only if dij (t) = kwi − sj (t)k ≤ ri . We then
assume that: (A1) pij (t) ∈ [0, 1] is monotonically nonincreasing in the value of dij (t) = kwi − sj (t)k, and (A2)
it satisfies pij (t) = 0 if dij (t) > ri . 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 [22]) 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) and (A4) only one agent at a time is

where αi is a weight factor for target i. 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).
Z
1 T
J1 (t)dt
(24)
P1 :
min J(T ) =
T 0
u(t),θ(t)
where 0 ≤ uj (t) ≤ 1, 0 ≤ θj (t) < 2π, and T is a given
finite mission time. This problem can be readily placed into
the general framework (2). In particular, the right hand side
of (24) is:"
#
Z T XZ
1
αi
E
2 xi (t)dwdt
T
0
C(wi ) πri
i
"Z Z
# (25)
T
X αi 1{w ∈ C(wi )}
1
xi (t)dwdt
= E
T
πri2
0
S i
This is now in the form of the general framework in (2) with
X αi 1{w ∈ C(wi )}
R(w, t) =
xi (t)
(26)
πri2
i
and
P (sj (t), w) = 1
(27)
Recalling the definition in (5), only points within the sensing
range of each target have non-zero values, while all other
point value are zero, which is the case in (26) above. In
addition, (27) simply shows that there is no meaningful
dynamic interaction between an agent and the environment.
Problem P1 is a finite time optimal control problem. In
order to solve this, following previous work in [19] we
proceed with a standard Hamiltonian analysis leading to a
Two Point Boundary Value Problem (TPBVP) [23]. We omit
this, since the details are the same as the analysis in [19]. The
main result of the Hamiltonian analysis is that the optimal
speed is always the maximum value, i.e.,
u∗j (t) = 1
(28)
Hence, we only need to calculate the optimal θj (t). This
TPBVP is computationally expensive and easily becomes
intractable when problem size grows. The ultimate solution
of the TPBVP is a set of agent trajectories that can be
put in a parametric form defined by a parameter vector

θ and then optimized over θ. If the parametric trajectory
family is broad enough, we can recover the true optimal
trajectories; otherwise, we can approximate them within
some acceptable accuracy. Moreover, adopting a parametric
family of trajectories and seeking an optimal one within it has
additional benefits: it allows trajectories to be periodic, often
a desirable property, and it allows one to restrict solutions to
trajectories with desired features that the true optimal may
not have, e.g., smoothness properties to achieve physically
feasible agent motion.
Parameterizing the trajectories and using gradient based
optimization methods, in light of the discussions from the
previous sections, enables us to make use of Infinitesimal
Perturbation Analysis (IPA) [20] to carry out the trajectory
optimization process. We represent each agent’s trajectory
through general parametric equations
sxj (t) = fx (θj , ρj (t)),
syj (t) = fy (θj , ρj (t))
(29)
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 trajectory. Let θ =
[θ1 , . . . , θN ]. We now revisit problem P1 in (24):
Z
1 T
(30)
J1 (θ, t)dt
min J(θ, T ) =
θ∈Θ
T 0
and will bring in the equations that were introduced in the
previous section in order to calculate an estimate of dJ(θ)
dθ
as in (10). For this problem due to the continuity of xi (t)
the last two terms in (20) vanish. From (23) we have:
Z τk+1 X
Z τk+1 X
M
M
d
αi xi (θ, t)dt =
αi x0i (θ, t)dt (31)
dθ τk
τk
i=1
i=1
In summary, the evaluation of (31) requires the state
derivatives x0i (t) explicitly and s0j (t) implicitly, (dropping the
dependence on θ for brevity). The latter are easily obtained
for any specific choice of f and g in (29). The former require
a rather laborious use of (15),(17),(18) which, reduces to a
simple set of state derivative dynamics as shown next.
Proposition 1. After an event occurrence at t = τk ,
the state derivatives x0i (τk+ ) with respect to the controllable
parameter θ satisfy the following:

if e(τk ) = ξi0
 0
0
−
+
0 +
xi (τk ) =
x0 (τ ) − µil (t)pil (τk )τk if e(τk ) = δij
 i0 k−
xi (τk )
otherwise
where l 6= j with pil (τk ) > 0 if such l exists and
0

τk =

∂dij (sj ) 0
∂sj sj

∂dij (sj )
∂sj ṡj (τk )

−1

.
Proof: The proof is omitted due to space limitations,
but it is very similar to the proofs of Propositions 1-3 in
[24].
As is obvious from Proposition 1, the evaluation of x0i (t)
is entirely dependent on the occurrence of events ξi0 and
+
+
δij
in a sample realization, i.e., ξi0 and δij
cause jumps in
this derivative which carry useful information. Otherwise,
x0i (τk+ ) = x0i (τk− ) is in effect and these gradients remain
unchanged. However, we can easily have realizations where
0
no events occur in the system (specifically, events of type δij
+
and δij ) if the trajectory of agents in the sample realization
does not pass through any target. This lack of event excitation

results in the algorithm in (11) to stall.
In the next section we overcome the problem of no event
excitation using the definitions in (6) and (7). We accomplish
this by adding a new metric to the objective function that
generates a non-zero sensitivity with respect to θ.
A. Event Excitation
Our goal here is to select a function hi (·) in (6) with the
property of “spreading” the value of xi (t) over all w ∈ S.
We begin by determining the convex hull produced by the
targets, since the trajectories need not go outside this convex
hull. Let T = {w1 , w2 , · · · , wM } be the set of all target
points. Then, the convex hull of these points is
X

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

i

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

i=1

so that R(w, t) can be viewed as a continuous density defined
for all points w ∈ C which results in a total value equivalent
to a weighted sum of the target states xi (t), i = 1, . . . , M . In
order to select an appropriate h(xi (t), di (w)) in (6), we first
define d+
i (w) = max(kw − wi k, ri ) where ri is the target’s
sensing range. We then define:
M
X
αi xi (t)
(34)
R(w, t) =
d+ (w)
i=1 i
Here, we are spreading a target’s reward (numerator) over
all w so as to obtain the “total weighted reward density” at
w. Note that d+
i (w) = max(kw − wi k, ri ) > 0 to ensure
that the target reward remains positive and fixed for points
w ∈ C(wi ). Moreover, following (7),
N
X
P (w, s(t)) =
ksj (t) − wk2
(35)
j=1

Using these definitions we introduce a new objective function
metric which is added to
the objective function in (24):
hZ
i
J2 (t) = E
P (w, s(t))R(w, t)dw
(36)
C

The expectation is a result of P (w, s(t)) and R(w, t) being
random variables defined on the same probability space as
xi (t).
Proposition 2. For R(w, t) in (34), there exist ci > 0,
i = 1, . . . , M , such that:
Z
M
X
R(w, t)dw =
ci xi (t)
(37)
C

i=1

Proof: We have
Z
Z X
M
αi xi (t)
R(w, t) =
dw
+
C
C i=1 di (w)
(38)
Z
M
X
xi (t)
=
αi
dw
+
C di (w)
i=1
R
(t)
We now need to find the value of C dx+i(w)
for each target
i
i. To do this we first look at the case of one target in a 2D
space and for now we assume C is just a disk with radius Λ

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

Λ 
= xi (t) 2π 1 + log( )
ri
In our case C is the convex
hull of all targets. We will use
the
idea to calculate the
R xsame
i (t)
dw
for the actual con+
C di (w)
vex hull. We do this for an interior target i.e., a target inside
the convex hull. Extending the
same to targets on the edge is
straightforward. Using the same Fig. 3. One Target R(w, t)
polar coordinate for each θ we Calculation
define Λ(θ) to be the distance of
the target to the edge of C in the direction of θ. (C shown
by a red polygon in Fig. 3).
Z 2π Z Λ
Z
xi (t)
xi (t)
dw
=
drdθ
+
+
0
0 di (r, θ)
C di (w)
Z 2π Z ri
Z 2π Z Λ(θ)
xi (t)
xi (t)
=
drdθ (40)
drdθ +
r
r
i
0
0
0
ri
Z 2π


Λ(θ)
= xi (t) 2π +
log(
)dθ
r
i
0
The second part in (40) has to be calculated knowing Λ(θ)
but since we assumed the target is inside the convex hull we
know Λ(θ) ≥ ri . This means log( Λ(θ)
ri ) > 0 and the xi (t)’s
multiplier is a positive value. We can define ci in (37) as:
Z 2π


Λ(θ)
ci = αi 2π +
)dθ
log(
(41)
r
i
0
The significance of J2 (t) is that it accounts for the
movement of agents through P (w, s(t)) and captures the
target state values through R(w, t). Introducing this term in
the objective function in the following creates a non-zero
gradient even if the agent trajectories are not passing through
any targets. We now combine the two metrics in (24) and
(36) and define problem P2:
Z

1 T
P2 :
min J(T ) =
J1 (t) + J2 (t) dt (42)
T 0
u(t),θ(t)
In this problem, the second term is responsible for adjusting
the trajectories towards the targets by creating a potential
field, while the first term is the original performance metric
which is responsible for adjusting the trajectories so as to
maximize the data collected once an agent is within a target’s
sensing range. It can be easily shown that the results in (28)
hold for problem P2 as well, through the same Hamiltonian
analysis presented in [19]. When sj (t) follows the parametric
functions in (29), the new metric simply becomes a function
of the parameter vector θ and we have:

Z

1 T
J1 (θ, t) + J2 (θ, t) dt
(43)
θ∈Θ
T 0
The new objective function’s derivative follows the same
procedure that was described previously. The first part’s
derivative can be calculated from (31). For the second part
we have:
Z τk+1 Z
min J(θ, T ) =

d
dθ

P (w, θ, t)R(w, θ, t)dw
τk

Z

C
τk+1 Z

=
τk

C

h dP (w, θ, t)
dθ

R(w, θ, t) + P (w, θ, t)

dR(w, θ, t) i
dw
dθ

(44)
In the previous section, we raised the problem of no events
being excited in a sample realization, in which case the total
derivative in (31) is zero and the algorithm in (11) stalls.
Now, looking at (44) we can see that if no events occur
the second part in the integration which involves dR(w,θ,t)
dθ
PM
will be zero, since i=1 x0i (t) = 0 at all t. However, the
first part in the integral does not depend on the events, but
calculates the sensitivity of P (w, s(t)) in (35) with respect
to the parameter θ. Note that the dependence on θ comes
through the parametric description of s(t) through (29). This
term ensures that the algorithm in (11) does not stall and
adjusts trajectories so as to excite the desired events.
V. S IMULATION R ESULTS
We provide some simulation results based on an elliptical parametric description for the trajectories in (29). The
elliptical trajectory formulation is:
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
(45)
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 an agent is moving
with constant speed of 1 on this trajectory, based on (28),
we have ṡxj (t)2 + ṡyj (t)2 = 1, which gives
h
2
ρ̇j (t) = a sin ρj (t) cos φj + bj cos ρj (t) sin φj

2 i− 21
+ a sin ρj (t) sin φj − bj cos ρj (t) cos φj
(46)
The first case we consider is a problem with one agent and
seven targets located on a circle, as shown in Fig. 4. We
consider a deterministic case with σi (t) = 0.5 for all i. The
other problem parameters are T = 50, µij = 100, ri = 0.2
and αi = 1. A target’s sensing range is denoted with solid
black circles with the target location at the center. The blue
polygon indicates the convex hull produced by the targets.
The direction of motion on a trajectory is shown with the
small arrow. Starting with an initial trajectory shown in light
blue, the on-line trajectory optimization process converges
to the trajectory passing through all targets in an efficient
manner (shown in dark solid blue). In contrast, starting with
this trajectory - which does not pass through any targets

using the event-based IPA calculus to estimate the objective
function gradient.
R EFERENCES

Fig. 4.

One agent and seven target scenario

Fig. 5.

Two agent and seven targets scenario

- problem P1 does not converge and the initial trajectory
remains unchanged. At the final trajectory, J1∗ = 0.0859
and J ∗ = 0.2128. Using the obvious shortest path solution,
the actual optimal value for J1 is 0.0739 that results from
moving on the edges of the convex hull (which allows for
shorter agent travel times).
In the second case, 7 targets are randomly distributed
and two agents are cooperatively collecting the data. The
problem parameters are σi = 0.5, µij = 10, ri = 0.5, αi =
1, T = 50. The initial trajectories for both agents are shown
in light green and blue respectively. We can see that both
agent trajectories converge so as to cover all targets, shown
in dark green and blue ellipses. At the final trajectories,
J1∗ = 0.1004 and J ∗ = 0.2979. Note that we may use these
trajectories to initialize the corresponding TPBVP, another
potential benefit of this approach. This is a much slower
process which ultimately converges to J1∗ = 0.0991 and
J ∗ = 0.2776.
VI. C ONCLUSIONS
We have addressed the issue of event excitation in a
class of multi-agent systems with discrete points of interest.
We proposed a new metric for such systems that spreads
the point-wise values throughout the mission space and
generates a potential field. This metric allows us to use eventdriven trajectory optimization for multi-agent systems. The
methodology is applied to a class of data collection problems

[1] C. G. Cassandras and S. Lafortune, Introduction to Discrete Event
Systems. Secaucus, NJ, USA: Springer-Verlag New York, Inc., 2006.
[2] W. Heemels, J. Sandee, and P. P. J. van den Bosch, “Analysis of eventdriven controllers for linear systems,” International journal of control,
vol. 81, no. 4, pp. 571–590, 2008.
[3] A. Anta and P. Tabuada, “To sample or not to sample: Self-triggered
control for nonlinear systems,” Automatic Control, IEEE Trans. on,
vol. 55, pp. 2030–2042, Sept 2010.
[4] S. Trimpe and R. D’Andrea, “Event-based state estimation with
variance-based triggering,” Automatic Control, IEEE Trans. on,
vol. 59, no. 12, pp. 3266–3281, 2014.
[5] M. Miskowicz, Event-Based Control and Signal Processing. CRC
Press, 2015.
[6] C. G. Cassandras, “The event-driven paradigm for control, communication and optimization,” Journal of Control and Decision, vol. 1,
no. 1, pp. 3–17, 2014.
[7] Y. Khazaeni and C. G. Cassandras, “A new event-driven cooperative receding horizon controller for multi-agent systems in uncertain
environments,” In Proceedings of IEEE 53rd Annual Conference on
Decision and Control, pp. 2770–2775, Dec 2014.
[8] M. Zhong and C. G. Cassandras, “Asynchronous distributed optimization with event-driven communication,” Automatic Control, IEEE
Trans. on, vol. 55, no. 12, pp. 2735–2750, 2010.
[9] M. Schwager, D. Rus, and J.-J. Slotine, “Decentralized, adaptive
coverage control for networked robots,” The International Journal of
Robotics Research, vol. 28, no. 3, pp. 357–375, 2009.
[10] 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.
[11] M. Cao, A. Morse, C. Yu, B. Anderson, and S. Dasgupta, “Maintaining
a directed, triangular formation of mobile autonomous agents,” Communications in Information and Systems, vol. 11, no. 1, p. 1, 2011.
[12] K.-K. Oh and H.-S. Ahn, “Formation control and network localization
via orientation alignment,” IEEE Trans. on Automatic Control, vol. 59,
pp. 540–545, Feb 2014.
[13] H. Yamaguchi and T. Arai, “Distributed and autonomous control
method for generating shape of multiple mobile robot group,” in Proc.
of the IEEE International Conf. on Intelligent Robots and Systems,
vol. 2, pp. 800–807 vol.2, Sep 1994.
[14] J. Desai, V. Kumar, and J. Ostrowski, “Control of changes in formation
for a team of mobile robots,” in Proc. of the IEEE International Conf.
on Robotics and Automation, vol. 2, pp. 1556–1561, 1999.
[15] M. Ji and M. B. Egerstedt, “Distributed coordination control of
multi-agent systems while preserving connectedness.,” IEEE Trans.
on Robotics, vol. 23, no. 4, pp. 693–703, 2007.
[16] J. Wang and M. Xin, “Integrated optimal formation control of multiple
unmanned aerial vehicles,” IEEE Trans. on Control Systems Technology, vol. 21, pp. 1731–1744, Sept 2013.
[17] D. L. Applegate, R. E. Bixby, V. Chvatal, and W. J. Cook, The traveling salesman problem: a computational study. Princeton University
Press, 2011.
[18] X. Lin and C. G. Cassandras, “An optimal control approach to the
multi-agent persistent monitoring problem in two-dimensional spaces,”
IEEE Trans. on Automatic Control, vol. 60, pp. 1659–1664, June 2015.
[19] Y. Khazaeni and C. G. Cassandras, “An optimal control approach for
the data harvesting problem,” in 54th IEEE Conf. on Decision and
Cont., pp. 5136–5141, 2015.
[20] 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.
[21] H. Kushner and G. Yin, Stochastic Approximation and Recursive
Algorithms and Applications. Springer, 2003.
[22] 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.
[23] A. E. Bryson and Y. C. Ho, Applied optimal control: optimization,
estimation and control. CRC Press, 1975.
[24] Y. Khazaeni and C. G. Cassandras, “An optimal control approach for
the data harvesting problem,” arXiv:1503.06133.

