Critical behaviour in charging of electric vehicles

4.6 MB·English

by  Rui Carvalho

Critical behaviour in charging of electric vehicles
Rui Carvalho,1,∗ Lubos Buzna,2,† Richard Gibbens,3,‡ and Frank Kelly4,§
1School of Engineering and Computing Sciences, Durham University,
Lower Mountjoy, South Road, Durham, DH1 3LE, UK
2University of Zilina, Univerzitna 8215/1, 01026 Zilina, Slovakia
5
1 3Computer Laboratory, University of Cambridge, William Gates Building,
0
2 15 JJ Thomson Avenue, Cambridge, CB3 0FD, UK

l
u 4Statistical Laboratory, Centre for Mathematical Sciences,
J

6 University of Cambridge, Wilberforce Road, Cambridge CB3 0WB, UK

]
C The increasing penetration of electric vehicles over the coming decades, taken together
O
with the high cost to upgrade local distribution networks and consumer demand for home
.
h
t charging, suggest thatmanagingcongestion onlowvoltagenetworks willbeacrucialcom-
a
m
ponent ofthe electric vehicle revolution andthe moveawayfrom fossil fuels intransporta-
[

tion. Here, we model the max-flow and proportional fairness protocols for the control of
2
v
7 congestion caused by a fleet of vehicles charging on two real-world distribution networks.
5
9 We show that the system undergoes a continuous phase transition to a congested state as
6
0 a function of the rate of vehicles plugging to the network to charge. We focus on the or-
.
1
0 der parameter and its fluctuations close to the phase transition, and show that the critical
5
1 point depends on the choice of congestion protocol. Finally, we analyse the inequality in
:
v
the charging times as the vehicle arrival rate increases, and show that charging times are
i
X
r considerably moreequitable inproportional fairnessthaninmax-flow.
a
PACSnumbers:88.85.Hj,88.80.Kg,89.75.-k,05.45.-a,
∗[email protected]
†[email protected]
‡[email protected]
§[email protected]
2
I. INTRODUCTION
Electric vehicles may become competitive, in terms of total ownership costs, with internal-
combustion engine vehicles overthe next couple of decades. Studies in the United States and the
UKsuggestthecurrentpowergridhasenoughgenerationcapacitytocharge70%ofcarsandlight
trucks overnight, during periods of low demand [1]. A recent survey suggests, however, vehicle
owners prefer home charging, would consider charging their vehicles during the day (typically
between 6 and 10 pm), and are unwilling to accept a charging time of 8 hours [2]. The time to
fullychargethebatteryofanelectricvehicleathomecurrentlyvariesfrom18hours(Level1,inthe
UnitedStatesat 110V and15 Awithachargepowerof1.4kW)to4hours(Level2, at220V, 30
Awithachargepowerof6.6kW).Alternatively,electricvehiclescouldchargeatpubliccharging
stations at Level 3 in less than 30 minutes [3]. Taken together, consumer behaviour and advances
inbatterytechnologymayleadtoariseinpeakdemandwiththeincreasingpenetrationofelectric
vehicles, overloading distribution networks and requiring local infrastructure reinforcement [4–
7]. To reduce the cost of upgrades to the last mile of cables, network operators may need to
coordinatechargingstrategiesinawaythatisbothtechnicallyandsociallyacceptable. Toachieve
this goal, network designers could implement charging protocols that prioritise the access of a
fleet of electric vehicles to electric power, thus simultaneouslymanaging network congestion and
accountingforthefairness ofuserallocations.
Through a series of papers, the power grid has recently gained increased visibility in the sci-
entificcommunity[8, 9], and physicistshavehelped to increase ourunderstandingofits synchro-
nization [10, 11] and stability [12, 13]. In parallel, recent advances in optimization and phase
transitions [14, 15] suggest that the tools of critical phenomena and optimization can be merged,
openingupnewhorizons. Fromthepointofviewofthedistributionnetworkoperator,theproblem
ofvehiclechargingistomanagecongestionondistributionnetworks,whilerespectingKirchhoff’s
laws and keeping voltage drops bounded. Here, we explore two congestion control mechanisms:
max-flow and proportional fairness. We show that if too many vehicles plug-in to the network,
charging takes too long, more cars arrive than leave fully charged, and the system undergoes a
continuous phase transition to a congested state [16, 17], where the critical point depends on the
choice of congestion control algorithm. By gaining insights into the critical behaviour that nat-
urally emerges with the increase of the number of vehicles, we hope to help network designers
decidewhichalgorithmsto implementinthereal-world.
3
II. THEMODEL
Physicists are familiar with simulated annealing, a global optimization method that can avoid
becomingtrappedinalocaloptimum. Inprinciple,itconvergestotheglobaloptimum,butinprac-
ticethisisnotguaranteed(seee.g.[18–21])becausetherequiredtheoreticalcoolingschedulesare
too slow to use in implementations. In contrast, convex optimization always finds the solution,
if it exists, independently of the starting point. Convex optimization problems can be solved effi-
ciently(typicallyinpolynomialtime),evenforproblemswithhundredsofvariablesandthousands
of constraints, using interior-point methods [22]. The burgeoning field of convex optimization in
electricity networks [23–25] is a good example of an application of the mathematical framework
developedoverthelast 20 years. Indeed, the extensivenumerical simulationswe present here are
only possible due to techniques developed since 2012 [24–26]. The networks that we study are
relatively small. The stochasticity of vehicle arrival times, however, implies solving an optimiza-
tion problem in each time step if the state of the system changes. Hence, to gain insights into the
steadystateofvehiclecharging,efficientalgorithmsareanecessityatthedesignstage. Ofcourse,
real-world implementationsalso depend on efficient algorithms, which would need to run online,
often inlargeurban distributionnetworks.
An optimizationproblem is determined by a function of a set of variables (the objectivefunc-
tion),for whichweseek a minimum,and a setof upperbound constraintsthat restrict thedomain
(or feasible set) of those variables [22]. A point is feasible if it belongs to the feasible set, and is
optimalifitistheminimumoftheobjectivefunctioninthefeasibleset. Anoptimizationproblem
isconvexifboththeobjectivefunctionandtheconstraintsareconvex,inwhichcasetheobjective
function has a global minimum. A convex relaxation of an optimization problem P is a convex
optimization problem P′ with an enlarged feasible set. If the optimum of P′ is feasible for P,
it is also the optimum for P and we say the relaxation is exact. Hence, convex relaxations are
more attractive than approximate methods, such as linearisations, because the feasibility of the
relaxed optimum of P′, which can be verified either analytically or numerically, is a certificate of
theexactnessoftherelaxation.
Consider a tree topology, such that electric power is distributed from a root node to electric
vehiclesthatchargeatthenodes. LetP(t)bethefeasiblesetofpowerallocationsattimet,i.e.the
set of all allocations of power to electric vehicles that do not violate the operational constraints
of the distribution network. Each feasible allocation P(t) ∈ P(t) is a vector P(t) = (P(t) : l =
l
4
1,...,N(t)), where N(t) is the number of vehicles in the network at time t. Vehicle l derives
a utility U (P(t)) from the allocated charging power P(t), and we wish to select the allocation
l l l
that maximises the sum of vehicle utilities [27]. This allocation acts as a network protocol that
distributesnetworkcapacity amongusers, and solvesthefollowingproblem:
N(t)
maximise U (P(t)) (1a)
l l
Xl=1
subjectto P(t) ∈ P(t). (1b)
Hereweexploretwouserutilityfunctions. First,weconsiderthenon-uniquemax-flowallocations
givenby U (P(t)) = P(t), i.e.wemaximisetheinstantaneousaggregatepowersent fromtheroot
l l l
nodetothevehicles,whichisabenchmarkofefficientnetworkthroughput[28]. Suchallocations,
however, can also leave users with zero power, which is considered unfair from the user point of
view. Hence, we next consider the proportionalfairness allocation. Mathematically, the problem
is to find the feasible allocation that maximises the sum of the logarithm of user rates, that is
U (P(t)) = log(P(t)). The proportional fairness allocation is especial, because the users and the
l l l
network operator simultaneouslymaximisetheir utility functions [27]. Furthermore, the problem
isconvex,andsocanbesolvedinpolynomialtime[22],anditcanbenaturallyextendedbyadding
positive weights to each term in the objective function Eq. (1a), to account for diversity in user
demandorformorethanoneuseratsomenodes[27]. ForthecompactandconvexsetP(t),itcan
beshownthat theallocation PPF(t)thatmaximisesEq. (1a), satisfies [27, 29]:
N(t) P(t)−PPF(t)
l l ≤ 0. (2)
PPF(t)
Xl=1 l
This allocation is known as proportionally fair, because the aggregate of proportional changes
with respect to all other feasible allocations is non-negative. In other words, Eq. (2) implies that
toincreasetheinstantaneouspowerallocatedtoavehiclebyapercentageǫ,wehavetodecreasea
setofotherpowerallocations,suchthatthesumofthepercentagedecreasesislargerorequaltoǫ.
In contrast, in max-flow, to increase the instantaneous power allocated to a vehicle by ǫ, we have
to decrease the sum of instantaneous powers allocated to other vehicles at least by ǫ. It turns out
thatfairnessandcongestioncontrolaretwosidesofthesamecoin,becausetheexistingalgorithms
for fair allocations also managenetwork congestion [27, 30–35]. Moreover, in the analysis of the
parallelproblemforcommunicationnetworks,proportionalfairnesshasemergedasacompromise
between efficiency and fairness with an attractive interpretation in terms of shadow prices and a
market clearingequilibrium[27, Section 7.2].
5
(a) (b)
i
P
i
j
V
i
V
j
P
j
P , Q
(j) (j)
FIG. 1. Schematic illustration of (a) a distribution network, (b) the circuit of a network edge. Electric
vehicles choose a charging node with uniform probability, and plug-in to the node until fully charged, as
illustrated by the electric vehicle icons on the network. Network edge e has impedance Z = R +iX .
ij ij ij ij
Thepowerconsumed bythesubtree⋔ (j)rootedatnode j(areashadedinpurple)isS⋔(j) = P⋔(j)+iQ⋔(j),
wherevehiclesconsumerealpoweronly,butnetworkedgeshavebothactive(real)andreactive(imaginary)
powerlosses.
The simplest flow model in electricity networks is the DC powerflow, which is a linearisation
of the AC power flow equations, and thus can be solved with tools of linear programming. It
assumes that voltage amplitude is constant on all nodes, a good approximation at the level of the
high-voltagetransmissionnetwork, butapooroneon localdistributionnetworks. Indeed, voltage
dropsare significantin distributionnetworks,and determinethenetwork capacity,which leads us
to using models of power flow specific to distribution networks [36]. We abstract the distribution
networktoarooteddirectedtree⋔ (r)withnode(oftencalledbus)setV,edge(alsocalledbranch)
set E, and a root noder (feeder) that injects powerinto thetree 1. Edgee ∈ E connects node i to
ij
node j, where i is closer to the root than j, and is characterised by the impedance Z = R +iX ,
ij ij ij
where R is the edge resistance and X the edge reactance. The power loss along edge e is
ij ij ij
given by S (t) = P (t)+iQ (t), where P (t) is the real power loss, and Q (t) the reactive power
ij ij ij ij ij
loss. Electric vehicle l receives active power P(t) until charged, but does not consume reactive
l
power[37]—seeFig.1(a). ThevectorV(t)denotesthevoltageallocatedtothenodes. Thevoltage
drop ∆V down theedge e is the difference between the amplitudeof thevoltagephasors V and
ij ij i
1Wewrite⋔(r)insteadof⋔(V,E,r)tosimplifythenotation.
6
V , assuming node i is closer to the root r than node j [36]. Let ⋔ (j) denote the subtree of the
j
distributionnetwork rooted in node j, with node set V and edge set E . Let P denote the
⋔(j) ⋔(j) ⋔(j)
activepower,and Q thereactivepowerconsumedbythesubtree⋔ (j)—seeFig.1. Kirchhoff’s
⋔(j)
voltagelawapplied tothecircuitin Fig.1(b)yields(seeAppendixA):
V(t)V (t)−V2(t)−P (t)R −Q (t)X = 0. (3)
i j j ⋔(j) ij ⋔(j) ij
Vehicle l has a battery with capacity B that charges with the instantaneous power P(t) from
l
empty(atarrivaltime)tofull(atdeparturetime),andthelevelofbatterychargeisthetimeintegral
of instantaneous power. Vehicles arrive to the network, choose a node to charge randomly with
uniformprobability,chargeuntiltheirbatteryisfull,anlastlyleavethenetwork. Ateachtimestep,
thenetworksolvesthecongestioncontrolproblemtoallocateinstantaneouspowertothevehicles.
The max-flow problem maximises the instantaneous aggregate power sent from the root node to
the electric vehicles, respecting the constraints of distribution networks: the voltage drop along
edges obeys Eq. (3), and node voltages are within ((1−α)V ,(1+α)V ) for α ∈ (0,1),
nominal nominal
with α = 0.1 typically [36]. Thus, to find the max-flow allocation of power to the vehicles, we
solvetheoptimizationproblemforfixed t:
N(t)
maximise U(t) = P(t) (4a)
l
V(t) Xl=1
subjectto (1−α)V ≤ V(t) ≤ (1+α)V , i ∈ V (4b)
nominal i nominal
V(t)V (t)−V (t)V (t)−P (t)R −Q (t)X = 0, e ∈ E. (4c)
i j j j ⋔(j) ij ⋔(j) ij ij
Constraint (4b) sets the limits on the node voltage. Equation (4c) is the physical law coupling
voltage to power, generalized from Eq. (3) for the subtree ⋔ (j), and encodes both Kirchhoff’s
voltagelaw on network edges and Kirchhoff’s current law applied recursively at each node of the
subtree (see Appendix B). We do not need to apply Kirchhoff’s voltage law on network loops,
however, because local distribution networks are approximately trees, and thus are cycle-free.
Constraint (4c) is quadratic, hence not convex in general, which implies that the problem is not
directly solvable by polynomial time methods. To overcome this limitation, we make a change
of variables in problem (4) by defining a weighted adjacency matrix W(t), such that edge e
ij
correspondsto the2×2 principalsubmatrixW(e ,t)defined by [38, 39]:
ij
V(t) V2(t) V(t)V (t) W (t) W (t)
W(eij,t) =  i  Vi(t) Vj(t) =  i i j  =  ii ij , (5)
Vj(t)(cid:18) (cid:19) Vj(t)Vi(t) V2j(t)  Wji(t) Wjj(t)
7
where W (t) = W (t), because V(t),V (t) ∈ R. The matrices W(e ,t) are positive semidefinite,
ij ji i j ij
because their eigenvalues (λ = 0 and λ = V2 + V2) are non-negative, and rank one because
1 2 i j
they are of the form vvT. Hence, constraint (4c) can be replaced by three constraints: the first
substitutesthequadratictermsinthevoltageswithlineartermsintheW(e ,t),andthesecondand
ij
thirdconstraintsguaranteethattheW(e ,t)arepositivesemidefiniteandrank one.
ij
The solutionof problem (4) is on the Pareto frontier 2, since we maximisean increasing func-
tion in the objective. The rank one constraint is nonconvex, but it does not change the Pareto
frontierortheoptimum[25, 40], and weremoveittorelax problem(4) to:
N(t)
maximise U(t) = P(t) (6a)
l
W(t) Xl=1
subjectto ((1−α)V )2 ≤ W (t) ≤ ((1+α)V )2, i ∈ V (6b)
nominal ii nominal
W (t)−W (t)−P (t)R −Q (t)X = 0, e ∈ E (6c)
ij jj ⋔(j) ij ⋔(j) ij ij
W(e ,t) (cid:23) 0, e ∈ E, (6d)
ij ij
where the generalized inequality in constraint (6d) means the W(e ,t) matrices are positive
ij
semidefinite[41].
Theproblemofallocatingpowertovehiclesinaproportionalfairwayhasthesameconstraints
as problem (6), however, the objective function is the sum of the logarithm of the power. It turns
out, however, that it is computationally more efficient to aggregate vehicles at the nodes, and to
maximise the sum of power allocated to the nodes, rather than the vehicles. To show this, we
observe that all vehicles are equivalent, and thus the power P(t) allocated to node i is divided
i
equallyamongthevehicleschargingonthenodeateachtimestep. Hence,ifoneormorevehicles
ischarging on nodei, each gets theinstantaneouspower:
P(t)
P(t) = i , (7)
l
w(t)
i
where w(t) = N(t)∆ (t) is the number of electric vehicles charging on node i at time t, and
i l=1 il
∆ (t) = 1 if elePctric vehicle l is charging on node i at time t and zero otherwise. Hence, the
il
2Wesaythatapowerallocation{P}forl = 1,...,N isbetterthananother{P′}if P ≥ P′ foralll,andforsomel,
l l l l
P > P′. ApowerallocationisParetooptimalorefficientifthereisnobetterpowerallocation. TheParetofrontier
l l
ofasetisthesetofallParetooptimalpoints.
8
proportionalfair allocationis givenby (seeAppendixC):
maximise U(t) = w(t)logP(t) (8a)
i i
W(t) iX∈V+
subjectto ((1−α)V )2 ≤ W (t) ≤ ((1+α)V )2, i ∈ V (8b)
nominal ii nominal
W (t)−W (t)−P (t)R −Q (t)X = 0, e ∈ E (8c)
ij jj ⋔(j) ij ⋔(j) ij ij
W(e ,t) (cid:23) 0, e ∈ E, (8d)
ij ij
where V+ is the subset of nodes with at least one vehicle, and we can recover the instantaneous
power allocated to electric vehicle l, located at node i, from Eq. (7). The complexity of the prob-
lem(8)thusscales withthenumber|V| ofnodes,which istypicallysmallerthanthenumber N(t)
of vehicles for large arrival rates λ. Similarly, we also aggregated vehicles in the implementation
ofproblem(6), butomittheproof.
To study the behaviour of max-flow and proportional fairness as a function of the number of
vehicles arriving at the network to be charged, we implement a discrete simulator that solves
the congestion control problem in discrete time steps, starting with no vehicles charging on the
network. Vehicles arriveat thenetwork in continuoustime(followingaPoisson process withrate
λ)andwithemptybatteries,chooseanodewithuniformprobabilityamongstallnodes(excluding
the root), and charge at that node until their battery is full, at which point in time they leave the
network. Once a vehicle plugs into a node, the congestion control algorithm will allocate it an
instantaneouspower, which is a function of the network topologyand electrical elements, as well
as thelocationofothervehicles.
Ateachtimestep,wefirstcheckwhetherthenumberofchargingvehicleschanged(i.e.vehicles
left the network fully charged, or new vehicles arrived to be charged), and if it has, we solve the
max-flow problem (6) and the proportional fairness problem (8), which allocate a constant power
during the time step to each of the charging vehicles. Next, we update the status of batteries at
the end of the time step. The simulation terminates when the simulation time reaches the time
horizon. Wesimulatedvehiclescharging ontherealisticSCE 47-busandSCE 56-busdistribution
networks [38], which are detailed in Fig. 2. To characterise the system behaviour in detail, we
varied the arrival rate λ from 0 to 1 in steps of 0.05 (0.005 close to the critical points), and for
each λ value we simulated an ensemble of 25 independent realisations of simulation runs, each
simulation running for 15,000 time units (150,000 time units close to the critical point). We ran
9
(a) SCE 47-bus 1 (b) SCE 56-bus
1 35
12 2 13 30 32 34 36 37
17 16 19 18 17 2 3 33 38 39 40 45
47 11 15 3 14 22 23
19 18 21 29 41 42
10 46 34 4 20 24 16 15 56 4 20 23 25 26 43 44 46
44 25 21 24
43 42 9 33 26 5 28 14 13 7 22 27 47 48
56
45 36 41 32 12 11 10 8 9 28 49 51 53 55
37 35 8 7 6 27 31 54
31 30 29 50 52
38 40
39
FIG.2. Topologyofthe(a)SCE47-busand(b)SCE56-busnetworks. Nodeindexesidentifytheedges,and
edge resistance and reactance is taken from [38]. Node 1isthe root node inboth networks. Nodes 13, 17,
19, 23 and 24 of the SCE 47-bus network (in lighter colour) are photovoltaic generators, and we removed
themfromthenetwork.
the simulationsusing CVXOPT [42] on theETHZ Brutus cluster3 due to thehigh computational
requirements. The computational time is comparable for max-flow and proportional fairness and
for the 47-bus and 56-bus networks, but it is grows with λ. For example, to simulate 5,000 time
unitsoftheproportionalfairnessalgorithmforλ = 1.0onthe47-busnetworktakesapproximately
40 hours,but4 minutesforλ = 0.05.
Wesetthebatterycapacity B = 1forallvehicles,andthenominalvoltageV = 1. Scaling
nominal
V by β, for β ∈ (0,∞), implies scaling both the the power delivered to vehicles and the
nominal
batterycapacitybyβ2. Toseethis,observethatproblems(4)and(8)areinvariantuponthescaling
V′ = βV , V′(t) = βV(t) for all nodal voltages, P′(t) = β2P(t), P′(t) = β2P (t) and
nominal nominal i i l l ⋔ ⋔
Q′(t) = β2Q (t), and B′ = β2B. Considering these scaling properties, our simulations can be
⋔ ⋔
extendedtovaluesofV , 1,providedthevehiclecapacity Bisrescaled accordingly,andwe
nominal
usethispropertytorescale theproblemwhenconvenient.
III. NUMERICALRESULTS
Wefindcriticalbehaviourthatresemblesresultsfoundincommunicationnetworks,inthatboth
systemsundergo a continuous phasetransition [43]. In order to characterize this phase transition,
weadopt theorderparameter η(λ)that represents theratio at thesteadystatebetween thenumber
3https://www1.ethz.ch/id/services/list/comp_zentral/cluster/index_EN
10
FIG.3. (a)Orderparameterηasafunctionofthearrivalrateλ,fortheSCE56-bus(filledsymbols)and47-
bus(unfilledsymbols)networks,whereweapplythemax-flow(circularsymbols)andproportional fairness
(squaresymbols)algorithmsforthesimulationhorizonof1.5×104timeunits. Weplotazoomofthecritical
region for the (b) 56-bus network and (c) 47-bus network for the longer horizon of 105 time units. Panel
(c)suggests thecritical arrival rateisdifferent forthemax-flowandproportional fairness algorithms inthe
47-busnetwork. Symbolsshowaveragevaluesoveranensembleof25runsandshadedareasrepresent95%
confidence intervals.
ofuncharged vehiclesand thenumberofvehiclesthat arriveat thenetwork tobecharged [43]:
1h∆N(t)i
η(λ) = lim , (9)
t→∞ λ ∆t
where ∆N(t) = N(t + ∆t) − N(t) and h...i indicates an average over time windows of width ∆t.
We calculate η(λ) in the steady state, that is lim ∆N(t) ∝ ∆t. For arrival rates λ < λ , all
t→∞ c