1

Random Access Communication for Wireless
Control Systems with Energy Harvesting Sensors

arXiv:1801.10141v1 [math.OC] 30 Jan 2018

Miguel Calvo-Fullana, Carles Antón-Haro, Javier Matamoros, and Alejandro Ribeiro

Abstract—In this paper, we study wireless networked control
systems in which the sensing devices are powered by energy
harvesting. We consider a scenario with multiple plants, where
the sensors communicate their measurements to their respective
controllers over a shared wireless channel. Due to the shared
nature of the medium, sensors transmitting simultaneously can
lead to packet collisions. In order to deal with this, we propose
the use of random access communication policies and, to this end,
we translate the control performance requirements to successful
packet reception probabilities. The optimal scheduling decision
is to transmit with a certain probability, which is adaptive to
plant, channel and battery conditions. Moreover, we provide
a stochastic dual method to compute the optimal scheduling
solution, which is decoupled across sensors, with only some of the
dual variables needed to be shared between nodes. Furthermore,
we also consider asynchronicity in the values of the variables
across sensor nodes and provide theoretical guarantees on the
stability of the control systems under the proposed random access
mechanism. Finally, we provide extensive numerical results that
corroborate our claims.
Index Terms—Energy harvesting, networked control systems,
random access communication.

I. I NTRODUCTION
The rapid pace of development of technologies such as
robotic automation, smart homes, autonomous transportation,
and the internet of things is causing a dramatic increase in the
average number of sensors in modern control systems. Usually,
the previously mentioned technologies rely on networked control systems, and tend to incorporate wireless sensing devices
to perform the monitoring of physical processes. These sensors
might be deployed in large quantities and over large areas,
making the replacement of their batteries a difficult and costly
task. This has led to an increasing interest in alternative ways
of powering wireless devices. An important technology that
has recently emerged as capable of alleviating the limitations
imposed by traditional battery operation is Energy Harvesting
(EH). The use of energy harvesting technologies allows the
devices to obtain energy from their environment (with common
sources being solar, wind or kinetic energy [2]). In turn, this
removes some of the limitations imposed by traditional battery
This work is supported by ARL DCIST CRA W911NF-17-2-0181 and the
Intel Science and Technology Center for Wireless Autonomous Systems.
M. Calvo-Fullana, and A. Ribeiro are with the Department of Electrical and
Systems Engineering, University of Pennsylvania, Philadelphia, PA 19104,
USA (e-mail: cfullana@seas.upenn.edu; aribeiro@seas.upenn.edu).
C. Antón-Haro, and J. Matamoros are with the Centre Tecnològic de Telecomunicacions de Catalunya (CTTC/CERCA), 08860 Castelldefels, Barcelona,
Spain (e-mail: carles.anton@cttc.cat; javier.matamoros@cttc.cat).
This work has been presented in part at the 2017 American Control
Conference (ACC) [1].

operation and grants an increase to the expected lifetime of the
devices.
The study of communication systems powered by energy
harvesting has recently received considerable attention. Current results available in the literature range from throughput
maximization [3]–[6], source-channel coding [7]–[10], estimation [11]–[13], and others (see [14] for a comprehensive
overview). However, in general, limited attention has been
given to the use of energy harvesting technologies in control
applications. Most of the works currently available deal with
the estimation of dynamical systems with sensors powered by
energy harvesting [15]–[18]. Nonetheless, the more explicit
study of closed-loop system stability under energy harvesting
constraints has been less studied, and only for single plant
scenarios [19], [20].
In this paper, we consider the multi-plant problem of
scheduling communication between sensor nodes and their respective controllers. For control systems with classically powered sensor nodes (i.e., not energy harvesting), the scheduling
problem in wireless networked control systems has been previously studied in several forms. The most common approach to
this problem is the design of centralized scheduling policies.
In such setup, in order to avoid packet collisions between the
transmissions of the nodes, there exists an overseeing entity
specifying which sensor is allowed to transmit at a given time
slot. These type of policies might be static [21], [22] or of a
more dynamic nature, where centralized decisions can be taken
based on plant state information [23] or others. Decentralized
policies have received less attention, with the authors in [24]
proposing a random access mechanism that adapts to channel
conditions.
In our case, we study the scenario in which the sensors are
powered by energy harvesting, and we focus on the design of
decentralized scheduling policies. We consider the coexistence
of multiple plants, with sensor nodes transmitting their measurements to their controllers over a shared wireless medium.
Due to this, multiple sensors accessing the medium at the same
time can cause collisions, leading to the unsuccessful reception
of the sensor measurements by the controller. To mitigate
this, we propose to use a random access communication
scheme. First, we abstract the required control performance
into a required successful packet reception probability. Under
this abstraction, a Lyapunov function of each control loop is
required to decrease at a given average rate. Then, we pose
the random access mechanism as an stochastic optimization
problem where the required successful reception probabilities
act as constraints of the optimization problem. Then, the
energy harvesting constraints are introduced into the problem

2

multiple sensors transmitting simultaneously can led to packet
collisions, with the consequential lack of packet delivery. It
is our objective to design transmission policies that adapt to
the wireless medium and the energy harvesting process of the
sensors, and are capable of stabilizing all control loops.
A. Control Model

Fig. 1. System model.

in an average manner and and we modify the formulation
to allow us to ensure time slot to time slot causality in the
stochastic framework. To solve the optimization problem, we
resort to a primal-dual stochastic subgradient method [25].
At a given time slot, the resulting scheduling decision is to
transmit with a certain probability, which is adaptive to the
plant, channel and battery conditions. The resulting policy
requires minimal coordination, with only some of the dual
variables being shared between the sensors. Furthermore, we
consider the possibility of asynchronism between the sensor
nodes (i.e., nodes with outdated dual information) and provide
theoretical guarantees that ensure the stability of all control
loops under these conditions when using our proposed scheme.
Finally, we validate our policy by means of simulations, which
illustrate its ability to adapt to environmental conditions and
satisfy the stability of all control loops.
The rest of the paper is organized as follows. In Section
II we introduce the system model and provide details on its
control, communication and control performance aspects. Section III develops the proposed random access communication
scheme and we discuss how to adapt it to deal with energy
harvesting. In Section IV we introduce the algorithm used to
obtain the random access communication policy. The stability
of the system under the proposed policy is studied in Section
V. After this, we devote Section VI to simulations assessing
the performance of the proposed random access mechanism.
Finally, we provide some concluding remarks in Section VII.
II. S YSTEM M ODEL
Consider the system model shown in Figure 1. This scenario
consists of M different plants, which have their system state
measured by sensor nodes powered by energy harvesting.
The energy harvesting process imposes causality constraints
on the transmission capabilities of the sensors, as sensors
cannot transmit if they have not harvested sufficient energy.
The measurements collected by the sensor nodes have then
to be wirelessly transmitted to their respective controllers in
order to ensure plant stability. However, the wireless medium
over which the sensors transmit is shared. This implies that

We consider a group of M plants and use xi [t] ∈ Rni to
denote the state of the i-th plant at time t. Plant dynamics
are dictated by a linear time-invariant system in which plant
control is contingent on the successful reception of information
from the sensors. Define then the indicator variable γi [t] ∈
{0, 1} to signify with the value γi [t] = 1 that the transmission
of the i-th sensor at the t-th time slot has been successfully
received by the i-th controller. If information is successfully
received, we have γi [t] = 1, in which case the controller closes
the loop and the state evolves according to the closed loop
dynamics described by the matrix Ac,i ∈ Rni ×ni . If, on the
contrary, γi [t] = 0, the state evolves in open loop as described
by the matrix Ao,i ∈ Rni ×ni . We then have that the state xi [t]
evolves according to the switched state linear dynamics

xi [t + 1] =

(

Ac,i xi [t] + wi [t], if γi [t] = 1,
Ao,i xi [t] + wi [t], if γi [t] = 0,

(1)

where wi [t] correspond to independent and identically distributed (i.i.d.) Gaussian noise with covariance Ci . The design
of the controllers is not the focus of this paper. The matrices
are assumed given and are such that the closed loop matrix
Ac,i produces stable dynamics. The open loop matrix Ao,i
may produce stable or unstable dynamics but the problem is
of most interest when the open loop dynamics are unstable.
B. Communication Model
In the control model we have defined the variables γi [t] to
signify the successful reception of the sensor measurements.
This is a random variable whose distribution is dependent on
the chosen communication policy. We consider a time-slotted
communication model. At every time slot t, the i-th sensor
node decides to transmit with probability zi [t] ∈ [0, 1], where
we denote zi [t] as the scheduling variable. Then, if multiple
sensors transmit during the same time slot, we consider that
a collision occurs with probability qc ∈ [0, 1]. If a collision occurs, then none of the colliding packets are received.
Therefore, the probability of the i-th sensor node transmitting
at time slot t andQnot colliding with any other transmission
is given by zi [t] j6=i (1 − qc zj [t]). Apart from collisions,
packet loss can also occur due to incorrect decoding. The
probability of successful decoding is dependent on the channel
conditions at the i-th link during time slot t, which we denote
by hi [t]. Channel states are considered independent across the
M systems. Further, we consider a block fading model [26],
whereby the channel states hi [t] are i.i.d. over time slots and
constant during a time slot. The probability of successfully
decoding a packet given the channel state is denoted by
q(hi [t]), which is a continuous and strictly increasing function
q : R+ → [0, 1] (We show in Fig. 2 a typical decoding

3

Since this condition needs to hold for all xi [t], we can
equivalently rewrite this condition as the following linear
matrix inequality

1
q(hi [t])

0.8
0.6

ATc,i Pi Ac,i Pr(γi [t] = 1)

0.4

+ ATo,i Pi Ao,i (1 − Pr(γi [t] = 1)) ≤ ρi Pi ,

0.2
0

0

1

2

3

4

5

hi [t]
Fig. 2. Probability of decoding as a function of the channel state.

function). Also, for notational compactness, we also define
qi [t] , q(hi [t]). Then, the probability of successful reception
γi [t] is given by
Y
Pr (γi [t] = 1) = qi [t]zi [t]
(1 − qc zj [t]) .
(2)
j6=i

This expression simply corresponds to the successful decoding probability multiplied by the probability of transmitting
without colliding. Also, we assume that sensor nodes have
knowledge of their channel state before transmitting (In practice, this is usually achieved with pilot signals [26]).
C. Control Performance

The control loop of each plant is closed with a probability
given by equation (2). Since it is our objective to design
communication policies that satisfy a desired control performance, we aim to establish a relationship between the control
performance and the probability of successful reception. We
can do so by the following proposition.
Proposition 1 (Control performance abstraction [24]). Consider the switched system described by (1) with γi [t] given
by a sequence of i.i.d. Bernoulli random variables, and the
quadratic Lyapunov function Vi (xi ) = xTi Pi xi , with Pi ∈
Rni ×ni positive definite. Then the function Vi (xi ) decreases
at an average rate ρi < 1, denoted by
E[Vi (xi [t + 1])|xi [t]] ≤ ρi V (xi [t]) + tr(Pi Ci )

(3)

if and only if Pr(γi [t] = 1) ≥ pi , where pi is given by

pi = min θATc,i Pi Ac,i + (1 − θ)ATo,i Pi Ao,i ≤ ρi Pi

(4)

θ≥0

Proof. By particularizing the function Vi (xi ) =
the system dynamics (1), we can write the equation

xTi Pi xi

with

E[Vi (xi [t + 1])|xi [t]] = xi [t]ATc,i Pi Ac,i xi [t] Pr(γi [t] = 1)
+ xi [t]ATo,i Pi Ao,i xi [t] Pr(γi [t] = 0)
+ tr(Pi Ci ).

(5)

Then, by substituting this expression in the left hand side
of the average decrease inequality (3) we have the following
inequality
xi [t]ATc,i Pi Ac,i xi [t] Pr(γi [t] = 1)
+xi [t]ATo,i Pi Ao,i xi [t] Pr(γi [t] = 0) ≤ ρi xi [t]Pi xi [t]. (6)

(7)

where we have also used the fact that Pr(γi [t] = 0) = (1 −
Pr(γi [t] = 1)). Then, the Pr(γi [t] = 1) values satisfying this
inequality define a convex set of which there is a minimum
value pi such that the condition is equivalent to Pr(γi [t] =
1) ≥ pi .


This proposition allows us to establish a connection between
the control performance and the packet transmisisons. By
solving the semidefinite program (4), we obtain the successful
reception probabilities pi that allow us to satisfy the required
control performance. Then, we simply need to design communication policies that satisfy Pr(γi [t] = 1) ≥ pi for all
systems.
III. R ANDOM ACCESS C OMMUNICATION

We aim to design communication policies that satisfy the
successful packet reception probabilities given by Proposition 1. Under an assumption of ergodic processes, the successful packet reception probabilities are given by the long
term behavior of expression (2). Hence, in order to stabilize
the control system to the required control performance, the
scheduling variables zi [t] need to satisfy the following long
term constraint
t

Y
1X
qi [l]zi [l]
(1 − qc zj [l]) .
t→∞ t

pi ≤ lim

l=1

(8)

j6=i

Since we are working under the assumption of ergodicity, we
can write the previous limit as the expected value over channel
realizations. That is,


Y
pi ≤ E qi zi
(1 − qc zj ) ,
(9)
j6=i

and, since scheduling decision are independent over nodes, we
can further rewrite the previous expression as
Y
pi ≤ E [qi zi ]
(1 − E [qc zj ]) .
(10)
j6=i

Aside from the stabilization of all control loops, we also
want to minimize the number of times that a sensor node
accesses the medium. P
We do this by the introduction of
M
2
the objective function
i=1 E zi . Then, we formulate the
following optimization problem
minimize
zi ∈Z

M
X

E zi2

i=1

subject to pi ≤ E [qi zi ]

(11a)
Y
j6=i

(1 − E [qc zj ]) , i = 1, . . . , M
(11b)

where Z := {zi : R+ → [0, 1]} is the set of functions R+ →
[0, 1] taking values on [0, 1]. Notice that, while scheduling

4

decisions are statistically independent across sensors, solving
the previous optimization problem requires it being done in
a centralized manner (as constraint (11b) is coupled across
sensors). Nonetheless, we can separate the problem in a per
sensor manner by taking the logarithm of constraint (11b) as
follows
minimize
zi ∈Z,
sij ∈[0,1]

M
X

E zi2

(12a)

i=1

subject to

log(pi ) ≤ log(sii ) +
sii ≤ E qi zi ,

sij ≥ E qc zj ,

X
j6=i

log (1 − sij ) ,

i = 1, . . . , M

(12b)

i = 1, . . . , M

(12c)

i = 1, . . . , M, j 6= i

(12d)

where we have introduced the auxiliary variables sii and sij
and converted the logarithm of the product into a sum of logarithms. Solving the optimization problem (12) is equivalent
to solving (11). Under the assumption that this problem is
strictly feasible, that is, that
zi capable

Q there exist schedules
of satisfying pi < E qi zi j6=i (1 − qc zj ) , the goal is then
to design an algorithm such
that

 the instantaneous scheduling
decisions zi [t] satisfy E zi [t] = zi .
A. Random Access Communication with Energy Harvesting

We have proposed a random access optimization problem
that allows us to stabilize all control loops. However, the
formulation previously introduced does not account for either
the energy consumption nor the energy harvesting process.
We consider that the i-th sensor at time slot t acquires ei [t]
units of energy and stores it in a battery of finite capacity
bmax
. Further, we assume the energy harvesting process to be
i
stationary with mean E ei [t] . We consider that the sensor
nodes consume one unit of energy per channel access, hence
the scheduling variable also represents the power consumption
of a medium access. Then, in order to ensure that the sensor
nodes only use the energy available in their batteries, we have
the following energy causality constraint
(13)

zi [t] ≤ bi [t].

where bi [t] is the battery state of node i at time t. Further,
the battery of the nodes evolves according to the following
dynamics

bmax
i
bi [t + 1] = bi [t] − zi [t] + ei [t]
,
(14)
0

bmax
[·]0i

denotes the projection to the interval [0, bmax
].
where
i
However, these constraints are coupled across time slots and
cannot be directly introduced into the stochastic optimization
problem (12). In order to circumvent this, we consider the
long-term behavior of the energy causality constraints (13),
which, by recursively substituting the battery dynamics (14)
in (13) can be written as follows
t

t

l=1

l=1

1X
1X
lim
zi [l] ≤ lim
ei [l].
t→∞ t
t→∞ t

(15)

That is, in the long term, the battery state is dominated
by the harvested energy. Then, due to the ergodicity of the
scheduling variables zi [t] and the energy harvesting process
ei [t], the previous expression (15) can be simply written as
the expectation with respect to the channel states hi [t] and the
energy harvesting process ei [t], as follows
 
 
E zi ≤ E ei .
(16)
This constraint simply implies that, on average, the energy
spent for transmitting has to be lower than the harvested energy. Then, by introducing constraint (16) into the optimization
problem (12) we have the following problem
minimize
zi ∈Z,
sij ∈[0,1]

M
X

E zi2

(17a)

i=1

subject to

log(pi ) ≤ log(sii ) +

X
j6=i

log (1 − sij ) ,

i = 1, . . . , M
sii ≤ E qi zi ,

sij ≥ E qc zj ,
E zi ≤ E ei

(17b)

i = 1, . . . , M

(17c)

i = 1, . . . , M, j 6= i

(17d)

i = 1, . . . , M

(17e)

However, substituting the time slot to time slot constraints
(13) by the average ones (16) does not ensure that they are
satisfied at each time slot. This means that solutions to the
optimization problem (17) do not necessarily satisfy the energy
causality constraints zi [t] ≤ bi [t]. To overcome this problem,
and ensure causality, we introduce the following modified
problem formulation
minimize
zi ∈Z,
sij ∈[0,1],
yij ∈[0,ȳij ]

subject to

M
X

E zi2 +

i=1

M
M X
X

E ν̄ij yij

log(pi ) ≤ log(sii ) +

X
j6=i

log (1 − sij ) ,

i = 1, . . . , M
sii ≤ E qi zi + yii ,

sij ≥ E qc zj − yij ,
E zi ≤ E ei

(18a)

i=1 j=1

i = 1, . . . , M

(18b)
(18c)

i = 1, . . . , M, j 6= i (18d)

i = 1, . . . , M

(18e)

This optimization problem has been modified by the introduction of the auxiliary variables yij in constraints (18c)
and (18d), as well as in the objective function. The auxiliary variable yij is forced to take values in the interval
yij ∈ [0, ȳij ], where ȳij is a system-dependent constant.
PM PM
The term i=1 j=1 E ν̄ij yij in the objective function has
the constant ν̄ij , where the ν̄ij value is an upper bound on
the Lagrange multipliers of constraints (18c) and (18d). This
modified problem formulation allows us to ensure that even
though the energy constraint (18e) is on average form, the
energy causality constraints are satisfied in a time slot to time
slot basis, as we will show in the upcoming sections.
IV. R ANDOM ACCESS A LGORITHM
In this section, we aim to solve the optimization problem
(18). For notational compactness, let us define the vector z =

5

{zi , sij , yij } collecting all the primal variables and the vector
λ = {φi , νij , βi } collecting the dual variables. Further, we
collect the implicit primal variable constraints in the set X ,
{zi ∈ Z, sij ∈ [0, 1], yij ∈ [0, ȳij ]}. Then, the Lagrangian of
problem (18) can be written as
L(z, λ) =
+

M
X
i=1

M
X
i=1

+
+

E zi2 +

M
X

φi log (pi ) − log (sii ) −

νii (sii
i=1
M X
X

M
X
i=1

In a similar way, the auxiliary variables yij [t] can be computed
as the solution of the minimization

E ν̄ij yij

i=1 j=1



i=1 j6=i

+

M X
M
X

X
j6=i



log (1 − sij )

− E qi zi − yii )

νij (E qc zj − yij − sij )
(19)

βi (E zi − E ei ) .

yij [t] := arg min yij (ν̄ij − νij [t]) ,

(26)

yij ∈[0,ȳij ]

which is a thresholding condition. The auxiliary variable yij [t]
takes the value yij [t] := 0 if νij [t] ≤ ν̄ij and yij [t] := ȳij if
νij [t] > ν̄ij . Next, since the dual function is concave, we
can perform a subgradient ascent on the dual domain. The
corresponding dual variable updates are given by


φi [t + 1] := φi [t] +  log (pi ) − log (sii [t])
+
X
log (1 − sij [t])
−
(27)
j6=i

The Lagrange dual function of this problem is given by
(20)

g(λ) = min L(z, λ).
z∈X

Note that, while the primal problem is infinite dimensional, the
dual problem has a finite number of variables (the dual variables). Furthermore, for this problem, the duality gap can be
shown to be zero [27]. Hence, we resort to a dual subgradient
method to solve the optimization problem. However, the sensor
nodes have no knowledge of the probability distribution over
which the expectation is taken. In order to overcome this, we
substitute the random variables by their instantaneous values,
which are known by the sensors. Finally, by reordering the
Lagrangian (19), the scheduling variables zi [t] are given by
the following minimization


X
zi [t] := arg min zi zi − νii [t]qi [t] + qc
νji [t] + βi [t] ,
zi ∈[0,1]

j6=i

(21)

which is separated across sensors and leads to the following
closed form solution

1
X
1
νii [t]qi [t] − qc
νji [t] − βi [t] .
(22)
zi [t] :=
2
0
j6=i

The resulting optimal scheduling policy is to transmit at time
slot t with the probability given by (22). This is a policy
that dynamically adapts to the time-varying conditions of the
system. Namely, the dual variables νii [t] and νji [t] depend on
the stability of all the plants, the qi [t] variable is dependent
on the channel state hi [t], and the dual variables βi [t] depend
on the energy harvesting process ei [t]. The rest of the primal
variables sii [t] and sij [t] can be found by the minimization
sii [t] := arg min −φi [t] log(sii ) + νii [t]sii ,

(23)

sij [t] := arg min −φi [t] log(1 − sij ) − νij [t]sij ,

(24)

sii ∈[0,1]

sij ∈[0,1]

which, again, are separated across sensors and have the following closed forms solutions

1

1
φi [t]
φi [t]
sii [t] :=
,
sij [t] := 1 −
.
(25)
νii [t] 0
νij [t] 0



+
νii [t + 1] := νii [t] +  sii [t] − zi [t]qi [t] − yii [t]
(28)


+
νij [t + 1] := νij [t] +  qc zi [t] − sij [t] − yij [t]


+
βi [t + 1] := βi [t] +  zi [t] − ei [t]

(29)

(30)

where, in order to have an algorithm than can be run in an
online manner, we have considered a fixed step sized . For
notational compactness, we also write the dual update
+ in a
concatenated vector form as λ[t + 1] := λ[t] + s[t] , where
si [t] corresponds to the stochastic subgradient. The steps in
the resulting random access mechanism are summarized in
Algorithm 1.
Also, it is important to note that we can establish a parallel
relationship between the dual variables βi [t] associated to the
energy constraint E zi ≤ E ei and the actual battery state bi [t].
Thisrelationship is given by the expression βi [t] =  bmax
−
i
bi [t] . Hence, a mirrored symmetry (scaled by the step size )
exists between these variables. This relationship will be crucial
in ensuring the energy causality of the algorithm, as we show
next.
A. Energy Causality
Now, we turn our attention to the study of the conditions
required to satisfy the causality constraints, i.e., zi [t] ≤ bi [t]
for all time slots. We have introduced the modified problem
formulation (18) to help achieve this. First, we show that
this problem formulation allows us to upper bound the dual
variables νij [t] over all time slots.
Proposition 2. Let the upper bound ȳij of the auxiliary
variables yij satisfy the inequality ȳij ≥ 1 (ν̄ij + 2). Then,
the dual variables νij [t] are upper bounded by νij [t] ≤ ν̄ij + 
for all time slots t.

6

Algorithm 1 Random access scheduling algorithm.
1: Initialize: Initialize the dual variables to φi [0] := 0,
− bi [0] .
νij [0] := 0, and βi [0] :=  bmax
i
2: Step 1: Medium access decision

1
P
3: zi [t] := 21 νii [t]qi [t] − qc
ν
[t]
−
β
[t]
i
j6=i ji
Step 2: Other
variables h
iprimal
h
1
φi [t]
and sij [t] := 1 −
5: sii [t] := ν [t]
ii
0
6: Step 3: Auxiliary variable
7: yij [t] := arg max yij (ν̄ij − νij [t])

0

4:

φi [t]
νij [t]

i1
0

yij ∈[0,ȳij ]

8:
9:

10:
11:
12:
13:

Step 4: The sensor updates the dual variables
φi [t + 1] := φi [t] +  log (pi ) − log (sii [t])

 +
P
−
log (1 − sij [t])
j6=i


 +
νii [t + 1] := νii [t] +  sii [t] − zi [t]qi [t] − yii [t]


 +
νij [t + 1] := νij [t] +  qc zi [t] − sij [t] − yij [t]


 +
βi [t + 1] := βi [t] +  zi [t] − ei [t]
Step 5: Set t := t + 1 and go to Step 1.

Proof. The dual variable νij [t] is updated according to the
following equations


+
νii [t + 1] := νii [t] +  sii [t] − zi [t]qi [t] − yii [t]
(31)


+
νij [t + 1] := νij [t] +  qc zi [t] − sij [t] − yij [t]
, (32)

where the subgradient terms are upper bounded by 1, namely
sii [t]−zi [t]qi [t]−yii [t] ≤ 1 and qc zi [t]−sij [t]−yij [t] ≤ 1 for
all time slots t. Hence, we have that the maximum increase of
these dual variables in a given time slot is |νij [t+1]−νij [t]| ≤
 for all i, j and t. Overall, the maximum value that the dual
variables νij [t] can take is controlled by the yij [t] term. As
long as yij [t] = 0 the dual variables can increase in value,
until νij [t] = ν̄ij +  and the auxiliary variable condition in
(26) is triggered, leading to the yij [t] term taking the value
yij [t] = ȳij . Then, the next update of the dual variable is
given by

+
νij [t + 1] ≤ ν̄ij +  − yij [t]


+
1
≤ ν̄ij +  +  − 
(ν̄ij + 2)

= 0.
(33)
Since after this event, the dual variables take the zero value, the
dual variables νij [t] are necessarily upper bounded by νij [t] ≤
ν̄ij +  for all time slots t.

This proposition states that by ensuring the correct value of
the parameter ȳij (which we can select freely), an upper bound
on νij [t] can be established. Then, by further appropriately
selecting the battery size bmax
of the nodes, we can ensure
i
that energy use is causal to the energy harvested.

Proposition 3 (Energy Causality). Let the battery capacity of the i-th sensor satisfy bmax
≥ 1 ν̄ii + 1 and let
i
1
ȳij ≥  (ν̄ij + 2). Then, Algorithm 1 satisfies the energy
consumption causality constraints zi [t] ≤ bi [t] for all time
slots.
Proof. In order to satisfy the energy causality constraints
zi [t] ≤ bi [t], it suffices to verify that no transmission occurs
when there is no energy left in the battery. This implies that the
scheduling variable zi [t] has to take the zero value when the
battery bi [t] is P
empty. By equation (22), it suffices to satisfy
νii [t]qi [t] − qc j6=i νji [t] − βi [t] ≤ 0, when bi [t] = 0. Note
that the battery state bi [t] and the battery multipliers
βi [t] are

−
b
[t]
.
Hence,
the
related by the expression βi [t] =  bmax
i
i
battery being empty, bi [t] = 0, implies the battery multipliers
taking the value βi [t] = bmax
i P. Therefore, the condition to
be satisfied is νii [t]qi [t] − qc j6=i νji [t] − bmax
≤ 0. Since
i
qi [t] ≤ 0 and qc ≥ 0, we can further rewrite this inequality as
≤ 0. Then, by Proposition 2, we have the upper
νii [t] − bmax
i
bound on the dual variables νij [t] ≤ ν̄ij + . This allows us
to further rewrite the inequality as ν̄ij +  − bmax
≤ 0. Then,
i
1
ν̄
+
1,
ensures
this
inequality,
and
the battery size bmax
≥
ii
i

hence, that the energy constraints zi [t] ≤ bi [t] are satisfied for
all time slots.

According to this proposition, by choosing a sufficiently
large battery size bmax
we can make the energy consumption
i
causal to the energy harvesting process. This is due to the
modified problem formulation proposed in (18). In the original
problem (17), a dual ascent algorithm can lead to the dual
variables becoming arbitrarily large. This is not the case
when introducing the auxiliary variables yij , as shown by
Proposition 2. Then, by Proposition 3, the bound on the dual
variables allows us to establish conditions on the battery size
that ensure energy causality.
B. Asynchronous Operation
In order for Algorithm 1 to properly function, the i-th
sensor node requires the dual variables νij [t] for j 6= i of
the other nodes. This is needed in order to compute the
optimal scheduling variable zi [t] (cf. equation (22)). In order
to ensure the robustness of our algorithm, we take into account
the notion of asynchronicity in the data shared across the
sensor nodes. This is to say that we consider the possibility of
different nodes having different (out of date) values of the dual
variables shared by the other nodes. Since nodes are powered
by energy harvesting, this might happen when a node is unable
to transmit or receive data due to lack of energy. Also, the
consideration of asynchronicity includes the practical case in
which the sensor nodes simply attach the value of their dual
variable to the packet containing the measurement. Therefore,
only sharing their dual variable when they need to transmit a
measurement to their controller. To consider this, we introduce
the asynchronicity model of [28] into our analysis.
Let us define the set T i ⊆ Z+ of all time slots in which
the i-th node is capable of receiving and sending information.
Then, we define a function π i [t], that for a given node and

7

time slot, returns the most recent time slot at which the node
was available. Namely,

π i [t] := max t̂ | t̂ < t, t̂ ∈ T i .
(34)

In a similar
manner, we then define the function πji [t] :=

j
i
π π [t] to denote the most recent time slot the i-th node
has received information sent by the j-th node. Then, at time
slot t, the i-th node has knowledge of a possibly outdated
vector of dual variables νij , given by

i
ν̃i [t] = νi1 [π1i [t]], . . . , νiM [πM
[t]] .
(35)

Further, we will denote by λ̃ the vector formed by the collection of the outdated duals together with the rest of the dual
variables. Following, we can write the asynchronism into the
dual variable update by defining the asynchronous stochastic
subgradient s̃νij [t], corresponding to the νij variable. We do
so as follows
(
sνij [t], if t ∈ T j ,
(36)
s̃νij [t] =
0,
otherwise.

Simply meaning that, if t ∈ T j , the ascent direction given
by the subgradient sνij [t], is available. Otherwise, the dual
variable is not updated. Then, we can concatenate all the subgradients of all dual variables into an asynchronous stochastic subgradient vector s̃[t]. Afterwards, the dual variable
update is simply given by the usual expression but with
the asynchronous stochastic subgradient, i.e., λ[t + 1] :=
+
[λ[t] + s̃[t]] .
V. S TABILITY A NALYSIS
In this section, we analyze the stability of the systems when
operating under the proposed random access communication
scheme. In order to do this, we leverage on the fact that
the proposed scheme is a stochastic subgradient algorithm.
Hence, we rely on duality theory arguments to show that the
iterates generated by Algorithm 1 satisfy the constraints of
the optimization problem (18) almost surely. Then, we further
show that if the constants ν̄ij are chosen to upper bound the
?
optimal Lagrange multipliers νij
, then the iterates generated
by Algorithm 1 also satisfy the constraints of the optimization
problem (17) (i.e., without the auxiliary variables). In turn, this
guarantees by Proposition 1 the stability of all control loops.
First, in order to ensure the convergence of Algorithm 1, we
need to assume an upper bound on the asynchronicity between
the sensor nodes.
Assumption 4. There exists an upper bound 0 < B < ∞ to
the asynchronicity between nodes, such that for all time t and
nodes i, j we have

Proposition 5. Given
dual

 variables λ[t], the conditional
expected value E s[t]|λ[t] of the stochastic subgradient s[t]
is a subgradient of the dual function. Namely, for any λ,



E sT [t]|λ[t] λ[t] − λ ≤ g(λ[t]) − g(λ).
(38)
Proof. We intend to show that the expected value of the
stochastic subgradient s[t] given λ[t] is a subgradient of the
dual function g(λ). To do this, we take the Lagrangian (19)
of optimization problem (18), given by
X
X
X
M
2
M
M
L(z, λ) =
i=1 E zi +
i=1
j=1 E ν̄ij yij
X
X

M
+
j6=i log (1 − sij )
i=1 φi log (pi ) − log (sii ) −
X

M
+
i=1 νii sii − E qi zi − yii
X
X

M
+
j6=i νij E qc zj − yij − sij
i=1
X

M
+
(39)
i=1 βi E zi − E ei .

Then, take the dual function at time t, denoted by g(λ[t])
and remember that the dual function is given by g(λ) =
minz∈X L(z, λ). The primal variables that minimize this dual
function are obtained by the primal minimization of Algorithm
1, namely, zi [t], sij [t] and yij [t], given by equations (22), (25),
and (26), respectively. Then, we write the dual function at time
t as
X
X
X
M
2
M
M
g(λ[t]) =
i=1 E zi [t] +
i=1
j=1 E ν̄ij yij [t]+
X
X

M
j6=i log (1 − sij [t])
i=1 φi [t] log (pi ) − log (sii [t]) −
X


M
+
i=1 νii [t] E sii [t] − qi [t]zi [t] − yii [t]
X
X


M
+
j6=i νij [t] E qc zj [t] − yij [t] − sij [t]
i=1
X


M
+
(40)
i=1 βi [t] E zi [t] − ei [t] ,
where, due to its linearity, we have moved the expectation E[·]
out of the subgradients. Now, by compacting the Lagrange
multipliers into a vector λ[t] and the subgradients to s[t], we
can rewrite the dual function at time t as
X
X
X
M
2
M
M
g(λ[t]) =
i=1 E zi [t] +
i=1
j=1 E ν̄ij yij [t]
 T

+ E s [t]|λ[t] λ[t].
(41)
Further, for any arbitrary λ, the dual function g(λ) can be
bounded as
X
X
X
M
2
M
M
g(λ) ≤
i=1 E zi [t] +
i=1
j=1 E ν̄ij yij [t]
 T

+ E s [t]|λ[t] λ.
(42)

(37)

Then, by substracting expression (42) from (41) we obtain



E sT [t]|λ[t] λ[t] − λ ≤ g(λ[t]) − g(λ),
(43)

This assumption simply implies that nodes are at most B
time slots out of synchronism and it is required to ensure the
convergence of the variables. Now, we proceed to show the
convergence of Algorithm 1. We start by recalling a common
property of the subgradient method.

The previous proposition states that, on average, the stochastic subgradient s[t] is an ascent direction of the dual function
g(λ[t]). Then, the next step is to quantify the average reduction
in distance to the optimal dual variables λ? that occurs in a
dual variable update step.

max {0, t − B + 1} ≤ πji [t] ≤ t.

which is the desired inequality.



8



Lemma 6. Let E ks[t]k2 |λ[t] ≤ S 2 be a bound on the second
moment of the norm of the stochastic subgradients s[t]. The
dual updates of Algorithm 1, satisfy the following inequality



2
2
E kλ? −λ[t + 1] |λ[t] ≤ 1 − m + 22 LB λ? − λ[t]

+ 2 S 2 + 22 LBS 2 −  g(λ? ) − g(λ[t]) ,
(44)

where the constant L > 0 corresponds to the L-Lipschitz
continuity of the gradients of the dual function g(λ) and
m > 0 to the strong concavity constant of the dual function
g(λ).
Proof. Let us consider the squared distance between the dual
iterates λ at time t+1 and their optimal value, i.e., λ? −λ[t+
2
+
1] . By means of the dual update λ[t + 1] = [λ[t] + s̃[t]] ,
we can rewrite this expression as

+ 2
2
λ? − λ[t + 1] = λ? − λ[t] + s̃[t]
2

≤ λ? − λ[t] − s̃[t]

(45)

where we have further upper bounded the expression by the
nonexpansive property of the nonnegative projection. Then,
we expand the square norm, yielding the expressions
λ? − λ[t + 1]

2

≤ λ? − λ[t]

2

+ 2 s̃[t]

2


− 2s̃T [t] λ? − λ[t] . (46)

We can further rewrite this inequality
by adding and subtract
ing the term 2sT [t] λ? − λ[t] to expression (46), leading to
the following
λ? − λ[t + 1]

2

2

2

≤ λ? − λ[t] + 2 s̃[t]
T ?


+ 2 s[t] − s̃[t]
λ − λ[t] − 2sT [t] λ? − λ[t] . (47)

By applying the Cauchy-Schwarz inequality to the third term
on the right hand side, we further rewrite the expression as
λ? − λ[t + 1]

2

+ 2 s[t] − s̃[t]

≤ λ? − λ[t]

2

+ 2 s̃[t]

2



λ? − λ[t] − 2sT [t] λ? − λ[t] . (48)

2

≤ λ? − λ[t]

+ 2L λ[t] − λ̃[t]

+ 2 s̃[t]

2


λ? − λ[t] − 2sT [t] λ? − λ[t] (49)

Now, we bound the difference λ[t] − λ̃[t] between the dual
variables λ[t] and their asynchronous counterpart λ̃[t]. Recall
that Assumption 4 states that there exists an asynchronicity
limit of B time slots between the global and asynchronous
variables. This means that the asynchronous dual variables
λ̃[t] are, at most, B subgradients steps out of synchronysm.
We can then bound the difference by λ[t] − λ̃[t] ≤
Pt−1
Pt−1

l=t−B−1 s̃[l] , where we have also
l=t−B−1 s̃[l] ≤ 
applied the triangle inequality. Then, we write
?

λ − λ[t + 1]

2

?

≤ λ − λ[t]
+ 22 L

2

t−1
X

+

l=t−B−1
?

2

s̃[t]

s̃[l]


− 2sT [t] λ − λ[t] .

λ? − λ[t + 1]

2

≤ λ? − λ[t]

2

t−1
X

+ 22 L

2

+ 2 s̃[t]

2
s̃[l] + λ? − λ[t]

l=t−B−1
?


− 2s [t] λ − λ[t] .
T

2

λ? − λ[t]

!
(50)

2



(51)

Rearranging the terms
λ? − λ[t + 1]

2

≤ λ? − λ[t]
+ 2 s̃[t]

2

2

+ 22 LB λ? − λ[t]

+ 22 L

t−1
X

s̃[l]

2

2

l=t−B−1

T

?

− 2s [t] λ − λ[t]



(52)

Then, separate the last term on the right hand side, and take
−sT [t] λ? − λ[t] and note that we can rewrite it as − s? −
T ?

s[t]
λ − λ[t] . Then, by strong concavity we can bound
2
this term by −m λ? − λ[t] . Now, we take the expectation
conditioned on λ[t] on both sides of the previous inequality



2
2
E kλ? − λ[t + 1] |λ[t] ≤ 1 − m + 22 LB λ? − λ[t]


+ 2 E ks̃[t]k2 |λ[t] + 22 L

t−1
X

l=t−B−1




−  E sT [t]|λ[t] λ? − λ[t]



E ks̃[l]k2 |λ[t]

(53)

And
 then, by applying the subgradient bound given by
E ks[t]k2 |λ[t] ≤ S 2 , and particularizing Proposition 5 with
λ = λ? , we have



2
2
E kλ? −λ[t + 1] |λ[t] ≤ 1 − m + 22 LB λ? − λ[t]

+ 2 S 2 + 22 LBS 2 −  g(λ? ) − g(λ[t]) ,
(54)

which gives us the desired inequality.

2

Then, we can further bound the expression by relying on the
L-Lipschitz continuity of the subgradients of the dual function
λ? − λ[t + 1]

The third term on the right hand side can be further expanded
by making use of the inequality kukkvk ≤ kuk2 + kvk2 ,
leading to the following expression



The previous lemma holds on average, while we are interested in establishing convergence almost surely. We leverage
on this lemma and resort to a supermartingale argument to
show that Algorithm 1 converges to a neighborhood of the
optimal solution of the dual function.


Lemma 7. Let E ks[t]k2 |λ[t] ≤ S 2 be a bound on the
second moment of the norm of the stochastic subgradients
s[t]. Further, consider the dual updates of Algorithm 1, with
step size  ≤ m/(2LB). Then, assume that the dual variable λ[T ] is given for an arbitrary time T and define as
λbest [t] := arg minλ[l] g(λ[l]) the dual variable leading to the
best value of the of the dual function for the interval l ∈ [T, t].
Then, we have
lim g(λbest [t]|λ[T ]) ≥ g(λ? ) − 2 S 2 − 22 LBS 2 a.s. (55)

t→∞

Proof. Let T = 0 for simplicity of exposition. Then, define
the sequence α[t] corresponding to a stopped process tracking
2
the dual distance λ? −λ[t] until the optimality gap g(λ? )−
g(λ[t]) falls bellow S 2 + 2LBS 2 . Namely,
2

α[t] : = λ? − λ[t]

I g(λ? ) − g(λbest [t]) > S 2 + 2LBS 2 ,

(56)

9

where I{·} is the indicator function. In a similar manner define
the sequence β[t] as follows



β[t] : =  g(λ? ) − g(λ[t]) − 2 S 2 − 22 LBS 2

I g(λ? ) − g(λbest [t]) > S 2 + 2LBS 2 .
(57)

Now, let F[t] be the filtration measuring α[t], β[t] and λ[t].
Since α[t] and β[t] are completely determined by λ[t], and
λ[t] is a Markov process, conditioning on F[t] is equivalent
to conditioning on λ[t]. Hence, we can write the expectation
E [α[t]|F[t]] = E [α[t]|λ[t]]. Now, consider this expectation at
time t + 1, given by



2
E α[t + 1]|λ[t] = E λ? − λ[t + 1]
 ?

I g(λ ) − g(λbest [t + 1]) > S 2 + 2LBS 2 |λ[t] . (58)

By noting that the indicator term is lower or equal than 1, we
can upper bound this expression as




2
(59)
E α[t + 1]|λ[t] ≤ E λ? − λ[t + 1] |λ[t] .
Then, by application of Lemma 6 we have



E α[t + 1]|λ[t] ≤ 1 − m + 22 LB λ? − λ[t]

+ 2 S 2 + 22 LBS 2 −  g(λ? ) − g(λ[t]) .

2

(60)

By making use of the definitions of α[t] and β[t], given by
equations (56) and (57), we can rewrite the previous expression
as



E α[t + 1]|λ[t] ≤ 1 − m + 22 LB α[t] − β[t].
(61)

Now, since the step size chosen is  ≤ m/(2LB), this means
that the factor multiplying
the process α[t] is upper bounded

by 1 − m + 22 LB ≤ 1. Therefore, we can further write
the expectation as


E α[t + 1]|λ[t] ≤ α[t] − β[t].
(62)
Since by definition the processes α[t] and β[t] are nonnegative,
the supermartingale convergence theorem [29, Theorem 5.2.9]
statesP
that the sequence α[t] converges almost surely, and the
∞
sum t=1 β[t] < ∞ is almost surely finite. This carries the
implication that lim inf t→∞ β[t] = 0. Given the definition
of the sequence β[t], this implies that limt→∞ g(λbest [t]) ≥
g(λ? ) − 2 S 2 − 22 LBS 2 almost surely.


Now, it suffices to show that the iterates generated by the
algorithm are almost surely feasible to the original problem
(17) if the constants ν̄ij are chosen to be an upper bound of
?
the optimal νij
multipliers.

Proposition 8 (Auxiliary Feasibility). Assume there exist
feasible primal variables zi , sij andPyij , such that for some
ξ > 0, we have log(pi ) − log(sii ) − j6=i log (1 − sij ) < −ξ,
sii − E qi zi − yii < −ξ, E qc zj − yij − sij < −ξ and
E zi −E ei < −ξ. Then, the sequences generated by Algorithm
1, satisfy the constraints (18b)−(18e) almost surely.
Proof. First, we start by upper bounding the value of the dual
variables. We collect the feasible primal variables in a vector
ẑ = {zi , sij , yij }. Then, given feasible primal variables ẑ we
bound the value of the dual function g(λ). Recall that the dual

function is defined as g(λ) = minz∈X L(z, λ), then for any
feasible ẑ, we necessarily have g(λ) ≤ L(ẑ, λ). Hence, we
can write
X
X
X
M
2
M
M
g(λ) ≤
i=1 E zi +
i=1
j=1 E ν̄ij yij
X
X

M
+
j6=i log (1 − sij )
i=1 φi log (pi ) − log (sii ) −
X

M
+
i=1 νii sii − E qi zi − yii
X
X

M
+
j6=i νij E qc zj − yij − sij
i=1
X

M
+
(63)
i=1 βi E zi − E ei .
Since
P we have a constant ξ > 0 such that log(pi ) − log(sii ) −
j6=i log (1 − sij ) < −ξ, sii − E qi zi − yii < −ξ, E qc zj −
yij − sij < −ξ and E zi − E ei < −ξ. We can simplify the
bound on the dual function to the following inequality
g(λ) ≤

M
X
i=1

E zi2 +

M X
M
X
i=1 j=1

E ν̄ij yij − ξλT 1.

(64)

Then, we simply reorder the previous expression to establish
an upper bound on the dual variables,

M
M
M X
X
1 X
2
λ≤
E ν̄ij yij − g(λ) ,
(65)
E zi +
ξ i=1
i=1 j=1
where the inequality is taken component-wise for all elements
of the vector λ with respect to the scalar on the right hand
side of the inequality. Then, by Lemma 7 we can certify the
existence of a time t ≥ T1 such that g(λ[t]) ≥ g(λ? ) − 2 S 2 −
22 LBS 2 . Hence, we can write
M
M
M X
X
1 X
λ[t] ≤
E ν̄ij yij
E zi2 +
ξ i=1
i=1 j=1

− g(λ? ) + 2 S 2 + 22 LBS 2
(66)
for some t ≥ T1 . Now, consider the feasibility conditions
of the optimization problem with the auxiliary variables
(18b)−(18e), which are given by the long term behavior of
the following inequalities
t
X

1X
log(pi ) − log(sii [l]) −
log (1 − sij [l]) ≤ 0
t→∞ t

lim

l=1

1
t→∞ t
lim

1
t→∞ t
lim

lim

t→∞

1
t

t
X
l=1
t
X
l=1
t
X
l=1

j6=i


sii [l] − qi [l]zi [l] − yii [l] ≤ 0

qc zj [l] − yij [l] − sij [l] ≤ 0

zi [l] − ei [l] ≤ 0.

(67)

(68)
(69)
(70)

These inequalities simply correspond to the stochastic subgradients of the optimization problem (18). Therefore, we
can write P
the feasibility conditions in a condensed form as
t
limt→∞ 1t l=1 s[l] ≤ 0. Now, we consider the dual updates

10


+
of the problem, given by λ[t + 1] = λ[t] + s̃[t] . Since the
projection is nonnegative, we can upper bound λ[t + 1] by
λ[t + 1] ≥ λ[t] + s̃[t] ≥ λ[1] + 

t
X
l=1

s̃[l] ≥ 

t
X

s̃[l], (71)

l=1

where we have removed the projection and further upper
bounded the expression by recursively substituting the dual updates. Now, we proceed to prove feasibility of the constraints
of the auxiliary problem (18). In order to to this, we follow
by contradiction. Assume that the conditions (67)−(70) are
infeasible. This means that there exists some time T2 where
there is a P
constant δ > 0 such that for t ≥ T2 we have
t
limt→∞ 1t l=1 s[l] ≥ δ. Substituting this expressions in the
dual update bound (71) we have that λ[t + 1] ≥ δt. Then, we
can freely choose a time t ≥ T2 such that
M
M X
M
X
1 X
λ[t] >
E zi2 +
E ν̄ij yij
ξ i=1
i=1 j=1

?
2 2
2
2
− g(λ ) +  S + 2 LBS
(72)
However, the upper bound established in (66) contradicts this
expression. Therefore, the inequalities (67)−(70) are satisfied
almost surely, which implies that the constraints (18b)−(18e)
of the auxiliary optimization problem (18) are almost surely
satisfied.


This proposition allows us to certify that the constraints of
the problem with the auxiliary variables are satisfied. However,
this does not ensure stability. Nonetheless, we show that if
the parameters ν̄ij are chosen as to upper bound the optimal
dual variables νij , then the optimal auxiliary variables are
zero, and the proposed algorithm also satisfies the constraints
(17b)−(17e) of the original problem without the auxiliary
variables.
Theorem 9 (Stability). Assume there exist feasible primal
variables zi and sijP
, such that for some ξ > 0, we have
log(pi ) − log(sii ) − j6=i log (1 − sij ) < −ξ, sii − E qi zi <
−ξ, E qc zj − sij < −ξ and E zi − E ei < −ξ. Further, let ν̄ij
be an upper bound to the optimal νij multipliers. Then, the
scheduling decisions zi [t] generated by Algorithm 1 satisfy the
successful packet transmission requirement
t

Y
1X
qi [l]zi [l]
(1 − qc zj [l]) > pi ,
lim
t→∞ t
l=1

(73)

j6=i

which ensures the stability of all control loops.
Proof. To verify this, subtract the Lagrangian of the optimization problem with the auxiliary variables (18) and the original
problem (17). This difference is given by,
L(z, λ) − L̂(z, λ) =
−

M X
M
X
i=1 j=1

M X
M
X
i=1 j=1


ν̄ij − νij − θij + µij yij
µij ȳij

(74)

where θij ≥ 0 and µij ≥ 0 are the Lagrange multipliers of
the implicit constraints yij ≥ 0 and yij ≤ ȳij , respectively. To
certify the equivalence between the optimization problems (17)
and (18) we need to certify that (74) is zero for the optimal
values. This implies that there must exist Lagrange multipliers
such that ν̄ij −νij −θij +µij = 0 and µij = 0 for all i, j. Since
?
≤ ν̄ij , we can find multipliers satisfying the constraints by
νij
?
?
. Then, L(z, λ) − L̂(z, λ) =
= ν̄ij − νij
letting µ?ij = 0 and θij
0, which implies
the
solution
of
both
problems is equal. Since
Pt
?
?
limt→∞ 1t l=1 yij [l] = yij
and yij
= 0, the primal variables
zi and sij almost surely satisfy the original constraints of the
optimization problem without the auxiliary variable, given by
(17b)−(17e). Since constraints

 are equivalent to
Q (17b)−(17d)
the constraint pi < E qi zi j6=i (1 − qc zj ) , by Proposition 1
we have that Algorithm 1 generates schedules zi [t] that ensure
the stability of all control loops.

Theorem 9 states that if there exist schedules zi [t] capable of
stabilizing the plants, then Algorithm 1 generates them. Note
that the optimal auxiliary variables yij take the zero value
and are therefore not necessary in the long run to satisfy the
constraints of the original optimization problem. As previously
discussed, they have been introduced in order to provide a way
to bound the dual variables and allow us to enforce causality
in the energy consumption.
VI. N UMERICAL R ESULTS
In this section, we study the performance of the proposed
random access scheme with energy harvesting sensors. We
consider a scalar control system, with M = 2 plants over
T = 10,000 time slots. The plant dynamics are given by
Ao,1 = 1.1 and Ac,1 = 0.15 for the first plant, and Ao,2 =
1.05 and Ac,2 = 0.1 for the second one. Hence, the first
system is slightly more unstable than the second one. Further,
we consider both systems to be perturbed by i.i.d. zeromean Gaussian noise. Also, we assume the same performance
requirement for both plants, given by the Lyapunov function
Vi (xi [t]) = x2i [t] and an expected decrease rate of ρi = 0.8.
With regards to the communication aspects, we consider a
system where the channel
hi [t] are i.i.d. exponential
 states

variables with mean E hi [t] = 2, and the decoding probability q(hi [t]) is given by the function shown in Figure 2.
Since the communication medium is shared, we consider that
packet collisions occur with probability qc = 0.25. Moreover,
we consider the sensing devices
  to be powered by an energy
harvesting process of rate E ei = 0.5, and that they store the
collected energy in batteries of size bmax
= 20. Finally, the
i
parameters of the algorithm are chosen as ȳij = 25, ν̄ij = 19,
and step size  = 1.
A. System Dynamics
We start by studying the evolution of the system dynamics.
In Figure 3, we plot the evolution of the plant state at each
time slot. As expected, since our proposed policy stabilizes
both plants, the system state oscillates around the zero value.
Furthermore, this figure illustrates that System 1 is slightly
more unstable than System 2. This is evidenced by the

11

20

10

System 1
System 2

8
kxi [l]k2

xi [t]

10

4

1
t

Pt

l=0

0

6

−10
−20

2
0

0.2

0.4

0.6

0.8

0

1
·104

Time (t)
Fig. 3. Evolution of the plant state at each time slot.

System 1
System 2
0

0.2

0.4

0.6

Time (t)

1
·104

Fig. 5. Evolution of the system control performance.

20

0.3

System 1
System 2

Energy Balance

15
bi [t]

0.8

10

5

0.2

0.1

System 1
System 2
0

0

0.2

0.4

0.6

Time (t)

0.8

1
·10

4

Fig. 4. Energy stored in the batteries at each time slot.

somewhat more pronounced peaks of instability, and higher
variance of the plant state x1 [t]. In a similar manner, this
behavior is also shown in Fig. 4. In this figure, we have plotted
the evolution of the battery state of the sensing devices of
both systems. Since the first plant is slightly more unstable
than the second plant, the energy consumption of the sensor
in the first system is slightly more pronounced. Also, by taking
a look at this figure, one can expect the first system to have
larger battery requirements than the second system. Intuitively,
since System 1 is more unstable, its sensor node requires a
larger battery to counteract its instability. The extent of this
requirement will become more apparent once we take a look
at the values of the dual variables.
B. Stability
In order to study the stability of the plants under our proposed policy, we look at the long term evolution of the system
states. First, we look at the evolution of the system control
performance. By our design, we require the Lyapunov function
Vi (xi [t]) = x2i [t] to decrease at an average rate of ρi = 0.8.
Iterating expression (3) in Proposition 1, we have that the limit

0

0

0.2

0.4

0.6

Time (t)

0.8

1
·10

4

Fig. 6. AveragePenergy balance in the system over time. Given by the
expression (1/t) tl=0 (ei [l] − zi [l]).

of the control performance is upper bounded in the long run by
the term tr(Pi Wi )/(1 − ρi ). By particularizing this expression
to our parameters, we expect the control performance in the
limit to be below tr(Pi Wi )/(1 − ρi ) = 1/(1 − 0.8) = 5. We
plot in Figure 5 the system control evolution. As expected,
both systems are asymptotically stable andP
the control perforT
mance converges to approximately (1/T ) t=0 kxi [t]k2 ≈ 4
for both plants.
Another interesting measure to consider is the energy
balance of the systems. We denote by energy balance the
difference between the harvested energy and the consumed
energy. Thus, the average energy
Pt balance of the i-th system
at time t is given by (1/t) l=0 (ei [l] − zi [l]). We plot this
measure in Figure 6. As expected, since both sensors are
powered
by energy harvesting processes of the same mean
 
E ei = 0.5 and System 1 is more unstable than System 2,
the energy balance of the first system is lower. Also, note
that the lower bound on the energy balance is zero, since
the total energy spent has to necessarily be lower or equal
to the total energy harvested. This allows us to interpret the

12

10

Transmission Probability

8
6

Pt

l=0

νij [l]

0.5

ν11
ν12
ν21
ν22

1
t

4
2
0

0

0.2

0.4

0.6

0.8

Time (t)

1

0.4
0.3
0.2
System 1
System 2

0.1
0

0

0.2

·104

0.4

0.6

Time (t)

Fig. 7. Average evolution of the dual variables νij over time.

Fig. 8. Average transmission probabilities.

energy balance as a measure of how much more control
performance can be obtained with the same energy harvesting
process. For example, System 1 has an energy balance of
approximately 0.05 units. Since we have assumed an unitary
power consumption, this means that System 1 has energy to
support an increase by 0.05 of its transmit probability. In the
same manner, System 2 can support an increase of around
0.14 of its transmit probability.
Now, we take a look at the dual variables. As we have
thoroughly discussed in previous sections, the selection of
the ν̄ij parameters is crucial to the proper operation of the
algorithm. Specifically, these parameters should be chosen
such that their corresponding optimal dual variables are upper
bounded by them. In Figure 7 we plot the time averaged dual
variables νij , where the average over time of these variables
converges to their optimal value. First, the choice of ν̄ij = 19
for all i, j satisfies the required condition, i.e., being an upper
bound of the optimal dual variables. Further, by taking a closer
look at these dual variables, we can gain some insight into the
behavior of the system. First, note that the variables νii are
associated to the constraint sii ≤ E qi zi + yii , and hence,
represent the requirement of plant i to transmit its state. By
looking at expression (22), corresponding to the closed-form
solution of the primal zi [t], a larger value of the dual variable
νii [t] leads to a larger value of the scheduling variable zi [t].
Thus, since System 1 is more unstable than System 2, we
have that ν11 > ν22 . In a similar way, the variables νij for
j 6= i are associated to the constraint sij ≥ E qc zj − yij and
represent, at node i, the interference-adjusted need to transmit
of the other plants j 6= i. Therefore, in our two-plant scenario,
a large value of ν21 [t] leads to a lower z1 [t]. And in a similar
manner as previously, since System 1 is more unstable than
System 2, and the systems interfere symmetrically, we have
that ν12 > ν21 . Also, when previously evaluating Fig. 4, we
expected System 1 to have higher battery requirements due
to its higher instability. Now, as per Proposition 3, which
states that the required battery size bmax
is proportional to
i
the dual variables ν̄ii by the inequality bmax
≥ 1 ν̄ii + 1, we
i
see that System 1 requires a larger ν̄ii value, and hence, a

larger battery.

Reception
Required
Transmission
0.8

1
·10

4

C. Communication
We turn our attention to the communication aspects of
the proposed policy. As we have discussed previously, System 1 is slightly more unstable than System 2. Since we
are requiring for both plants an expected decrease rate of
ρi = 0.8 for a Lyapunov function Vi (xi [t]) = x2i [t], by
Proposition 1 this translates to required successful transmission probabilities of p1 ≈ 0.3453 and p2 ≈ 0.2769,
respectively. As expected, the less stable system requires a
higher successful transmission probability. In Figure 8 we
have plotted the resulting transmission probabilities during
our simulation. We look at three different probabilities, (i)
the required transmission probabilities, given P
by pi ; (ii) the
t
actual transmission probabilities, pTi X , (1/t) l=0 zi [l]; and
RX
(iii) the
probabilities,
,
 given by pi
Ptsuccessful reception
Q
(1/t) l=0 qi [l]zi [l] j6=i (1 − qc zj [l]) .
While the required probabilities are p1 ≈ 0.3453 and
p2 ≈ 0.2769, we have that the actual successful reception
probabilities are pRX
= 0.3607 and pRX
= 0.2827. These
1
2
probabilities are over the required ones and, hence, ensure the
stability of the systems. However, the probability at which the
sensors try to access the medium are higher, pT1 X = 0.4446
and pT2 X = 0.3558, respectively. This is due to the effects
of the transmission medium. Packets can be lost if collisions
occur, and they might not be decoded properly if the channel
conditions are not sufficiently favorable (cf. Figure 2).
The overall effect of the transmission medium is better displayed in Figure 9, where we show the transmission schedules
from t = 1050 to t = 1100. In this figure, the bars represent
the probability of successful decoding for a given time slot,
and the stems represent an access to the medium. Further,
collisions are represented by a red dot. From this plot, it is
clear that the sensor node tends to access the medium when the
channel conditions are favorable (i.e., the decoding function
qi [t] takes values closer to 1). Also, this figure evidences that
collisions happen with a sufficiently low chance, since due to

13

Plant 1

Idle

Plant 2

Tx

Idle

1,050

1,060

1,070

1,080

1,090

1,100

Time (t)
Fig. 9. Transmission schedules from t = 1050 to t = 1100. Bars represent
the decoding probability qi [t], stems denote a transmission taking place, and
a red dot denotes the occurrence of a collision.

the independence assumption of the channel states between
sensors, access does not happen at the same time very often.
VII. C ONCLUSIONS
In this work, we have designed a random access communication scheme for energy harvesting sensors in wireless control
systems. We have considered a scenario with multiple plants
sharing a wireless communication medium. Under these conditions, we have shown that the optimal scheduling decision is
to transmit with a certain probability, which is adaptive to the
time-varying channel, battery and plant conditions. In order to
compute the optimal policy, we have provided an algorithm
based on a stochastic dual method. The proposed algorithm is
decentralized, where the sensors only need to share some of
their dual variables. Furthermore, we have provided theoretical
guarantees on the stability of all control loops under the
proposed policy, including the consideration of asynchroniticty
between the information shared between the nodes. Finally, we
have validated by means of simulations the performance of the
proposed scheme. The numerical results show that the random
access policy is capable of stabilizing all control loops while
also satisfying the energy constraints imposed by the energy
harvesting process.
R EFERENCES
[1] M. Calvo-Fullana, C. Antón-Haro, J. Matamoros, and A. Ribeiro,
“Random access policies for wireless networked control systems with
energy harvesting sensors,” in American Control Conference (ACC),
2017. IEEE, 2017, pp. 3042–3047.
[2] R. J. Vullers, R. Schaijk, H. J. Visser, J. Penders, and C. V. Hoof,
“Energy harvesting for autonomous wireless sensor networks,” IEEE
Solid-State Circuits Magazine, vol. 2, no. 2, pp. 29–38, 2010.
[3] J. Yang and S. Ulukus, “Optimal packet scheduling in an energy harvesting communication system,” IEEE Transactions on Communications,
vol. 60, no. 1, pp. 220–230, 2012.
[4] K. Tutuncuoglu and A. Yener, “Optimum transmission policies for
battery limited energy harvesting nodes,” IEEE Transactions on Wireless
Communications, vol. 11, no. 3, pp. 1180–1189, 2012.
[5] O. Ozel, K. Tutuncuoglu, J. Yang, S. Ulukus, and A. Yener, “Transmission with energy harvesting nodes in fading wireless channels: Optimal
policies,” IEEE Journal on Selected Areas in Communications, vol. 29,
no. 8, pp. 1732–1743, 2011.

[6] C. K. Ho and R. Zhang, “Optimal energy allocation for wireless
communications with energy harvesting constraints,” IEEE Transactions
on Signal Processing, vol. 60, no. 9, pp. 4808–4818, 2012.
[7] M. Calvo-Fullana, J. Matamoros, and C. Antón-Haro, “Reconstruction
of correlated sources with energy harvesting constraints in delayconstrained and delay-tolerant communication scenarios,” IEEE Transactions on Wireless Communications, vol. 16, no. 3, pp. 1974–1986,
2017.
[8] P. Castiglione and G. Matz, “Energy-neutral source-channel coding with
battery and memory size constraints,” IEEE Transactions on Communications, vol. 62, no. 4, pp. 1373–1381, 2014.
[9] O. Orhan, D. Gunduz, and E. Erkip, “Source-channel coding under
energy, delay, and buffer constraints,” IEEE Transactions on Wireless
Communications, vol. 14, no. 7, pp. 3836–3849, July 2015.
[10] P. Castiglione, O. Simeone, E. Erkip, and T. Zemen, “Energy management policies for energy-neutral source-channel coding,” IEEE Transactions on Communications, vol. 60, no. 9, pp. 2668–2678, 2012.
[11] G. Yang, V. Y. Tan, C. K. Ho, S. H. Ting, and Y. L. Guan, “Wireless
compressive sensing for energy harvesting sensor nodes,” IEEE Transactions on Signal Processing, vol. 61, no. 18, pp. 4491–4505, 2013.
[12] S. Knorn, S. Dey, A. Ahlén, and D. E. Quevedo, “Distortion minimization in multi-sensor estimation using energy harvesting and energy
sharing.” IEEE Trans. Signal Processing, vol. 63, no. 11, pp. 2848–2863,
2015.
[13] M. Calvo-Fullana, J. Matamoros, and C. Antón-Haro, “Sensor selection
and power allocation strategies for energy harvesting wireless sensor
networks,” IEEE Journal on Selected Areas in Communications, vol. 34,
no. 12, pp. 3685–3695, 2016.
[14] S. Ulukus, A. Yener, E. Erkip, O. Simeone, M. Zorzi, P. Grover, and
K. Huang, “Energy harvesting wireless communications: A review of
recent advances,” IEEE Journal on Selected Areas in Communications,
vol. 33, no. 3, pp. 360–381, 2015.
[15] A. Nayyar, T. Başar, D. Teneketzis, and V. V. Veeravalli, “Optimal
strategies for communication and remote estimation with an energy
harvesting sensor,” IEEE Transactions on Automatic Control, vol. 58,
no. 9, pp. 2246–2260, 2013.
[16] M. Nourian, A. S. Leong, and S. Dey, “Optimal energy allocation for
kalman filtering over packet dropping links with imperfect acknowledgments and energy harvesting constraints,” IEEE Transactions on
Automatic Control, vol. 59, no. 8, pp. 2128–2143, 2014.
[17] Y. Li, F. Zhang, D. E. Quevedo, V. Lau, S. Dey, and L. Shi, “Power
control of an energy harvesting sensor for remote state estimation,” IEEE
Transactions on Automatic Control, vol. 62, no. 1, pp. 277–290, 2017.
[18] J. Huang, D. Shi, and T. Chen, “Event-triggered state estimation with
an energy harvesting sensor,” IEEE Transactions on Automatic Control,
2017.
[19] N. J. Watkins, K. Gatsis, C. Nowzari, and G. J. Pappas, “Battery
management for control systems with energy harvesting sensors,” in
IEEE Conference on Decision and Control. IEEE, 2017.
[20] ——, “Stability of control systems with feedback from energy harvesting
sensors,” arXiv preprint arXiv:1712.02847, 2017.
[21] W. Zhang, M. S. Branicky, and S. M. Phillips, “Stability of networked
control systems,” IEEE Control Systems, vol. 21, no. 1, pp. 84–99, 2001.
[22] D. Hristu-Varsakelis, “Feedback control systems as users of a shared
network: Communication sequences that guarantee stability,” in Decision
and Control, 2001. Proceedings of the 40th IEEE Conference on, vol. 4.
IEEE, 2001, pp. 3631–3636.
[23] G. C. Walsh, H. Ye, and L. G. Bushnell, “Stability analysis of networked
control systems,” IEEE transactions on control systems technology,
vol. 10, no. 3, pp. 438–446, 2002.
[24] K. Gatsis, A. Ribeiro, and G. J. Pappas, “Random access design for
wireless control systems,” arXiv preprint arXiv:1605.00627, 2016.
[25] A. Ribeiro, “Ergodic stochastic optimization algorithms for wireless
communication and networking,” IEEE Transactions on Signal Processing, vol. 58, no. 12, pp. 6369–6386, 2010.
[26] A. Goldsmith, Wireless communications. Cambridge University Press,
2005.
[27] A. Ribeiro, “Optimal resource allocation in wireless communication
and networking,” EURASIP Journal on Wireless Communications and
Networking, vol. 2012, no. 1, p. 272, 2012.
[28] D. P. Bertsekas and J. N. Tsitsiklis, Parallel and distributed computation:
numerical methods. Prentice hall Englewood Cliffs, NJ, 1989, vol. 23.
[29] R. Durrett, Probability: theory and examples. Cambridge university
press, 2010.

