1

A Stochastic Resource-Sharing Network for Electric
Vehicle Charging

arXiv:1711.05561v2 [math.OC] 13 Jan 2018

Angelos Aveklouris, Maria Vlasiou, and Bert Zwart, Member, IEEE

Abstract—We consider a distribution grid used to charge
electric vehicles such that voltage drops stay bounded. We
model this as a class of resource-sharing networks, known
as bandwidth-sharing networks in the communication network
literature. We focus on resource-sharing networks that are driven
by a class of greedy control rules that can be implemented in
a decentralized fashion. For a large number of such control
rules, we can characterize the performance of the system by
a fluid approximation. This leads to a set of dynamic equations
that take into account the stochastic behavior of EVs. We show
that the invariant point of these equations is unique and can be
computed by solving a specific ACOPF problem, which admits an
exact convex relaxation. For the class of weighted proportional
fairness control, we show additional appealing properties under
the linearized Distflow model, such as fairness and a productform property of the stochastic model.
Index Terms—Electric vehicle charging, distribution network,
AC power flow model, linearized Distflow, queueing theory,
stochastic processes, fluid approximation.

I. I NTRODUCTION

T

HE rise of electric vehicles (EVs) is unstoppable due
to factors such as the decreasing cost of batteries and
various policy decisions [1]. These vehicles need to be charged
and will therefore cause congestion in distribution grids in the
(very near) future. This paper proposes to model and analyze
such congestion by the use of a class of resource-sharing networks, which in the queueing network community are known
as bandwidth-sharing networks. Bandwidth-sharing networks
are a specific class of queueing networks, where customers
(in our context: EVs) need to be served simultaneously. Their
service requires the concurrent usage of multiple “servers” (in
our case: all upstream lines between the location of the EV and
the feeder of the distribution grid). Determining how fast to
charge each EV taking into account network stability and the
randomness of future arriving EVs is one of the key problems
in the analysis of distribution grids, leading to challenging
mathematical problems.
Similar questions have been successfully answered in communication networks, where the set of feasible schedules is
The research of Angelos Aveklouris is funded by a TOP grant of the
Netherlands Organization for Scientific Research (NWO) through project
613.001.301. The research of Maria Vlasiou is supported by the NWO MEERVOUD grant 632.003.002. The research of Bert Zwart is partly supported by
the NWO VICI grant 639.033.413.
A. Aveklouris (a.aveklouris@tue.nl) is with the Department of Mathematics
and Computer Science, Eindhoven University of Technology, The Netherlands.
M. Vlasiou (m.vlasiou@tue.nl) and B. Zwart (bert.zwart@cwi.nl) are with
the Centrum Wiskunde en Informatica and the Department of Mathematics and
Computer Science, Eindhoven University of Technology, The Netherlands.

determined by the maximum amount of data a communication channel can carry per time unit, leading to the powerful concept of bandwidth-sharing networks [2]. Bandwidthsharing networks couple the important fields of network utility
maximization with stochastic process dynamics [3, 4]. Apart
from yielding various qualitative insights, they have been
instrumental in the comparison and optimization of various
data network protocols and even in the design of new protocols
[5]. The stochastic analysis of bandwidth-sharing networks
was initially restricted to some specific examples [6, 7]. By
now, fluid and diffusion approximations are available, which
are computationally tractable [8, 9, 10, 11, 12, 13] and hold
for a large class of networks.
From a mathematical viewpoint, the present paper is influenced by [11]. However, in the setting of charging EVs,
an important constraint that needs to be satisfied is to keep
the voltage drops bounded, making the bandwidth-sharing
network proposed in this paper different from the abovementioned works. This also causes new technical issues, as the
capacity set can be non-polyhedral or even non-convex. The
first paper to suggest the class of bandwidth-sharing networks
in the content of EV charging is [14], where simulation studies
were conducted to explore stability properties, assuming that
the arrival times are exponentially distributed. Our work is
a significant extension of [14], both in terms of models and
results: we allow for load limits, finitely many parking spaces,
deadlines (associated with parking times), and do not make
any assumption on the joint distribution of the parking time
and the demand for electricity. More importantly, despite our
assumptions leading to an intricate class of measure-valued
processes, we obtain a number of mathematical results that
are computationally tractable and in some cases explicit.
We develop a fluid approximation for the number of uncharged EVs (for the single-node Markovian case see [15]),
allowing the dynamics of the stochastic model to be approximated with a deterministic model. This model is still quite
rich, as it depends on the joint distribution of the charging
requirements and parking times. We show that the invariant
point of this dynamical system is unique and can be characterized in a computationally friendly manner by formulating a
nontrivial AC optimal-power-flow problem (ACOPF), which
is tractable as its convex relaxation is exact in many cases.
When we replace the AC load flow model with the simpler
linearized Distflow [16], we obtain more explicit results, as
the capacity set becomes polyhedral. For each load flow
model, we compare the fluid approximation with its original
stochastic model and also compare across load flow models.
Our computational study, though certainly not exhaustive,

2

indicates that our fluid models can accurately predict the
performance of a specific control rule, with relative errors of
no more than 10% in most cases.
For the class of weighted proportional fairness controls
where the weights are chosen as function of the line resistances, we derive a fairness property. In this case, all EVs
are charged at the same rate, independent of the location in
the network, while keeping voltage drops bounded. When the
weights are instead chosen equally, we can even derive an
explicit formula for the invariant distribution of the original
stochastic process. Specifically, we show that under certain
assumptions, our network behaves like a multiclass processorsharing queue. Such properties have proven quite fruitful in
other areas of engineering, particularly in the analysis of
computer systems [17], communication networks [6], and
wireless networks [18].
EVs can be charged in several ways. Our setup can be
seen as an example of slow charging, in which drivers typically park their EV and are not physically present during
charging (but are busy shopping, working, sleeping, etc).
For queueing models focusing on fast charging, we refer to
[19, 20]. Both papers consider a gradient scheduler to control
delays, but do not consider physical load flow models as is
done here. [21] presents a queueing model for battery swapping. [22] is an early paper on a queueing analysis of EV
charging, focusing on designing safe (in term of voltage drops)
control rules with minimal communication overhead.
We can only provide a small additional (i.e., non-queueing)
sample of the already vast but still emerging literature on EV
charging. The focus of [23] is on a specific parking lot and
presents an algorithm for optimally managing a large number
of plug-in EVs. Algorithms to minimize the impact of plugin EV charging on the distribution grid are proposed in [24].
In [25], the overall charging demand of plug-in EVs is considered. Mathematical models where vehicles communicate
beforehand with the grid to convey information about their
charging status are studied in [26]. In [27], EVs are the central
object and a dynamic program is formulated that prescribes
how EVs should charge their battery using price signals.
The present paper aims to illustrate how state-of-the-art
methods from the applied probability and queueing communities can contribute to the analysis of the interplay between
EV charging and the analysis of congestion in distribution
networks. Our analysis does not take into account other
important features in distribution networks touched upon in
some of the above-cited works. In particular, it would be useful
to incorporate smart appliances/buildings/meters, rooftop solar
panels, and other sources of electricity demand. We think that
the tractability and generality of our formulation and the established connection with an OPF problem makes our framework
promising towards a comprehensive stochastic network model
of a distribution grid.
The paper is organized as follows. In Section II, we provide
a detailed model description — in particular we introduce our
stochastic model, the class of charging controls, and the load
flow models. The section concludes with a fairness property
of a specific class of controls. Section III shows how the
proportional fairness control mechanism leads to a product-

form property for the stochastic model. In Section IV, we
consider our model in full generality. We present a fluid model
of our system and show that the associated dynamic equations
have a unique invariant point, which is shown to be stable
under an additional assumption. In Section V, we apply our
results to design a control rule that maximizes the fraction
of successful charges. Numerical validations are presented in
Section VI. All proofs are gathered in the Appendix.
II. M ODEL DESCRIPTION
In this section, we provide a detailed formulation of our
model and explain various notational conventions that are used
in the remainder of this work. A notational convention is that
all vectors and matrices are denoted by bold letters.
A. Network and infrastructure
Consider a radial distribution network described by a rooted
tree graph G = (I, E), where I = {0, 1, . . . , I}, denotes
its set of nodes (buses) and E is its set of directed edges,
assuming that node 0 is the root node (feeder). Further, we
denote by ik ∈ E the edge that connects node i to node k,
assuming that i is closer to the root node than k. Let I(k)
and E(k) be the node and edge set of the subtree rooted in
node k ∈ I. The active and reactive power consumed by
the subtree (I(k), E(k)) are PI(k) and QI(k) . The resistance,
the reactance, and the active and reactive power losses along
Q
edge ik are denoted by Rik , Xik , LP
ik , and Lik , respectively.
Moreover, Vi is the voltage at node i and V0 is known. At
any node, except for the root node, there is a charging station
with Ki > 0, i ∈ I \{0} parking spaces (each having an
EV charger). Further, we assume there are J = {1, . . . , J}
different types of EVs indexed by j.
B. Stochastic model for EVs
Type-j EVs arrive at node i according to a Poisson process
with rate λij . If all spaces are occupied, a newly arriving EV
does not enter the system, but is assumed to leave immediately.
Each EV has a random charging requirement and a random
parking time. These depend on the type of the EV, but
are independent between EVs and are denoted by Bj and
Dj , respectively, for type-j EVs. In queueing terminology,
these quantities are called service requirements and deadlines,
respectively. The joint distribution of (Bj , Dj ) is given by a
bivariate probability distribution: Fj (x, y) = P (Bj ≤ x, Dj ≤
y) for x, y ≥ 0.
Our framework is general enough to distinguish between
types. For example, we can classify types according to intervals of ratio of the charging requirement and parking time
and/or according to the contract they have. An EV leaves
the system after its parking time expires. It may be not fully
charged. If an EV finishes its charge, it remains at its parking
space without consuming power until its parking time expires.
EVs that have finished their charge are called “fully charged”.
C. State descriptor
We denote by Qij (t) ∈ {0, 1, . . . , Ki } the total number of
type-j EVs at node i at time t ≥ 0, where Qij (0) is the
PJ
initial number of EVs. Thus, Qi (t) := j=1 Qij (t) denotes

A STOCHASTIC RESOURCE-SHARING NETWORK FOR ELECTRIC VEHICLE CHARGING

the total number of EVs at node i. Further, we denote by
Zij (t) ∈ {0, 1, . . . , Qij (t)} the number of type-j EVs at node
i with a not-fully-charged battery at time t and by Zij (0) the
number of vehicles initially at node i. Last, we write Zij (∞)
or simply Zij to represent the process in steady-state.
For some fixed time, let z = (zij : i ∈ I \{0}, j ∈
J ) ∈ RI×J
be the vector giving the number of uncharged
+
EVs for all types and nodes. Note that although the vector that
gives the number of uncharged EVs should have integer-valued
coordinates, we allow real values in order to accommodate
fluid analogues later. Moreover, we assume that EVs receive
only active power during their charge and do not absorb
reactive power; see [28] for a justification.
D. Charging control rule
An important part of our framework is the way we specify
how the charging of EVs takes place. Given that the state of
uncharged vehicles is equal to z, we assume the existence
of a function p(z) = (pij (z) : i ∈ I \{0}, j ∈ J )
that specifies the instantaneous rate of power each uncharged
vehicle receives. Moreover, we assume that this function is
obtained by optimizing a “global” function. Specifically, for
a type-j EV at node i we associate a function uij (·), which
is strictly increasing and concave in R+ , twice differentiable
in (0, ∞) with limx→0 u0ij (x) = ∞. The charging rate p(z)
PI PJ
is then determined by maxp i=1 j=1 zij uij (pij ) subject
to a number of constraints, which take into account physical
limits on the charging of the batteries, load limits, and most
importantly voltage drop constraints.
Before we describe these constraints in detail, we first
provide some comments on this charging protocol. An important example is the choice uij (pij ) = wij log pij , which
is known as weighted proportional fairness. Note that this
scheme assumes the existence of a virtual agent that is
capable of optimizing the global function. In practice, this
control may be implemented in a decentralized fashion, using
mechanisms such as the alternating direction method of
multipliers [29, 30]. It is even possible to come up with
decentralized allocation schemes that achieve this control if
the functions uij (·) are unknown, which dictates the use of
proportional fairness with a specific choice of the weights
wij . For background, we refer to [3, 5]. A limitation of our
formulation is that it does not take into account the remaining
time until the deadline expires and the remaining charging
requirement. Our multiclass framework allows to at least
partially overcome this, for example by letting the functions
uij (·) depend on the joint distribution of (Bj , Dj ) and by
allowing routing between classes.
We now turn to a discussion of the constraints. We assume
that the highest power that the parking lot at node i can
consume is Mi > 0 and that the maximum power rate which
a type-j EV can be charged is cmax
j . That is,
J
X

zij pij ≤ Mi

and 0 ≤ pij ≤

cmax
.
j

(1)

j=1

We refer to (1) as load constraints. In addition, we impose
voltage drop constraints. These constraints rely on load flow
models. Two of these models we consider are described next.

3

E. AC voltage model
We first consider a minor simplification of the full AC power
flow equations. The angle between voltages in distribution
networks are small and hence they can be chosen so that
the phasors have zero imaginary components [31, Chapter 3].
Under this assumption, Kirchhoff’s law [32, Eq. 1] for our
model takes the form, for pk ∈ E,
Vp Vk − Vk Vk − PI(k) Rpk − QI(k) Xpk = 0,

(2)

where p ∈ I denotes the unique parent of node k. The previous
equations are non-linear. Applying the transformation




Vp2 Vp Vk
Wpp Wpk
W (pk ) =
=:
Wkp Wkk
Vk Vp Vk2
leads to the following linear equations (in terms of W (pk )):
Wpk − Wkk − PI(k) Rpk − QI(k) Xpk = 0, pk ∈ E .

(3)

Note that W (pk ) are positive semidefinite matrices (denoted
by W (pk )  0) of rank one. The active and reactive power
consumed by the subtree (I(k), E(k)) are given by
PI(k) =

J
X X

zlj plj +

l∈I(k) j=1

QI(k) =

X

X

X

X

LP
ls ,

(4)

l∈I(k) ls ∈E(k)

LQ
ls ,

l∈I(k) ls ∈E(k)

where by [14, Appendix B],
2
2
LP
ls = (Wll − 2Wls + Wss )Rls /(Rls + Xls ),
2
2
LQ
ls = (Wll − 2Wls + Wss )Xls /(Rls + Xls ).

Note that Wkk are dependent on vectors p and z. We write
Wkk (p, z), when we wish to emphasize this. If we use this
model to describe voltages, the function p(z) is then given by
maxp

I X
J
X

zij uij (pij )

i=1 j=1

subject to

(5)

(1), (3), υ i ≤ Wii ≤ υ i ,
W (ik )  0, rank(W (ik )) = 1, ik ∈ E,

for zij > 0. If zij = 0, then take pij = 0. In addition,
0 < υ k ≤ W00 ≤ υ k are the voltage limits. Observe that
the optimization problem (5) is non-convex and NP hard due
to rank-one constraints. Removing the non-convex constraints
yields a convex relaxation, which is a second-order cone
program. Further, by Remark 1 (see below) and [33, Theorem
5], we obtain that the convex-relaxation problem is exact.
F. Linearized Distflow model
Though the previous voltage model is tractable enough for a
convex relaxation to be exact, it is rather complicated. Assuming that the active and reactive power losses on edges are small
relative to the power flows, but now allowing the voltages
to be complex numbers, we derive a linear approximation
of the previous model, called the linearized (or simplified)

4

Distflow model [16]. In this case, the voltage magnitudes
lin
Wkk
:= |Vklin |2 have an analytic expression [32, Lemma 12]:
lin
Wkk
(p, z) = W00 − 2

X
ls ∈P(k)

Rls

J
X X

zmj pmj , (6)

m∈I(s) j=1

where the P(k) is the unique path from the feeder to node k.
lin
Remark 1. Note that Wkk
≤ W00 for all nodes k, as
we assume that the nodes only consume power, and by [32,
lin
Lemma 12] we obtain Wkk (p, z) ≤ Wkk
(p, z).

To derive the representation of the power allocation mechanism p(z) in this setting, one replaces the constraints in (5)
lin
(p, z) ≤ υ k .
by (1) and υ k ≤ Wkk
When adding stochastic dynamics, the resulting model is
still rather complicated. Even Markovian assumptions yield a
multidimensional Markov chain of which the transition rates
are governed by solutions of nonlinear programming problems.
In the following proposition, we identify a special case of our
setting in which p(z) is explicit and leads to a fair allocation
of power to all users, in the sense that the charging rate of
a battery does not depend on the location where it is parked.
Specifically, if we remove the load constraints by setting Mi =
∞ and cmax
= ∞, we can not directly control the power each
j
node receives after the functions uij (·) have been set. As a
consequence, this may lead to the situation that the nodes close
to the feeder consume almost all the power.
An important question is under what assumptions all EVs in
the system are charged at the same power rate. The following
proposition
P gives a partial positive answer to this question. Let
Ri = ls ∈P(i) Rls and δ = W002−υ .
Proposition 1 (Fairness property). Let υ i = υ and the network
be a line. Moreover, assume that Mi = ∞ and cmax
= ∞.
j
Take uij (pij ) = wij u(pij ) and the power flow model (6). If
wi = Ri , then we have that pij (z) = p(z) > 0. Moreover,
P
−1
PJ
I
p(z) = δ
.
i=1 Ri
j=1 zij
A similar result can be shown for more general trees, under
the assumption that the root node has only one child.
III. A N ALLOCATION MECHANISM WITH THE
PRODUCT- FORM PROPERTY
The dynamics of the high-dimensional non-Markovian
stochastic process (Qij (t), Zij (t)), t ≥ 0, are in general
not tractable from a probabilistic viewpoint. To obtain a
Markovian description, we would also have to keep track of all
residual parking times and charging requirements, leading to
a measure-valued process as in [11]. Therefore, we consider
fluid approximations of (Qij (t), Zij (t)) in the next section,
which are more tractable.
Despite the complexity of our stochastic model, we are
able to identify a special case for which the entire network
behaves like a multiclass processor-sharing queue, of which
the invariant distribution is explicit, and for which even timedependent properties are known [4, 34].
We take J = 1 for convenience and drop all indices
j from the notation in this
PIsection. For every node i, let
ρi = λi E[B]Ri /δ and ρ = i=1 ρi .

Theorem 1. Assume Ki = Mi = ∞ for all i, and a
line network. Assume also cmax = ∞. Furthermore, consider
the power allocation rule p(z) under the linearized Distflow
model (6). If the power allocation rule is proportional fairness,
i.e., ui (pi ) = log pi , then, for every n ∈ NI+ ,
I
I
X
Y

ρni i
,
lim P Z(t) = n = (1 − ρ)(
ni )!
t→∞
n!
i=1
i=1 i

(7)

provided ρ < 1.
Note that this result is valid for arbitrary distributions of the
charging requirements and as such it provides an insensitivity
property. This result is another exhibition of the appealing
nature of proportional fairness, which has also been shown to
give similar nice properties in communication network models
[4, 6, 7, 18].
The proof of this theorem follows from a similar argument
as in [6], making a connection with the class of Whittle
networks, by showing that a specific local balance property
called balanced fairness is satisfied. We explain this procedure
for the case of exponential charging times in the Appendix.
IV. F LUID A PPROXIMATION
In this section, we develop a fluid approximation for the
stochastic model defined in Section II, of which the invariant
point is characterized through an OPF problem. To do so, we
follow a similar approach as in [11] and [35].
The fluid approximation, which is deterministic, can be
thought of as a formal law of large numbers approximation.
More precisely, consider a family of models as defined in
Section II, indexed by
a scaling parameter n ∈ N. The fluid
n
Zij
(·)
scaling is given by n . To obtain a non-trivial fluid limit, we
choose the following scaling for the node parameters in the nth
system. The maximum power that node i can consume is nMi ,
the arrival rate is nλij , the number of parking spaces is nKi ;
all other parameters remain unchanged. A mathematically
rigorous justification of this scaling is beyond the scope of
this study, and will be pursued elsewhere. If the set of feasible
power allocations is polyhedral, the methods from [11] can be
applied directly to achieve this justification. This scaling gives
rise to the following definition of a fluid model.
Definition 1 (Fluid model). A nonnegative continuous vectorvalued function z(·) is a fluid-model solution if it satisfies the
functional equations for i, j ≥ 1
Z t

zij (t) = zij (0)P Bj0 >
pij (z(u))du, Dj0 ≥ t
0
Z t
Z t

+
γij (s)P Bj >
pij (z(u))du, Dj ≥ t − s ds,
0
s
P
where γij (t) := λij 1{qi (t)<Ki } , and qi (t) = j qij (t), with
Z t


0
γij (s)P Dj ≥ t − s ds.
qij (t) = qij (0)P Dj ≥ t +
0

Bj0

Dj0

Further,
and
are the charging requirement and the
parking time for the initial population in the system.
The time-dependent fluid model solution can be used directly to approximate the evolution of the system at time t;

A STOCHASTIC RESOURCE-SHARING NETWORK FOR ELECTRIC VEHICLE CHARGING

e.g., one may take E[Zij (t)] ≈ zij (t). This set of equations can
be extended to time-varying arrival rates by replacing λij by
λij (t). Also, one can consider schemes in which blocked EVs
are not lost, but routed to adjacent parking lots, which leads
to further modifications to γij (·). The fluid model equations,
though still rather complicated,
have an intuitive
 meaning; e.g.,
Rt
the term P Bj > s pij (z(u))du, Dj ≥ t − s resembles the
fraction of EVs of type-j admitted to the system at time s
at node i that are still in the system at time t (for this to
happen, their deadline needs to exceed t − s and their service
requirement
R t needs to be bigger than the service allocated,
which is s pij (z(u))du). To solve the fluid model equations,
one can proceed numerically by Picard iteration. In some
cases, the set of equations can also be solved explicitly, as
the next proposition illustrates.
Proposition 2. Assume that Ki = ∞, J = 1, and the
charging requirements and parking times are exponential and
independent. Under the assumptions of Proposition 1 with
u(pi ) = log(pi ), we have that
zi (t) = zi∗ + (zi (0) − zi∗ )e−t/E[D] ,
with
zi∗ = E[D] (λi − Λ∗i /E[B]) , Λ∗i = PI

λi δ

k=1

Rk λk

.

The previous proposition continues to hold for Ki big
enough such that the parking lots are never full.
We now turn to the behavior of our fluid model in equilibrium, i.e., for t = ∞. In this case, we obtain a computationally
tractable characterization through a particular OPF problem.
Before we state our main theorem, we introduce some notation. Let

−1
J
n
o
X
0
λij
λ
ij
γij :=
min λi , Ki 
E[Dj 0 ]
, (8)
λi
λi
0
j =1

where λi :=

PJ

j=1

λij . Furthermore, define the functions

gij (x) := γij E[min{Dj x, Bj }],

(9)

and the node allocation (the power which type-j EVs consume
at node i), Λij (z) := zij pij (z). Also, for a random variable
Y , denote by inf(Y ) the leftmost point of its support.
∗
Theorem 2 (Characterization of invariant point). (i) If zij
is
an invariant point for the fluid model, it is given by the solution
of the fixed-point equation
∗
zij
= γij E[min{Dj ,

Bj
}].
pij (z ∗ )

(10)

∗
(ii) Let inf (Dj /Bj ) ≤ 1/cmax
j . The solution z of (10) is
Λ∗

∗
ij
∗
unique and is given by zij
= g−1 (Λ
∗ ) , where Λ is the unique
ij
ij
solution of the optimization problem

maxΛ

I X
J
X

Gij (Λij )

i=1 j=1

subject to Wik − Wkk − PI(k) Rik − QI(k) Xik = 0, (11)
υ i ≤ Wii ≤ υ i , Λij ≤ Mi ,
0 ≤ Λij ≤ gij (cmax
), ik ∈ E .
j

5

Gij (·) is a strictly concave function such that G0ij (·) =
−1
u0ij (gij
(·)).
By (4), observe that Wkk depends on z through the products
zij pij (z). By the definition of Λ, we have that Wkk depends
only on Λ. That is, the previous optimization problem is
indeed independent of the fixed point z ∗ .
The intuition behind (10) follows from a result in queueing
theory known as Little’s law [17], which says that the expected
number of customers equals the arrival rate times the expected
sojourn time of a particular customer. The sojourn time of a
customer of type-j at node i is the minimum of its deadline and
its potential service time. The latter quantity is approximated
by the total service time divided by the service rate, which will
be approximately constant in equilibrium for large systems
(this is called the snapshot principle in queueing theory [36]).
To arrive at (11) the essential idea is to add Little’s law (10) to
the set of KKT conditions that characterize p(z) and rewrite
all equations in such a way that they form the KKT conditions
for the problem (11). We refer to the Appendix for more
details.
To show that a fluid model solution converges towards the
equilibrium point as t → ∞, we also assume that the power
rate that an EV gets does not increase if a newly arriving EV
enters the system.
Definition 2. An allocation mechanism is called “monotone”
if 0 < y ≤ z implies that pij (y) ≥ pij (z).
By adapting techniques from [10], it can be shown that
this assumption holds in the linearized Distflow model. We
conjecture it is true for the more general AC allocation
mechanism as well.
Proposition 3. If the allocation mechanism is monotone and
∗
Ki = ∞, we have that zij (t) → zij
, as t → ∞.
We conjecture this result to be true for general AC allocation
mechanisms, without imposing the condition Ki = ∞, but the
mathematical techniques required to establish this are beyond
the scope of this paper.
V. D ESIGNING CONTROL SCHEMES
A natural question is how to choose suitable functions uij (·)
in order to optimize the behavior of the system. In this section,
we illustrate the scope of this question by focusing on a
special case. In particular, we consider the weighted proportional fairness allocation, line network topology, the linearized
Distflow model, a single type of EVs, and Mi = cmax = ∞.
The main goal is to describe a procedure of choosing the
weights wi , i ≥ 1, such that the log-run fraction of EVs that
gets successfully charged is maximized. This is a particularly
relevant problem in anPoverloaded regime and to this end we
I
make the assumption i=1 E[B]γi Ri > δ.


wi
The solution of (11) is given by Λ∗i = gi hR
, where
i
gi (·) is given by (9). Under the assumptions of this section,
(11) has only a single constraint.
P Its associated
P Lagrange multiplier h > 0 satisfies δ = i Ri Λ∗i = i Ri gi (wi /(hRi )).
hRi Λ∗
i
It follows by Theorem 2 that zi∗ =
wi . The probability an EV gets successfully charged at node i is given

6

by P(D > Bzi∗ /Λ∗i ) = P(wi D > hRi B). Note that the
last probability depends on the Langrange multiplier which
depends on the weights in an intricate fashion. To overcome
this difficulty we scale w such that h(w) = 1. This can always
be achieved as multiplying all weights with the same factor
scales h accordingly but leads to the same solution Λ∗ . Thus,
we can formulate the problem of optimally selecting weights
as follows:
I
X
maxw
γi P(wi D > Ri B)
subject to

i=1
I
X

(12)
γi E[min{wi D, Ri B}] ≤ δ.

i=1

This problem can be transformed into a recourse allocation
problem as considered in [37]. It is in general non-convex
and its solution depends on the joint distribution of B and D.
In the sequel, consider the special case B = HD, with H
a random variable independent of D. Below, we study some
special cases of H.
A. H is random with decreasing hazard rate
Here, we assume that H is a random variable where its
distribution has decreasing hazard rate. The latter ensures
that transformation of (12) leads to a convex optimization
problem. A well-known example of such a distribution is a
Pareto distribution with parameters a > 1 and κ > 0, i.e.,
κ a
P(H > x) = ( x+κ
) , x ≥ 0. Note that E[H] = κ/(a − 1).
Setting ci = wi /Ri , (12) reduces to
I
X

maxc

γi [1 −

i=1

subject to

κ a
]
ci + κ

I
X
E[D]κγi Ri
i=1

1−a

κ a−1
− 1] ≤ δ.
ci + κ

[

After setting yi = ( ciκ+κ )a−1 , we obtain ci =
I
X

maxy
subject to

i=1
I
X
i=1

a/(a−1)

γi (1 − yi

κ
1/(a−1)

yi

− κ and

)
0 ≤ yi ≤ 1,

I
X

maxc
subject to

i=1
I
X

γi 1{H≤ci }
E[D]γi Ri min{ci , H} ≤ δ.

i=1

There exists a non-empty set I ∗ ⊆ I \{0} such that ci < H
for each i ∈ I ∗ . To see this observe that if I ∗ = ∅, then

i=1

subject to

i=1
I
X

γi xi
(13)
E[D]γi Ri Hxi ≤ δ,

xi ∈ {0, 1},

i=1

which is the well-known knapsack problem.
We now show that the case of deterministic H is the worstcase situation in overload by showing that the value of (13) is
a lower bound for the value of (12). Assume
PI w.l.o.g. E[H] = 1
and use the Markov inequality to obtain i=1 γi P(H < ci ) ≥
PI
In addition, an upper bound for
i=1 γi (1 − 1/ci )1{ci >1} . P
I
the constraint is given by
i=1 E[D]γi Ri E[min{H, ci }] ≤
PI
E[D]γ
R
1
due
to
the fact that E[min{H, ci }] =
i i {ci >0}
i=1
0 if i ∈ I ∗ and E[min{H, ci }] ≤ E[H] = 1, otherwise. These
bounds leads to the problem
maxc
subject to

I
X
i=1
I
X

γi (1 − 1/ci )1{ci >1}
E[D]γi Ri 1{ci >0} ≤ δ.

i=1

It is clear that the feasible set of the last problem is a subset
of the feasible set of (12). Taking now ci = ∞ if ci > 0 yields
(13) with H = 1.
The numerical results provided in the next section confirm
that the case of deterministic H indeed leads to a lower
fraction of successful charges than the case of a random H.
This in itself results in a selection of weights where nodes
far away from the feeder will not receive any power if H is
deterministic, while this is not the case when H is random.

E[D]γi Ri min{ci , H} =

In this section, we give a numerical illustration of our
results.
A. Stochastic model and fluid approximation

B. H is deterministic
Let H be deterministic and strictly positive. In this case,
setting ci = wi /Ri , (12) reduces to

I
X

I
X

maxx

VI. N UMERICAL EVALUATION
E[D]κγi Ri
(yi − 1) ≤ δ,
1−a

which is a convex problem.

δ≥

which contradicts the overload assumption. Moreover, observe
that the nodes in I ∗ do not increase the objective function and
so we can set ci = 0 for i ∈ I ∗ . On the other hand, if i ∈
/ I ∗,
then ci ∈ [H, ∞) and minimization of the constraint leads to
ci = H. So, the last optimization problem is equivalent to

I
X
i=1

E[B]γi Ri ,

First, we examine the stochastic model and the fluid approximation. We consider the case J = 1, and assume I = 2
nodes on a line. Fix V0 = 1 in the per-unit system, take
R01 = X01 = 0.01, R12 = X12 = 0.005, and υ 1 =
υ 2 = (0.9)2 . We consider a weighted proportional fairness
allocation and weights as chosen in Proposition 1. Further, the
charging requirements and the parking times are exponentially
distributed with unit mean, and we take K1 = K2 = K in all
our computations.
As a first experiment, we consider this model without load
constraints, i.e., Mi = ∞ and cmax = ∞. We take λi = 1.2Ki
and compare E[Zi ] between the AC and Distflow allocation
mechanisms, as well as their fluid approximations zi∗ .
First, we compare the two power flow models for the
stochastic model computing the expected number of uncharged

TABLE I
K
10
20

E[Z1 ], E[Z2 ] (Distflow)
4.5336, 4.6179
14.0174, 14.0385

E[Z1 ], E[Z2 ] (AC)
4.6801, 4.6924
14.1725, 14.1948

Rel. error
3.13 %, 1.58 %
1.09 %, 1.10 %

7

15

Relative error %

EVs at any node; see Table I. We note that the expectations
for the AC model are higher as this model takes into account
power losses leading to lower service rates.

Relative error %

A STOCHASTIC RESOURCE-SHARING NETWORK FOR ELECTRIC VEHICLE CHARGING

10
5
0
10

8

6

M2

4

2

0

0

2

4

6

M1

8

10

15
10
5
0
10

8

6

4

2

M2

0

0

2

4

6

8

10

M1

Fig. 2. AC in case K = (10, 10) and λ = (12, 12)

Next, we evaluate the fluid approximation for the two
load models; see Table II. Observe again that the number of
uncharged EVs for the AC model is higher. The relative error
between the two load flow models is similar to what we saw
in Table I.
TABLE II
z1∗ , z2∗ (Distflow)
4.5769, 4.5769
14.0300, 14.0300
23.6820, 23.6820
33.4293, 33.4293
43.2330, 43.2330

K
10
20
30
40
50

z1∗ , z2∗ (AC)
4.7356, 4.7513
14.1849, 14.2069
23.8357, 23.8597
33.5823, 33.6073
43.3857, 43.4112

Rel. error
3.35%, 3.67%
1.09%, 1.25%
0.64%, 0.74%
0.45%, 0.53%
0.35%, 0.41%

Not surprisingly, the original stochastic model is not numerically tractable for high values of K. Our results though show
that this is not a problem, as our fluid approximation performs
very well:
TABLE III
Relative error of fluid approximations
Distflow
AC
0.95%, 0.86%
1.18%, 1.25%
0.09%, 0.06%
0.08%, 0.08%

K
10
20

15

Relative error %

Relative error %

Now, we move to the case that Mi < ∞ and cmax < ∞. In
this case, we fix cmax = 1 and we plot the relative error between the fluid approximation and the stochastic model for the
fraction of EVs that get successfully charged. Using arguments
from queueing theory, it can be shown that this fraction equals
1 − E[Zi ]/E[Qi ]. For the fluid approximation, we replace the
numerator with zi∗ , while E[Qi ] can be computed explicitly,
using the Erlang Loss formula [4]. Figures 1 and 2 show the
results for all possible values of Mi and for both load models.
Though the quality of the fluid approximation deteriorates, the
relative error is generally below 10% and for reasonably high
values of Mi , it is even smaller. We also expect these results
to improve for bigger K. For higher values of K one needs
to solve millions of OPF problems of type (5) to obtain the
steady-state behavior of the stochastic model, while the fluid
model only requires the solution of a single OPF problem.

10
5
0
10

8

6

M2

4

2

0

0

2

4

6

M1

8

10

15
10
5
0
10

8

6

M2

4

2

0

0

Fig. 1. Distflow in case K = (10, 10) and λ = (12, 12)

2

4

6

M1

8

10

B. Optimal weights
Now, we move to the choice of optimal weights; see Section V. We take ten nodes and parameters of the system such
P10
i Ri
= 1.17; i.e., the system is 17% overloaded.
that i=1 E[B]γ
δ
We present three examples. First, for H deterministic and
equal to 1. Second, for H being Pareto distributed with
E[H] = 1 and two shape parameter a, namely a = 1.1 and
a = 3.
For the deterministic case, we have that the optimal weights
are (0.0610, 0.0915, 0.1110, 0.1268, 0.1390, 0.1488, 0.1573,
0.1646, 0, 0), and the number of successful charges per unit
of time equals 8.
For the Pareto case, the optimal weights wi , for a = 1.1
are given by the vector (0.1105, 0.1031, 0.1006, 0.0994,
0.0987, 0.0981, 0.0978, 0.0975, 0.0972, 0.0971), and for
a = 3, (0.1293, 0.1173, 0.1093, 0.1033, 0.0985, 0.0945,
0.0911, 0.0881, 0.0854, 0.0830). The number of successful
charges per unit of time equals 10 and 9.49, respectively.
The numerical experiments agree with the results in Section V. It turns out that the Pareto case works better and
also the more variability (smaller a), the better the system
performance. For H being Pareto distributed, the weights are
decreasing and all the nodes receive strictly positive power in
contrast to the deterministic case. In the latter case, the weights
are increasing and some nodes do not receive power. This can
be explained intuitively as follows: the nodes away from the
feeder receive more power in order the successful charges in
the system being maximized. The last statement also agrees
with the choice of the weights in Proposition 1.
VII. C ONCLUDING REMARKS
This paper has proposed a queueing network model for
electric vehicle charging. The main result is a fluid approximation of the number of uncharged vehicles in the system,
which is derived by combining key ideas from network utility
maximization (using a class of utility-maximizing scheduling disciplines), from queueing theory (Little’s law and the
snapshot principle), and from power systems engineering
(AC and linearized Distflow load flow models). Our fluid
approximation explicitly captures the interaction between these
three elements, as well as physical network parameters and
can be computed using convex programming techniques. Our
approach can easily be extended to impose other reliability
constraints, such as line limits, and our fluid approach can be
extended to deal with superchargers, by allowing additional
queues at each station.

8

We focused on the specific class of weighted proportional
fairness protocols. Our optimization framework allows for
a further comparison between charging protocols, which is
a natural next step for further research. As mentioned in
the introduction, another important problem is to extend our
model to allow for batteries to discharge and to include other
features, such as smart appliances/buildings/meters, rooftop
solar panels, and other sources of electricity demand. This
naturally leads to various questions about pricing schemes.
Again, we think that our characterization of the performance of
the system for a fixed control in terms of the OPF problem (11)
can be a useful starting point to design economic mechanisms.
A full mathematical examination of our system is beyond
the scope of the present work. In a follow-up work, we will
∗
rigorously show that zij
is a good approximation of Zij (·)
by developing a fluid limit theorem using the framework of
measure-valued processes, extending the framework of [11].
A PPENDIX
P ROOF OF P ROPOSITION 1
Under the assumptions of Proposition 1, and observing that
Wiilin is decreasing in i, (5) takes the following form:
I X
J
X

maxp

I
X
k=1

k=1

Rk−1k

I X
J
X

zmj pmj ≤ δ.

m=k j=1

By the Karush-Kuhn-Tucker (KKT) conditions, there exists an
h ∈ R+ such that for any i ∈ I \{0} and j ∈ E,

and

zij wij u0 (pij ) = hzij Ri ,

(14)

I
I X
J
X

X
h
Rk−1k
zmj pmj − δ = 0.

(15)

k=1

P ROOF OF P ROPOSITION 2
For simplicity take E[B] = 1/µ and E[D] = 1/ν, and
zi (0) = 0. By the assumption that Ki = ∞, it follows that
γi (t) ≡ γi ≡ λi (for the definition of γi , see (8)). Before we
move to the main part of the proof, we show some helpful
relations.
Under the assumptions of Proposition 1, we can solve (11)
−µΛi
. Writing
explicitly. To see this, note that G0i (Λi ) = wi γiνΛ
i
the KKT conditions for (11), we have that there exists an
h ∈ R+ , such that for
any i ∈ I \{0}, the optimal solution
γ −µΛ∗
Λ∗ satisfies wi iνΛ∗ i = hRi , and
i
I
X
Rk Λ∗k = δ.
(17)

zij wij u(pij )

i=1 j=1

subject to

i
, the transition rates become the
if zi > 0. Setting µi = E[B]R
δ
same as those of a multiclass Markovian processor-sharing
queue with I classes of customers and mean service times µi .
By [39], the stationary distribution of a multiclass processorsharing queue is given by (7).
Further, it is shown in [38], that the stationary distribution
of a multiclass processor-sharing queue is insensitive to the
distribution of service times, if the balanced property satisfied.
That is, (7) holds for general charging requirements B. The
construction in [38] can also be carried out in our setting.

m=k j=1

Note that0 h cannot be zero due to (14). Again, by (14), we have
w u (p )
that ij R ij should be constant for any i, j. In particular,
i

wij u0 (pij )
w11 u0 (p11 )
=
.
Ri
R1

(16)

Choosing wij = Ri , and noting that u0 (·) is a strictly
decreasing function, we have that pij = p11 , for all i, j ≥ 1.
Combining (15) and (16), we derive the expression for p(z).
P ROOF OF T HEOREM 1
Using similar analysis as in the proof of Proposition 1, it
can be shown that pi (z) = R PδI z . Next, observe that the
i
l=1 l
so-called balanced property [38] holds: for i, k ∈ I \{0},
pi (z + ek )pk (z) = pi (z)pk (z + ei ),
where ei , i = 1, . . . , I, denote unit vectors.
If B is an exponential random variable, the process
Z(t), t ≥ 0 is Markov. Let Z(t) = z denote the state of
the process Z(t) at time t ≥ 0. The process can move to state
z + ei with rate λi and to state z − ei with rate zi pi (z)/E[B],

The first
equation∗ (taking into account wi = Ri ) implies that
γi −µΛ∗
γ1 −µΛ1
i
=
. Finally, we get the following expressions
νΛ∗
νΛ∗
1
i
∗
for Λi ,
γi
Λ∗i = Λ∗1 ,
(18)
γ1
and by (17),
γi δ
.
(19)
Λ∗i = PI
k=1 Rk γk
Next, by (10), we compute the invariant point, which is
given by
γi − µΛ∗i
.
(20)
zi∗ =
ν
It is helpful to note the following relation for the invariant
point. Combining the last equation and (18), we get
γi − µ(γi /γ1 )Λ∗1
γi γ1 − µΛ∗1
γi
zi∗ =
=
= z1∗ .
(21)
ν
γ1
ν
γ1
Now, we move to the main part of the proof. Under the
Markovian assumptions, zi (·) is given (alternatively) by the
following ODE:
δ
.
zi0 (t) = γi − νzi (t) − µzi (t) PI
k=1 Rk zk (t)
The last ODE has a unique solution z(·) for given initial point.
So, it is enough to show that the function zi (t) = zi∗ (1−e−νt )
satisfies the previous ODE. Plugging in it into the ODE, we
have that
µz ∗ (1 − e−νt )δ
νzi∗ e−νt = γi − νzi∗ (1 − e−νt ) − PI i
,
∗
−νt )
k=1 Rk zk (1 − e
which can be simplified to γi = νzi∗ +
we derive Λ∗i =
PI

γi δ
,
R k γk

k=1

z∗ δ
PI i
∗.
k=1 Rk zk

µzi∗ δ
∗.
k=1 Rk zk

PI

By (20),

Now, we apply (21) to get Λ∗i =

which holds by (19).

A STOCHASTIC RESOURCE-SHARING NETWORK FOR ELECTRIC VEHICLE CHARGING

P ROOF OF T HEOREM 2
First, we observe some properties of the function gij (·).
Define aij = inf{x ≥ 0 : gij (x) = γij /E[Dj ]}. By [11,
Lemma 6], we have that gij (·) is continuous, strictly increasing
in [0, aij ], and constant in (aij , ∞]. Further, by [11, Lemma
7], we have that if aij < ∞, then inf (Dj /Bj ) = 1/aij and
if aij = 0, then inf (Dj /Bj ) = 0.
Now, we move to the main part of the proof. We note that
for i ∈ I \{0}, we have that
∂
∂
Wkk (z, p) = zij
Wkk (Λ).
∂pij
∂Λij

(22)

To see this, observe that by (2) and the definition of the
bandwidth allocation Λ, as Vk (z, p) depends on the vectors
z and p only through the product zij pij , we get
Vk (z, p) = Vk (Λ).
Applying the transformation Λij (z, p) = zij pij , we have
∂
∂
∂
∂
Vk (z, p) =
Vk (Λ)
Λij = zij
Vk (Λ).
∂pij
∂Λij
∂pij
∂Λij
Recalling that Wkk = Vk Vk , we obtain that
∂
∂
Wkk (z, p) =2Vk (z, p)
Vk (z, p)
∂pij
∂pij
∂
∂
Vk (Λ) = zij
Wkk (Λ).
=2zij Vk (Λ)
∂Λij
∂Λij
By the KKT conditions for (5) there exist Lagrangian
2I+2(I×J)
multipliers (h1i , h2i , h3ij , h̄4ij ) ∈ R+
such that for
i ∈ I \{0}, j ∈ J,
I

∗
zij

X ∂
∂
uij (pij (z ∗ )) =
Wkk (p(z ∗ ))(h1k − h2k )
∂pij
∂pij
k=1

∗
+h3ij zij
+ h̄4ij ,

and

9

observations and (22), the above equations can be rewritten as
follows
I

X ∂
∂
−1
uij (gij
(Λij (z ∗ ))) =
Wkk (Λ)(h1k − h2k )
∂Λij
∂Λij
k=1

+h3ij + h4ij ,
and
h1i (Wii (Λ(z ∗ )) − υ i ) = 0,
h3ij (Λij (z ∗ ) − Mi ) = 0,

h2i (Wii (p(z ∗ )) − υ i ) = 0,

∗
h3ij (zij
pij (z ∗ ) − Mi ) = 0,

h̄4ij (pij (z ∗ ) − cmax
) = 0.
j

∗
Setting h4ij = h̄4ij /zij
, the previous equations take the following (equivalent) form
I
∂
1 X ∂
∗
uij (pij (z )) = ∗
Wkk (p(z ∗ ))(h1k − h2k )
∂pij
zij
∂pij
k=1

+h3ij + h4ij ,
and
h1i (Wii (p(z ∗ )) − υ i ) = 0,

h2i (Wii (p(z ∗ )) − υ i )) = 0,

∗
h3ij (zij
pij (z ∗ ) − Mi ) = 0,

h4ij (pij (z ∗ ) − cmax
) = 0.
j

By the definition of gij (·) and (10), we have that Λij (z ∗ ) =
gij (pij (z ∗ )). Moreover, by the assumption in Theorem 2, we
max
have that 1/aij = inf (Dj /Bj ) ≤ 1/cmax
≤
j . That is, cj
max
aij . Thus, gij (·) is strictly increasing in [0, cj ]. This implies
−1
pij (z ∗ ) = gij
(Λij (z ∗ )). We note that (pij (z ∗ ) − cmax
) =
j
0 if and only if (Λij (z ∗ ) − gij (cmax
))
=
0.
Using
the
last
i

h4ij (Λij (z ∗ ) − gij (cmax
)) = 0.
i

Now, we observe that the last equations are the KKT conditions for (11). To complete the proof, it remains to be
shown that the function Gij (·) is strictly concave. To this end,
−1
observe that gij
(·) is strictly increasing and u0ij (·) is strictly
decreasing since uij (·) is a strictly concave function. It follows
that G0ij (·) is strictly decreasing and hence Gij (·) is a strictly
concave function. The latter implies that the optimization
problem (11) has a unique solution Λ∗ , which is independent
of the invariant point z ∗ . Further, the unique invariant point
is given by
Λ∗ij
Λ∗ij
∗
.
=
zij
=
−1
pij (z ∗ )
gij
(Λ∗ij )
P ROOF OF P ROPOSITION 3
As there similarities with [11, Theorem 3, Corollary 1],
we present only a short proof sketch. Under the assumption
that Ki = ∞, we have that 1{qi (t)<Ki } = 1. By following the same argument as in [11, Theorem 3], there exist
l, h ∈ (0, ∞)I×J , such that
lij ≤ lim inf zij (t) ≤ lim sup zij (t) ≤ hij .
t→∞

(23)

t→∞

Further, using the assumption of the monotone network, l and
h satisfy the following relations
lij = γij E[min{Dj ,

h1i (Wii (p(z ∗ )) − υ i ) = 0,

h2i (Wii (Λ(z ∗ )) − υ i ) = 0,

Bj
Bj
}], hij = γij E[min{Dj ,
}].
pij(l)
pij(h)

By the uniqueness of the invariant point, we have that lij =
hij . Now, by (23), the result follows.
R EFERENCES
[1] “International energy agency. Global EV outlook 2017.” [Online].
Available: https://www.iea.org/publications/freepublications/publication/
GlobalEVOutlook2017.pdf/
[2] L. Massoulié and J. Roberts, “Bandwidth sharing: objectives and algorithms,” in Proc. INFOCOM’99., vol. 3. IEEE, 1999, pp. 1395–1403.
[3] F. Kelly, “Charging and rate control for elastic traffic,” Trans. Emerging
Telecommun. Technol., vol. 8, no. 1, pp. 33–37, 1997.
[4] F. Kelly and E. Yudovina, Stochastic networks. Cambridge University
Press, 2014, vol. 2.
[5] Y. Yi and M. Chiang, “Stochastic network utility maximisationa tribute
to kelly’s paper published in this journal a decade ago,” Trans. Emerging
Telecommun. Technol., vol. 19, no. 4, pp. 421–442, 2008.
[6] T. Bonald and A. Proutiere, “Insensitive bandwidth sharing in data
networks,” Queueing Systems, vol. 44, no. 1, pp. 69–100, 2003.
[7] T. Bonald, L. Massoulié, A. Proutiere, and J. Virtamo, “A queueing analysis of max-min fairness, proportional fairness and balanced fairness,”
Queueing Systems, vol. 53, no. 1, pp. 65–84, 2006.
[8] W. N. Kang, F. P. Kelly, N. H. Lee, and R. J. Williams, “State space
collapse and diffusion approximation for a network operating under a fair
bandwidth sharing policy,” The Annals of Applied Probability, vol. 19,
no. 5, pp. 1719–1780, 2009.

10

[9] H.-Q. Ye and D. D. Yao, “A stochastic network under proportional fair
resource control-diffusion limit with multiple bottlenecks,” Operations
Research, vol. 60, no. 3, pp. 716–738, 2012.
[10] S. Borst, R. Egorova, and B. Zwart, “Fluid limits for bandwidth-sharing
networks in overload,” Mathematics of Operations Research, vol. 39,
no. 2, pp. 533–560, 2014.
[11] M. Remerova, J. Reed, and B. Zwart, “Fluid limits for bandwidth-sharing
networks with rate constraints,” Mathematics of Operations Research,
vol. 39, no. 3, pp. 746–774, 2014.
[12] J. Reed and B. Zwart, “Limit theorems for markovian bandwidth sharing
networks with rate constraints,” Operations Research, vol. 62, no. 6, pp.
1453–1466, 2014.
[13] M. Vlasiou, J. Zhang, and B. Zwart, “Insensitivity of proportional
fairness in critically loaded bandwidth sharing networks,” preprint
arXiv:1411.4841, 2014.
[14] R. Carvalho, L. Buzna, R. Gibbens, and F. Kelly, “Critical behaviour in
charging of electric vehicles,” New Journal of Physics, vol. 17, no. 9,
p. 095001, 2015.
[15] A. Aveklouris, Y. Nakahira, M. Vlasiou, and B. Zwart, “Electric vehicle
charging: a queueing approach,” SIGMETRICS Performance Evaluation
Review, vol. 45, no. 2, pp. 33–35, 2017.
[16] M. Baran and F. F. Wu, “Optimal sizing of capacitors placed on a radial
distribution system,” IEEE Trans. Power Del., vol. 4, no. 1, pp. 735–743,
1989.
[17] L. Kleinrock, Queueing systems, volume 2: Computer applications.
Wiley New York, 1976, vol. 66.
[18] S. Borst, “User-level performance of channel-aware scheduling algorithms in wireless data networks,” IEEE/ACM Trans. Netw., vol. 13,
no. 3, pp. 636–647, 2005.
[19] I. S. Bayram, G. Michailidis, M. Devetsikiotis, and F. Granelli, “Electric
power allocation in a network of fast charging stations,” IEEE J. Sel.
Areas Commun., vol. 31, no. 7, pp. 1235–1246, 2013.
[20] E. Yudovina and G. Michailidis, “Socially optimal charging strategies
for electric vehicles,” IEEE Trans. Autom. Control, vol. 60, no. 3, pp.
837–842, 2015.
[21] P. You, Y. Sun, J. Pang, S. Low, and M. Chen, “Battery swapping
assignment for electric vehicles: A bipartite matching approach,” SIGMETRICS Performance Evaluation Review, vol. 45, no. 2, pp. 85–87,
2017.
[22] K. Turitsyn, N. Sinitsyn, S. Backhaus, and M. Chertkov, “Robust
broadcast-communication control of electric vehicle charging,” in 2010
1st IEEE Int. Conf. Smart Grid Commun., 2010, pp. 203–207.
[23] W. Su and M.-Y. Chow, “Performance evaluation of an EDA-based largescale plug-in hybrid electric vehicle charging algorithm,” IEEE Trans.
Smart Grid, vol. 3, no. 1, pp. 308–315, 2012.
[24] E. Sortomme, M. M. Hindi, S. J. MacPherson, and S. Venkata, “Coordinated charging of plug-in hybrid electric vehicles to minimize
distribution system losses,” IEEE Trans. Smart Grid, vol. 2, no. 1, pp.
198–205, 2011.
[25] G. Li and X.-P. Zhang, “Modeling of plug-in hybrid electric vehicle
charging demand in probabilistic power flow calculations,” IEEE Trans.
Smart Grid, vol. 3, no. 1, pp. 492–499, 2012.
[26] D. Said, S. Cherkaoui, and L. Khoukhi, “Queuing model for EVs
charging at public supply stations,” Wireless Commun. and Mobile
Computing Conf., vol. 1, no. 5, pp. 65–70, 2013.
[27] P. Kempker, N. v. Dijk, W. Scheinhardt, H. v. d. Berg, and
J. Hurink, “Optimization of charging strategies for electric vehicles in
PowerMatcher-driven smart energy grids,” in Proc. 9th EAI Int. Conf.
Perf. Eval. Meth. and Tools, 2015.
[28] S. Sojoudi and S. H. Low, “Optimal charging of plug-in hybrid electric
vehicles in smart grids,” IEEE Power and Energy Society General
Meeting, pp. 1–6, 2011.
[29] A. Hassan, Y. Dvorkin, D. Deka, and M. Chertkov, “Chance-constrained
ADMM approach for decentralized control of distributed energy resources,” preprint arXiv:1710.09738, 2017.
[30] F. Moret, T. Baroche, E. Sorina, and P. Pinson, “Chanceconstrained admm approach for decentralized control of distributed
energy resources.” [Online]. Available: http://pierrepinson.com/docs/
Moretetal17PSCC.pdf
[31] W. H. Kersting, Distribution system modeling and analysis. CRC press,
2012.
[32] S. H. Low, “Convex relaxation of optimal power flow–part I: Formulations and equivalence,” IEEE Trans. Control Netw. Syst., vol. 1, no. 1,
pp. 15–27, 2014.
[33] ——, “Convex relaxation of optimal power flow–part II: Exactness,”
IEEE Trans. Control Netw. Syst., vol. 1, no. 2, pp. 177–189, 2014.
[34] B. Fralix and B. Zwart, “Time-dependent properties of symmetric

queues,” Queueing Systems, vol. 67, no. 1, pp. 33–45, 2011.
[35] H. C. Gromoll, P. Robert, B. Zwart, and R. Bakker, “The impact
of reneging in processor sharing queues,” SIGMETRICS Performance
Evaluation Review, vol. 34, no. 1, pp. 87–96, 2006.
[36] M. I. Reiman, “The heavy traffic diffusion approximation for sojourn
times in jackson networks,” Applied probability–computer science: the
interface II, pp. 409–421, 1982.
[37] N. Katoh and T. Ibaraki, Resource allocation problems.
Kluwer
Academic Publishers, Dordrecht, The Netherlands, 1998, vol. 2.
[38] T. Bonald and A. Proutiere, “Insensitivity in processor-sharing networks,” Performance Evaluation, vol. 49, no. 1, pp. 193–209, 2002.
[39] J. W. Cohen, “The multiple phase service network with generalized
processor sharing,” Acta informatica, vol. 12, no. 3, pp. 245–284, 1979.

Angelos Aveklouris studied Mathematics at the
National and Kapodistrian University of Athens,
where he graduated in 2012. In 2015, he earned a
masters degree in Mathematical Modeling in Modern
Technologies and Economics from the department
of Mathematics at the National Technical University of Athens. Later that year, he began his PhD
research project in the department of Mathematics
and Computer Science at the Eindhoven University
of Technology. His research interests are mainly in
stochastic networks and probability theory.

Maria Vlasiou (BSc Mathematics 2002, Hons.; PhD
Mathematics 2006) is an Associate Professor at the
Eindhoven University of Technology, a scientific
staff member at CWI Amsterdam, and a research
fellow of EURANDOM. Her research focuses on the
performance of stochastic processing networks with
layered architectures and on perturbation analysis for
heavy-tailed risk models. Dr. Vlasiou serves as the
President of the Dutch association of female professional mathematicians (EWM-NL), as a member
of the scientific advisory board in mathematics of
the Lorenz Center, and as associate editor for IISE Transactions and for
Mathematical Methods of Operations Research.

Bert Zwart (MA Econometrics 1997, PhD Mathematics 2001, Member of IEEE since 2017) is leader
of the CWI Stochastics group (Amsterdam) and Professor at the Eindhoven University of Technology.
He is Stochastic Models area editor for Operations
Research, the flagship journal of his profession,
since 2009. His research expertise is in applied
probability and stochastic networks. His work on
power systems is focusing on the applications of
probabilistic methods to reliability issues and rare
events such as cascading failures and blackouts. He
is co-organizer of a special semester on the Mathematics of energy systems
taking place in Cambridge UK, spring 2019.

