---
abstract: 'The increasing penetration of electric vehicles over the coming decades, taken together with the high cost to upgrade local distribution networks and consumer demand for home charging, suggest that managing congestion on low voltage networks will be a crucial component of the electric vehicle revolution and the move away from fossil fuels in transportation. Here, we model the max-flow and proportional fairness protocols for the control of congestion caused by a fleet of vehicles charging on two real-world distribution networks. We show that the system undergoes a continuous phase transition to a congested state as a function of the rate of vehicles plugging to the network to charge. We focus on the order parameter and its fluctuations close to the phase transition, and show that the critical point depends on the choice of congestion protocol. Finally, we analyse the inequality in the charging times as the vehicle arrival rate increases, and show that charging times are considerably more equitable in proportional fairness than in max-flow.'
author:
- Rui Carvalho
- Lubos Buzna
- Richard Gibbens
- Frank Kelly
title: Critical behaviour in charging of electric vehicles
---

Introduction
============

Electric vehicles may become competitive, in terms of total ownership costs, with internal-combustion engine vehicles over the next couple of decades. Studies in the United States and the UK suggest the current power grid has enough generation capacity to charge $70\%$ of cars and light trucks overnight, during periods of low demand [@Service09]. 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 [@Deloitte10]. The time to fully charge the battery of an electric vehicle at home currently varies from 18 hours (Level 1, in the United States at 110 V and 15 A with a charge power of 1.4 kW) to 4 hours (Level 2, at 220 V, 30 A with a charge power of 6.6 kW). Alternatively, electric vehicles could charge at public charging stations at Level 3 in less than 30 minutes [@Dickerman10]. Taken together, consumer behaviour and advances in battery technology may lead to a rise in peak demand with the increasing penetration of electric vehicles, overloading distribution networks and requiring local infrastructure reinforcement [@Clement10; @Green11; @Tran12; @Keshav12]. To reduce the cost of upgrades to the last mile of cables, network operators may need to coordinate charging strategies in a way that is both technically and socially acceptable. To achieve this goal, network designers could implement charging protocols that prioritise the access of a fleet of electric vehicles to electric power, thus simultaneously managing network congestion and accounting for the fairness of user allocations.

Through a series of papers, the power grid has recently gained increased visibility in the scientific community [@DSouza13; @Baptista14], and physicists have helped to increase our understanding of its synchronization [@Motter13; @Rohden12] and stability [@Sole08; @Kurths14]. In parallel, recent advances in optimization and phase transitions [@Scala14; @Seoane14] suggest that the tools of critical phenomena and optimization can be merged, opening up new horizons. From the point of view of the distribution network operator, the problem of vehicle charging is to manage congestion on distribution networks, while respecting Kirchhoff’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 [@Guimera02; @Lai05], where the critical point depends on the choice of congestion control algorithm. By gaining insights into the critical behaviour that naturally emerges with the increase of the number of vehicles, we hope to help network designers decide which algorithms to implement in the real-world.

The Model {#sec:model}
=========

Physicists are familiar with simulated annealing, a global optimization method that can avoid becoming trapped in a local optimum. In principle, it converges to the global optimum, but in practice this is not guaranteed (see  [@Hajek88; @Boyd01; @Donetti05; @Arenas08]) because the required theoretical cooling schedules are 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 efficiently (typically in polynomial time), even for problems with hundreds of variables and thousands of constraints, using interior-point methods [@Boyd04]. The burgeoning field of convex optimization in electricity networks [@TaylorBook15; @Low14b; @Low14c] is a good example of an application of the mathematical framework developed over the last 20 years. Indeed, the extensive numerical simulations we present here are only possible due to techniques developed since 2012 [@Lavaei12; @Low14b; @Low14c]. The networks that we study are relatively small. The stochasticity of vehicle arrival times, however, implies solving an optimization problem in each time step if the state of the system changes. Hence, to gain insights into the steady state of vehicle charging, efficient algorithms are a necessity at the design stage. Of course, real-world implementations also depend on efficient algorithms, which would need to run online, often in large urban distribution networks.

An optimization problem is determined by a function of a set of variables (the objective function), for which we seek a minimum, and a set of upper bound constraints that restrict the domain (or *feasible set*) of those variables [@Boyd04]. A point is feasible if it belongs to the feasible set, and is optimal if it is the minimum of the objective function in the feasible set. An optimization problem is convex if both the objective function and the constraints are convex, in which case the objective function has a global minimum. A convex relaxation of an optimization problem $\mathsf{P}$ is a convex optimization problem $\mathsf{P}^\prime$ with an enlarged feasible set. If the optimum of $\mathsf{P}^\prime$ is feasible for $\mathsf{P}$, it is also the optimum for $\mathsf{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 $\mathsf{P}^\prime$, which can be verified either analytically or numerically, is a certificate of the exactness of the relaxation.

Consider a tree topology, such that electric power is distributed from a root node to electric vehicles that charge at the nodes. Let $\mathcal{P}(t)$ be the feasible set of power allocations at time $t$,  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) \in \mathcal{P}(t)$ is a vector $P(t) = (P_l(t): l = 1, \dots, N(t))$, where $N(t)$ is the number of vehicles in the network at time $t$. Vehicle $l$ derives a utility $U_l(P_l(t))$ from the allocated charging power $P_l(t)$, and we wish to select the allocation that maximises the sum of vehicle utilities [@Kelly14]. This allocation acts as a network protocol that distributes network capacity among users, and solves the following problem:

[align]{} \[eq:max\_flow\_final\_a\] & \_[l=1]{}\^[N(t)]{} U\_l(P\_l(t)) &&\
\[eq:max\_flow\_final\_b\] & P(t) (t).

Here we explore two user utility functions. First, we consider the non-unique *max-flow* allocations given by $U_l\left( P_l(t) \right) = P_l(t)$,  we maximise the instantaneous aggregate power sent from the root node to the vehicles, which is a benchmark of efficient network throughput [@Bertsimas10]. Such allocations, however, can also leave users with zero power, which is considered unfair from the user point of view. Hence, we next consider the *proportional fairness* allocation. Mathematically, the problem is to find the feasible allocation that maximises the sum of the logarithm of user rates, that is $U_l\left( P_l(t) \right)= log(P_l(t))$. The proportional fairness allocation is especial, because the users and the network operator simultaneously maximise their utility functions [@Kelly14]. Furthermore, the problem is convex, and so can be solved in polynomial time [@Boyd04], and it can be naturally extended by adding positive weights to each term in the objective function Eq. (\[eq:max\_flow\_final\_a\]), to account for diversity in user demand or for more than one user at some nodes [@Kelly14]. For the compact and convex set $\mathcal{P}(t)$, it can be shown that the allocation $P^{PF}(t)$ that maximises Eq. (\[eq:max\_flow\_final\_a\]), satisfies [@Kelly14; @Luss_2012]: \[eq:proportional\_fairness\_frac\] \_[l=1]{}\^[N(t)]{} 0. 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. (\[eq:proportional\_fairness\_frac\]) implies that to increase the instantaneous power allocated to a vehicle by a percentage $\epsilon$, we have to decrease a set of other power allocations, such that the sum of the percentage decreases is larger or equal to $\epsilon$. In contrast, in max-flow, to increase the instantaneous power allocated to a vehicle by $\epsilon$, we have to decrease the sum of instantaneous powers allocated to other vehicles at least by $\epsilon$. It turns out that fairness and congestion control are two sides of the same coin, because the existing algorithms for fair allocations also manage network congestion [@Bertsekas92; @Kelly98; @Tan99; @Srikant04; @Carvalho12; @Carvalho14; @Kelly14]. Moreover, in the analysis of the parallel problem for communication networks, proportional fairness has emerged as a compromise between efficiency and fairness with an attractive interpretation in terms of shadow prices and a market clearing equilibrium [@Kelly14 Section 7.2].

![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_{ij}$ has impedance $Z_{ij}=R_{ij}+iX_{ij}$. The power consumed by the subtree $\pitchfork(j)$ rooted at node $j$ (area shaded in purple) is $S_{\pitchfork(j)} = P_{\pitchfork(j)} + iQ_{\pitchfork(j)}$, where vehicles consume real power only, but network edges have both active (real) and reactive (imaginary) power losses.[]{data-label="Fig:tree_network"}](fig1.eps){width="\textwidth"}

The simplest flow model in electricity networks is the DC power flow, 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-voltage transmission network, but a poor one on local distribution networks. Indeed, voltage drops are significant in distribution networks, and determine the network capacity, which leads us to using models of power flow specific to distribution networks [@Kersting01]. We abstract the distribution network to a rooted directed tree $\pitchfork(r)$ with node (often called *bus*) set $\mathcal{V}$, edge (also called *branch*) set $\mathcal{E}$, and a root node $r$ (*feeder*) that injects power into the tree [^1]. Edge $e_{ij}\in \mathcal{E}$ connects node $i$ to node $j$, where $i$ is closer to the root than $j$, and is characterised by the impedance $Z_{ij}=R_{ij}+iX_{ij}$, where $R_{ij}$ is the edge resistance and $X_{ij}$ the edge reactance. The power loss along edge $e_{ij}$ is given by $S_{ij}(t) = P_{ij}(t) + i Q_{ij}(t)$, where $P_{ij}(t)$ is the real power loss, and $Q_{ij}(t)$ the reactive power loss. Electric vehicle $l$ receives active power $P_l(t)$ until charged, but does not consume reactive power [@Low11] —see Fig. \[Fig:tree\_network\](a). The vector $V(t)$ denotes the voltage allocated to the nodes. The voltage drop $\Delta V_{ij}$ down the edge $e_{ij}$ is the difference between the amplitude of the voltage phasors $V_i$ and $V_j$, assuming node $i$ is closer to the root $r$ than node $j$ [@Kersting01]. Let $\pitchfork(j)$ denote the subtree of the distribution network rooted in node $j$, with node set $\mathcal{V}_{\pitchfork(j)}$ and edge set $\mathcal{E}_{\pitchfork(j)}$. Let $P_{\pitchfork (j)}$ denote the active power, and $Q_{\pitchfork(j)}$ the reactive power consumed by the subtree $\pitchfork(j)$—see Fig. \[Fig:tree\_network\]. Kirchhoff’s voltage law applied to the circuit in Fig. \[Fig:tree\_network\](b) yields (see Appendix \[appendix:voltage\_drop\]): \[eq:voltage\_drop\] V\_i(t)V\_j(t) - V\_j\^2(t) -P\_[(j)]{}(t) R\_[ij]{} - Q\_[(j)]{}(t) X\_[ij]{} = 0.

Vehicle $l$ has a battery with capacity $B$ that charges with the instantaneous power $P_l(t)$ from empty (at arrival time) to full (at departure time), and the level of battery charge is the time integral of instantaneous power. Vehicles arrive to the network, choose a node to charge randomly with uniform probability, charge until their battery is full, an lastly leave the network. At each time step, the network solves the congestion control problem to allocate instantaneous power to the vehicles. 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. (\[eq:voltage\_drop\]), and node voltages are within $\left ( (1-\alpha)V_{nominal}, (1+\alpha)V_{nominal} \right )$ for $\alpha\in(0,1)$, with $\alpha=0.1$ typically [@Kersting01]. Thus, to find the max-flow allocation of power to the vehicles, we solve the optimization problem for fixed $t$: \[eq:problem\_max-flow\]

[align]{} & U(t) = \_[l=1]{}\^[N(t)]{} P\_l(t) &&\[eq:problem\_max-flow\_a\]\
& (1-) V\_[nominal]{} V\_i(t) (1+) V\_[nominal]{}, && i \[eq:problem\_max-flow\_b\]\
& V\_i(t) V\_j(t) - V\_j(t)V\_j(t) - P\_[(j)]{}(t) R\_[ij]{} - Q\_[(j)]{}(t) X\_[ij]{} = 0 , && e\_[ij]{} . \[eq:problem\_max-flow\_c\]

Constraint (\[eq:problem\_max-flow\_b\]) sets the limits on the node voltage. Equation (\[eq:problem\_max-flow\_c\]) is the physical law coupling voltage to power, generalized from Eq. (\[eq:voltage\_drop\]) for the subtree $\pitchfork (j)$, and encodes both Kirchhoff’s voltage law on network edges and Kirchhoff’s current law applied recursively at each node of the subtree (see Appendix \[appendix:subtree\]). 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 (\[eq:problem\_max-flow\_c\]) 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 (\[eq:problem\_max-flow\]) by defining a weighted adjacency matrix $W(t)$, such that edge $e_{ij}$ corresponds to the $2\times 2$ principal submatrix $W(e_{ij},t)$ defined by [@Low14a; @BoseThesis14]: \[eq:def\_W\] W(e\_[ij]{},t) = V\_i(t)\
V\_j(t) V\_i(t) & V\_j(t) = V\_i\^2(t) & V\_i(t)V\_j(t)\
V\_j(t)V\_i(t) & V\_j\^2(t)\
= W\_[ii]{}(t) & W\_[ij]{}(t)\
W\_[ji]{}(t) & W\_[jj]{}(t)\
, where $W_{ij}(t)=W_{ji}(t)$, because $V_i(t),V_j(t) \in \mathbb{R}$. The matrices $W(e_{ij},t)$ are positive semidefinite, because their eigenvalues ($\lambda_1=0$ and $\lambda_2=V_i^2+V_j^2$) are non-negative, and rank one because they are of the form $v v^T$. Hence, constraint (\[eq:problem\_max-flow\_c\]) can be replaced by three constraints: the first substitutes the quadratic terms in the voltages with linear terms in the $W(e_{ij},t)$, and the second and third constraints guarantee that the $W(e_{ij},t)$ are positive semidefinite and rank one.

The solution of problem (\[eq:problem\_max-flow\]) is on the *Pareto frontier* [^2], since we maximise an increasing function in the objective. The rank one constraint is nonconvex, but it does not change the Pareto frontier or the optimum [@Lavaei14; @Low14c], and we remove it to relax problem (\[eq:problem\_max-flow\]) to: \[eq:max\_flow\_final\]

[align]{} & U(t) = \_[l=1]{}\^[N(t)]{} P\_l(t) &&\
& ( (1-) V\_ )\^2 W\_[ii]{}(t) ( (1+) V\_ )\^2, && i\
& W\_[ij]{}(t) - W\_[jj]{}(t) - P\_[(j)]{}(t) R\_[ij]{} - Q\_[(j)]{}(t) X\_[ij]{} = 0, && e\_[ij]{}\
& W(e\_[ij]{},t) 0, && e\_[ij]{} , && \[eq:psd\_W\]

where the generalized inequality in constraint (\[eq:psd\_W\]) means the $W(e_{ij},t)$ matrices are positive semidefinite [@Strang06].

The problem of allocating power to vehicles in a proportional fair way has the same constraints as problem (\[eq:max\_flow\_final\]), 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 $\mathrm{P_i}(t)$ allocated to node $i$ is divided equally among the vehicles charging on the node at each time step. Hence, if one or more vehicles is charging on node $i$, each gets the instantaneous power: \[eq:P\_l\] P\_l(t)=, where $w_i(t)= \sum_{l=1}^{N(t)} \Delta_{il}(t)$ is the number of electric vehicles charging on node $i$ at time $t$, and $\Delta_{il}(t)=1$ if electric vehicle $l$ is charging on node $i$ at time $t$ and zero otherwise. Hence, the proportional fair allocation is given by (see Appendix \[appendix:aggregation\_vehicles\]): \[eq:problem\_proportional\_fairness\]

[align]{} & U(t) = \_[i \^+]{} w\_i(t) (t) && \[eq:problem\_proportional\_fairness\_a\]\
& ( (1-) V\_ )\^2 W\_[ii]{}(t) ( (1+) V\_ )\^2, && i \[eq:problem\_proportional\_fairness\_b\]\
& W\_[ij]{}(t) - W\_[jj]{}(t) - P\_[(j)]{}(t) R\_[ij]{} - Q\_[(j)]{}(t) X\_[ij]{} = 0 , && e\_[ij]{} \[eq:problem\_proportional\_fairness\_c\]\
& W(e\_[ij]{},t) 0, && e\_[ij]{} , \[eq:problem\_proportional\_fairness\_d\]

where $\mathcal{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. (\[eq:P\_l\]). The complexity of the problem (\[eq:problem\_proportional\_fairness\]) thus scales with the number $\lvert \mathcal{V} \rvert$ of nodes, which is typically smaller than the number $N(t)$ of vehicles for large arrival rates $\lambda$. Similarly, we also aggregated vehicles in the implementation of problem (\[eq:max\_flow\_final\]), but omit the proof.

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 arrive at the network in continuous time (following a Poisson process with rate $\lambda$) and with empty batteries, choose a node with uniform probability amongst all nodes (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 instantaneous power, which is a function of the network topology and electrical elements, as well as the location of other vehicles.

At each time step, we first check whether the number of charging vehicles changed ( vehicles left the network fully charged, or new vehicles arrived to be charged), and if it has, we solve the max-flow problem (\[eq:max\_flow\_final\]) and the proportional fairness problem (\[eq:problem\_proportional\_fairness\]), 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. We simulated vehicles charging on the realistic SCE 47-bus and SCE 56-bus distribution networks [@Low14a], which are detailed in Fig. \[Fig:networks\]. To characterise the system behaviour in detail, we varied the arrival rate $\lambda$ from 0 to 1 in steps of 0.05 (0.005 close to the critical points), and for each $\lambda$ 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 the simulations using CVXOPT [@cvxopt] on the ETHZ Brutus cluster [^3] due to the high 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 $\lambda$. For example, to simulate 5,000 time units of the proportional fairness algorithm for $\lambda = 1.0$ on the 47-bus network takes approximately 40 hours, but 4 minutes for $\lambda=0.05$ .

We set the battery capacity $B=1$ for all vehicles, and the nominal voltage $V_{nominal}=1$. Scaling $V_{nominal}$ by $\beta$, for $\beta\in (0,\infty)$, implies scaling both the the power delivered to vehicles and the battery capacity by $\beta^2$. To see this, observe that problems (\[eq:problem\_max-flow\]) and (\[eq:problem\_proportional\_fairness\]) are invariant upon the scaling $V_{nominal}^\prime=\beta V_{nominal}$, $V_i^\prime(t)=\beta V_i(t)$ for all nodal voltages, $P_l^\prime(t)=\beta^2P_l(t)$, $P_{\pitchfork}^\prime(t)=\beta^2P_{\pitchfork}(t)$ and $Q_{\pitchfork}^\prime(t)=\beta^2Q_{\pitchfork}(t)$, and $B^\prime=\beta^2 B$. Considering these scaling properties, our simulations can be extended to values of $V_{nominal}\neq 1$, provided the vehicle capacity $B$ is rescaled accordingly, and we use this property to rescale the problem when convenient.

Numerical results
=================

\

We find critical behaviour that resembles results found in communication networks, in that both systems undergo a continuous phase transition [@Guimera01]. In order to characterize this phase transition, we adopt the order parameter $\eta(\lambda)$ that represents the ratio at the steady state between the number of uncharged vehicles and the number of vehicles that arrive at the network to be charged [@Guimera01]: () = \_[t]{} , \[eq:eta\] where $\Delta N(t)=N(t+\Delta t) - N(t)$ and $\left\langle\ldots\right\rangle$ indicates an average over time windows of width $\Delta t$. We calculate $\eta(\lambda)$ in the steady state, that is $\lim_{t\rightarrow\infty} \Delta N(t)\propto  \Delta t$. For arrival rates $\lambda < \lambda_c$, all vehicles that plug-in to the network with empty batteries within a large enough time window leave fully charged within that period ([*free flow*]{} phase), but for $\lambda > \lambda_c$ some vehicles have to wait for increasingly long times to fully charge ([*congested*]{} phase). The order parameter characterises the phase transition: $\eta(\lambda)=0$ in the free-flow regime, and $\eta(\lambda)>0$ in the congested phase, a higher order parameter meaning that queues of charging vehicles build up more rapidly.

Figure \[fig:phase\_diagram\] is a plot of the order parameter for the 47 and 56-bus networks and the two congestion control methods, as a function of $\lambda$. Simulation results shown in Figs. \[fig:phase\_diagram\](a) suggest that $\lambda_c$ depends on several factors (the network topology, the complex impedance on the edges, battery capacity, $V_{nominal}$, as well as the position of vehicles on the network). At this resolution of the control parameter, it is unclear, however, whether the critical point is the same for max-flow and proportional fairness in both networks. To clarify this, we studied the order parameter with higher resolution close to the critical points—see Figs. \[fig:phase\_diagram\](b) and (c). The critical point is numerically indistinguishable for max-flow and proportional fairness in the 56-bus network. In the 47-bus network, however, we find that $\lambda_c$ is larger for proportional fairness than for max-flow.

The number $N(t)$ of charging vehicles at time $t$ fluctuates widely close to the critical point, and thus it is difficult to determine $\lambda_c$ from Fig. \[fig:phase\_diagram\]. To overcome this limitation, we adopt the susceptibility-like function [@Guimera01]: ()=\_[t ]{}t \_(t), where $\Delta t$ is the length of a time window, and $\sigma_\eta (\Delta t)$ is the standard deviation of the order parameter $\eta$. To compute $\chi(\lambda)$, we first consider a long time series and split it into windows with length $\Delta t$. We next determine the value of the order parameter in each window, and finally calculate the standard deviation of these values. The susceptibility displays a singular point at $\lambda_c$ (see Fig. \[fig:susceptibility\]) , allowing us to study the dependencies of the critical arrival rates on the congestion control algorithm, as well as network topology and size.

Similarly to our analysis of $\eta(\lambda)$, the values of $\lambda_c$ are indistinguishable in the 56-bus network. In contrast, however, in the 47-bus network the singular point of $\chi(\lambda)$ is smaller for max-flow than for proportional fairness. This suggests that proportional fairness charges a slightly larger number of vehicles than max-flow, and is thus marginally more efficient, on a neighbourhood of its critical point. To support this conclusion, we show in Figs. \[fig:susceptibility\](c) and (d) four representative instances of the time series of the number of vehicles charging on the 47-bus network at $\lambda=0.39$. The number $N(t)$ of vehicles grows linearly with time in max-flow in all four cases, suggesting that the critical point is below $\lambda=0.39$ for this algorithm. In contrast, $N(t)$ oscillates in proportional fairness, suggesting that the critical point is above $\lambda=0.39$, in agreement with the analysis of $\chi(\lambda)$.

The two congestion control algorithms lead to different allocations of instantaneous power, with vehicles charging in different order and over different time intervals. If there are vehicles on a path $p$ between the root and a leaf node, the voltage drops with increasing distance from the root, the lower limit voltage constraint (\[eq:problem\_max-flow\_b\]) is fulfilled at equality for one node on $p$, and nodes further away than that will not receive any power. The objective function of proportional fairness guarantees that each vehicle gets a positive power allocation, thus the lower limit voltage constraint is satisfied at equality on the occupied node that is the most distant from the root on $p$. In max-flow, however, to maximise the aggregate power allocated to vehicles that can take all instantaneous power they are allocated (elastic demand), on a network with bounded voltage drops ( capacity), implies also minimizing the power losses, and this is achieved by allocating all power on $p$ to the closest occupied node from the root on that path. For max-flow, this implies vehicles on the path $p$ further away from the root than the closest occupied node will only receive power after all vehicles on this node have left the network fully charged. In other words, under max-flow, users experience a charging time that depends strongly on their location on the network: vehicles close to the root charge faster, and vehicles on the tree leaves may take a very long time to charge. In contrast, under proportional fairness, the charging times are more homogeneous, because vehicles receive instantaneous powers that are also more uniform.

![Gini coefficient $G$ of the charging time as a function of the electric vehicle arrival rate $\lambda$, for the SCE 47-bus (unfilled symbols) and 56-bus (filled symbols) networks, where we apply the max-flow (circular symbols) and proportional fairness (square symbols) algorithms. We run the simulation for 15,000 times units, and compute the Gini coefficient from the charging time of vehicles that have charged fully during the simulation. To reduce the effect of a transient regime, we consider only vehicles that are fully charged after iteration $1000$. Vertical lines show the value of the critical points for max-flow (MF) and proportional fairness (PF) identified from the susceptibility $\chi(\lambda)$ for both networks. Symbols show average values over an ensemble of 25 runs and shaded areas represent $95\%$ confidence intervals.[]{data-label="Fig:Gini"}](fig5.eps){width="\textwidth"}

To characterise inequalities in the user experience, we analyse the Gini coefficient of charging time. Originally devised as a measure of inequality in income distributions, the Gini coefficient is defined as [@Ullah98]: G= =\_[0]{}\^\_[0]{}\^u-vf(u)f(v) du dv, \[eq:gini\_coefficient\_continuous\] where $u$ and $v$ are independent identically distributed random variables with probability density $f$ and mean $\mu$. In other words, the Gini coefficient is one half of the mean difference in units of the mean. The difference between the two variables receives a small weight in the tail of the distribution, where $f(u)f(v)$ is small, but a relatively large weight near the mode. Hence, $G$ is more sensitive to changes near the mode than to changes in the tails. For a random sample ($x_i$, $i=1,2, \ldots, n$), the empirical Gini coefficient, $\widehat{G}$, may be estimated by a sample mean =. \[eq:gini\_coefficient\_discrete\] The Gini coefficient is used as a measure of inequality, because a sample where the only non-zero value is $x$ has $\mu=x/n$ and hence $\widehat{G}=(n-1)/n \to 1$ as $n\to\infty$, whereas $\widehat{G}=0$ when all data points have the same value.

We observe in Fig. \[Fig:Gini\] that the Gini coefficient of the charging time is larger in max-flow than in proportional fairness, for each of the networks. Moreover, the Gini coefficient increases faster in max-flow than in proportional fairness in the non-congested regime, showing that, when the system is stable, vehicles will experience a faster increase in the inequality of charging times in max-flow than in proportional fairness, with the increase of the vehicle arrival rate $\lambda$. For comparison with well-known measures of income inequality, Sweden has a Gini of 0.26, the United States has a Gini of 0.41 and the Seychelles has the highest Gini of 0.66 [@WorldBank15]. The proportional fairness algorithm reaches a maximum Gini of 0.45, which is comparable with the level of inequality in the US society, and thus may be judged sociable acceptable. The max-flow algorithm, however, reaches a Gini of 0.91, which measures a level of inequality considerably higher than present in any contemporary society.

Discussion
==========

In conclusion, we modelled the max-flow and proportional fairness protocols for the control of congestion caused by a fleet of vehicles charging on distribution networks. We analysed the second order phase transition that occurs with the increase of the number of electric vehicles that arrive at the network with empty batteries to be charged, and found that the critical arrival rate $\lambda_c$ depends on the congestion control method. Indeed, we showed numerically on the 47-bus bus network that the onset of congestion takes place for larger values of $\lambda$ in proportional fairness than in max-flow. This result is surprising, because one would expect that, for a chosen arrival rate $\lambda$, the maximisation of the aggregate instantaneous power would also lead to a maximisation of the energy (the time integral of power), and hence to a maximisation of the number of charged vehicles. This discovery, illustrates how the greediness of max-flow can be sub-optimal in relation to proportional fairness, which is an example of a fair allocation of instantaneous power.

We analysed the inequality in the charging times as the vehicle arrival rate increases, and showed that charging times are considerably more equitable in proportional fairness than in max-flow. Indeed, vehicles close to the root get all the power allocation in max-flow, leaving other vehicles excluded from the network and unable to charge. Hence, proportional fairness is preferable to max-flow, not only because it does not exclude users from the network, but also because the charging times are more equitable, and it can serve a higher number of vehicles. In conclusion, proportional fairness is a promising candidate protocol to manage congestion in the charging of electric vehicles.

We thank Janusz Bialek, Chris Dent and Emmanouil Loukarakis for help with modelling distribution networks. We also thank Dirk Helbing for granting access to the ETHZ Brutus high-performance cluster. This work was supported by the Engineering and Physical Sciences Research Council under grant number EP/I016023/1, by APVV (project APVV-0760-11) and by VEGA (project 1/0339/13).

Voltage drop on one edge {#appendix:voltage_drop}
========================

The angle $\theta$ between $V_i$ and $V_j$ is small in distribution networks [@Kersting01] (see Fig. \[Fig:phasor\]), and hence the phases of $V_i$ and $V_j$ are approximately the same, and can be chosen so the phasors have zero imaginary components. Since the phasors are real, we can derive the voltage drop from Kirchhoff’s voltage law applied to the circuit in Fig. \[Fig:tree\_network\](b), $$\begin{aligned}
\label{eq:voltage_dropone}
 \Delta V_{ij} &= \rvert V_i\rvert- \rvert V_j\rvert \simeq V_i - V_j \nonumber=\\
& =\Re\left(I_{ij}\text{ }Z_{ij}\right) = \Re\left(\frac{I_{ij}\text{ }Z_{ij}\text{ }V_j^{*}}{V_j^{*}}\right)  =\Re\left(\frac{S_{\pitchfork(j)}^*\text{ }Z_{ij}} {V_j^{*}} \right)= \nonumber\\
& = \Re\left(\frac{(P_{\pitchfork(j)} - iQ_{\pitchfork(j)})(R_{ij} + iX_{ij})}{V_j^{*}}\right)  \simeq \frac{P_{\pitchfork(j)} R_{ij} + Q_{\pitchfork(j)} X_{ij} }{V_j},\end{aligned}$$ where the superscript asterisk denotes the complex conjugate transpose.

![The difference $I_{ij} Z_{ij}$ between the $V_i$ and $V_j$ phasors, decomposed along the $V_j$ vector and its orthogonal direction. The phase angle $\theta$ difference between $V_i$ and $V_j$ is small, and hence the voltage drop can be approximated by $\Delta V_{ij}\simeq \Re\left(I_{ij}\text{ }Z_{ij}\right)$.[]{data-label="Fig:phasor"}](fig6.eps){width="\textwidth"}

Active and reactive loads on a subtree {#appendix:subtree}
======================================

From Kirchhoff’s current law, the active and reactive power consumed by the loads in the subtree rooted in node $k$ can be computed as: $$\begin{aligned}
\label{P_tree}
& P_{\pitchfork (k)} = \sum_{i \in \mathcal{V}_{\pitchfork(k)}} \sum_{l=1}^{N(t)} \Delta_{il}(t) P_l(t) + \sum_{i \in \mathcal{V}_{\pitchfork(k)}} \sum_{j: e_{ij}\in \mathcal{E}_{\pitchfork(k)}} P_{ij}(t),\end{aligned}$$ and $$\begin{aligned}
\label{Q_tree}
& Q_{\pitchfork (k)} =  \sum_{i \in \mathcal{V}_{\pitchfork(k)}} \sum_{j: e_{ij}\in \mathcal{E}_{\pitchfork(k)}} Q_{ij}(t),\end{aligned}$$ where $P_{ij}(t)$ is the active and $Q_{ij}(t)$ the reactive power dissipated on a cable connecting nodes $i$ and $j$. The complex power is given by: $$\begin{aligned}
& S_{ij}(t) = P_{ij}(t) + iQ_{ij}(t) = (V_i(t) - V_j(t)) I^{*} =  \nonumber \\ 
& = (V_i(t) - V_j(t)) \left( \frac{V_i(t) - V_j(t)}{R_{ij} + iX_{ij}} \right)^{*} = \nonumber \\
& = \left (V_i(t) - V_j(t)\right ) \left (V_i(t) - V_j(t)\right )^{*}\left( \frac{R_{ij} - iX_{ij}}{R_{ij}^{2} + X_{ij}^{2}} \right)^{*} = \nonumber \\
& = \left (W_{ii}(t)-W_{ij}(t)-W_{ji}(t)+W_{jj}(t)\right ) \frac{R_{ij} + iX_{ij}}{R_{ij}^{2} + X_{ij}^{2}}.\end{aligned}$$ Since, the voltages are real, $W_{ij}(t)=W_{ji}(t)$, and thus P\_[ij]{}(t) = (W\_[ii]{}(t) -2W\_[ij]{}(t) + W\_[jj]{}(t)), and Q\_[ij]{}(t) = (W\_[ii]{}(t) -2W\_[ij]{}(t) + W\_[jj]{}(t)).

Aggregation of vehicles at the nodes {#appendix:aggregation_vehicles}
====================================

In proportional fairness, we maximise the sum of the logarithm of the instantaneous power allocated to electric vehicles: \[eq:utility\_1\] \_[l=1]{}\^[N(t)]{} P\_l(t) = \_[i ]{} \_[l=1]{}\^[N(t)]{} \_[il]{}(t) , where $P_l(t)$ is the instantaneous power allocated to electric vehicle $l$, and $\mathrm{P_i}$ the instantaneous power allocated to node $i$. To maximise Eq. (\[eq:utility\_1\]), we solve a problem with gradient and Hessian matrices that grow in size with the number of electric vehicles on the network. A more efficient way to approach the problem is to aggregate cars for each node $i$, then solve the optimization problem for the nodes (as if they were ‘super-cars’), and finally distribute the power allocated to each node among the cars on the node. To do this, we remove constant terms in the objective function Eq. (\[eq:utility\_1\]), yielding: \[eq:utility\_2\] U(t) = \_[i \^+]{} \_[l=1]{}\^[N(t)]{} \_[il]{}(t) (t) = \_[i \^+]{} w\_i(t) (t).

[10]{} url \#1[[\#1]{}]{}urlprefix\[2\]\[\][[\#2](#2)]{} Service R F 2009 [*Science*]{} [**324**]{} 1257–1259

Deloitte 2010 Gaining traction: A customer view of electric vehicle mass adoption in the US automotive market. US survey of vehicle owners Tech. Rep. Deloitte Development LLC

Dickerman L and Harrison J 2010 [*IEEE Power & Energy Magazine*]{} [**8**]{} 55–61

Clement-Nyns K, Haesen E and Driesen J 2010 [*IEEE Trans. Power Syst.*]{} [**25**]{} 371–380

Green R C, Wang L F and Alam M 2011 [*Renew. Sust. Energ. Rev.*]{} [**15**]{} 544–553

Tran M, Banister D, Bishop J D K and McCulloch M D 2012 [*Nat. Clim. Chang.*]{} [**2**]{} 328–333

Ardakanian O, Rosenberg C and Keshav S 2012 [*ACM SIGMETRICS Performance Evaluation Review*]{} [**40**]{} 38–42

Brummitt C D, Hines P D H, Dobson I, Moore C and D’Souza R M 2013 [*Proc. Natl. Acad. Sci. U. S. A.*]{} [**110**]{} 12159–12159

Nardelli P H J, Rubido N, Wang C W, Baptista M S, Pomalaza-Raez C, Cardieri P and Latva-aho M 2014 [*Eur. Phys. J.-Spec. Top.*]{} [**223**]{} 2423–2437

Motter A E, Myers S A, Anghel M and Nishikawa T 2013 [*Nat. Phys.*]{} [**9**]{} 191–197

Rohden M, Sorge A, Timme M and Witthaut D 2012 [*Phys. Rev. Lett.*]{} [ **109**]{} 5

Sol[é]{} R V, Rosas-Casals M, Corominas-Murtra B and Valverde S 2008 [*Phys. Rev. E*]{} [**77**]{} 026102–7

Menck P J, Heitzig J, Kurths J and Schellnhuber H J 2014 [*Nat. Commun.*]{} [**5**]{} 3969

Pahwa S, Scoglio C and Scala A 2014 [*Sci Rep*]{} [**4**]{} 9

Seoane L F and Sol[é]{} R V 2014 A multiobjective optimization approach to statistical mechanics http://arxiv.org/abs/1310.6372

Guimera R, Diaz-Guilera A, Vega-Redondo F, Cabrales A and Arenas A 2002 [ *Phys. Rev. Lett.*]{} [**89**]{} 248701

Zhao L, Lai Y C, Park K and Ye N 2005 [*Phys. Rev. E*]{} [**71**]{} 026125

Hajek B 1988 [*Math. Oper. Res.*]{} [**13**]{} 311–329

Hershenson M D, Boyd S P and Lee T H 2001 [*IEEE Trans. Comput-Aided Des. Integr. Circuits Syst.*]{} [**20**]{} 1–21

Donetti L, Hurtado P I and Munoz M A 2005 [*Phys. Rev. Lett.*]{} [**95**]{} 4

Arenas A, D[í]{}az-Guilera A, Kurths J, Moreno Y and Zhou C S 2008 [*Physics Reports*]{} [**469**]{} 93–153

Boyd S and Vandenberghe L 2004 [*Convex Optimization*]{} (New York: Cambridge University Press)

Taylor J A 2015 [*Convex Optimization of Power Systems*]{} (Cambridge University Press)

Low S 2014 [*IEEE Transactions on Control of Network Systems*]{} [**1**]{} 15–27

Low S 2014 [*IEEE Transactions on Control of Network Systems*]{} [**1**]{} 177–189

Lavaei J and Low S H 2012 [*IEEE Trans. Power Syst.*]{} [**27**]{} 92–107

Kelly F and Yudovina E 2014 [*Stochastic Networks*]{} (Cambridge University Press)

Bertsimas D, Farias V F and Trichakis N 2011 [*Oper. Res.*]{} [**59**]{} 17–31

Luss H 2012 [*Equitable resources allocation: Models, Algorithms and Applications*]{} (New Jersey: John Wiley & Sons)

Bertsekas D P and Gallager R 1992 [*Data Networks*]{} (Prentice Hall)

Kelly F P, Maulloo A K and Tan D K H 1998 [*Journal of the Operational Research Society*]{} [**49**]{} 237–252

Tan D K H 1999 [*Mathematical Models of Rate Control for Communication Networks*]{} Ph.D. thesis Statistical Laboratory, University of Cambridge

Srikant R 2003 [*The Mathematics of Internet Congestion Control*]{} (Boston: Birkh[ä]{}user)

Carvalho R, Buzna L, Just W, Helbing D and Arrowsmith D K 2012 [*Phys. Rev. E*]{} [**85**]{} 14

Carvalho R, Buzna L, Bono F, Masera M, Arrowsmith D K and Helbing D 2014 [ *PLoS ONE*]{} [**9**]{} e90265

Kersting W H 2001 [*Distribution System Modeling and Analysis*]{} (CRC Press)

Sojoudi S and Low S H 2011 [*IEEE Power and Energy Society General Meeting*]{}

Gan L, Li N, Topcu U and Low S H to appear, 2014 [*IEEE Trans. on Automatic Control*]{}

Bose S 2014 [*An integrated design approach to power systems: from power flows to electricity markets*]{} Ph.D. thesis Electrical Engineering, California Institute of Technology

Lavaei J, Tse D and Zhang B S 2014 [*IEEE Trans. Power Syst.*]{} [**29**]{} 572–583

Strang G 2006 [*Linear Algebra and its Applications*]{} (WB Saunders)

M. S. Andersen, J. Dahl, and L. Vandenberghe [*CVXOPT: A Python package for convex optimization*]{}, version 1.1.6

Arenas A, Diaz-Guilera A and Guimera R 2001 [*Phys. Rev. Lett.*]{} [**86**]{} 3196–3199

Ullah A and Giles D E A 1998 [*Handbook of Applied Economic Statistics*]{} (New York: CRC Press)

World Bank 2015 World development indicators Tech. Rep.

[^1]: We write $\pitchfork(r)$ instead of $\pitchfork(\mathcal{V},\mathcal{E},r)$ to simplify the notation.

[^2]: We say that a power allocation $\{ P_l \}$ for $l=1,\ldots,N$ is better than another $\{ P^\prime_l \}$ if $P_l\ge P^\prime_l$ for all $l$, and for some $l$, $P_l>P^\prime_l$. A power allocation is *Pareto optimal* or efficient if there is no better power allocation. The Pareto frontier of a set is the set of all Pareto optimal points.

[^3]: <https://www1.ethz.ch/id/services/list/comp_zentral/cluster/index_EN>
