Critical behaviour in charging of electric vehicles

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 [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 h [2]. The time to fully charge the battery of an electric vehicle at home currently varies from 18 h (Level 1, in the United States at 110 V and 15 A with a charge power of 1.4 kW) to 4 h (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 min [3]. 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 [4–7]. 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 prioritize 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 [8, 9], and physicists have helped to increase our understanding of its synchronization [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, 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 [16, 17], 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.

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 e.g., [18–21]) 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 [22]. The burgeoning field of convex optimization in electricity networks [23–25] 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 [24–26]. 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 [22]. 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 linearizations, 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, 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)\in {\mathcal{P}}(t)$ is a vector $P(t)=({P}_{l}(t):l=1,...,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 maximizes the sum of vehicle utilities [27]. This allocation acts as a network protocol that distributes network capacity among users, and solves the following problem:

$$\begin{eqnarray}&&\mathrm{maximize}\;\;\;\displaystyle \sum _{l=1}^{N(t)}{U}_{l}({P}_{l}(t))\end{eqnarray} \tag{ 1a }$$

$$\begin{eqnarray}&&\mathrm{subject}\;\;\mathrm{to}\quad P(t)\in {\mathcal{P}}(t).\end{eqnarray} \tag{ 1b }$$

Here we explore two user utility functions. First, we consider the non-unique max-flow allocations given by ${U}_{l}({P}_{l}(t))={P}_{l}(t)$, i.e. we maximize the instantaneous aggregate power sent from the root node to the vehicles, which is a benchmark of efficient network throughput [28]. 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 maximizes the sum of the logarithm of user rates, that is ${U}_{l}({P}_{l}(t))=\mathrm{log}({P}_{l}(t))$. The proportional fairness allocation is especial, because the users and the network operator simultaneously maximize their utility functions [27]. Furthermore, the problem is convex, and so can be solved in polynomial time [22], and it can be naturally extended by adding positive weights to each term in the objective function equation (1a), to account for diversity in user demand or for more than one user at some nodes [27]. For the compact and convex set ${\mathcal{P}}(t)$, it can be shown that the allocation ${P}^{\mathrm{PF}}(t)$ that maximizes equation (1a), satisfies [27, 29]:

$$\begin{eqnarray}&&\displaystyle \sum _{l=1}^{N(t)}\displaystyle \frac{{P}_{l}(t)-{P}_{l}^{\mathrm{PF}}(t)}{{P}_{l}^{\mathrm{PF}}(t)}\leqslant 0.\end{eqnarray} \tag{ 2 }$$

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, equation (2) implies that to increase the instantaneous power allocated to a vehicle by a percentage , we have to decrease a set of other power allocations, such that the sum of the percentage decreases is larger or equal to . 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 that fairness and congestion control are two sides of the same coin, because the existing algorithms for fair allocations also manage network congestion [27, 30–35]. 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 [27, section 7.2].

The simplest flow model in electricity networks is the dc power flow, which is a linearization 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 [36]. 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 ⁵ . Edge ${e}_{{ij}}\in {\mathcal{E}}$ connects node i to node j, where i is closer to the root than j, and is characterized by the impedance ${Z}_{{ij}}={R}_{{ij}}+{\rm{i}}{X}_{{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)+{\rm{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 [37] —see figure 1(a). The vector V(t) denotes the voltage allocated to the nodes. The voltage drop ${\rm{\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 [36]. 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 figure 1. Kirchhoff's voltage law applied to the circuit in figure 1(b) yields (see appendix A):

$$\begin{eqnarray}&&{V}_{i}(t){V}_{j}(t)-{V}_{j}^{2}(t)-{P}_{\pitchfork(j)}(t){R}_{{ij}}-{Q}_{\pitchfork(j)}(t){X}_{{ij}}=0.\end{eqnarray} \tag{ 3 }$$

Zoom In Zoom Out Reset image size

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, and 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 maximizes 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 equation (3), and node voltages are within $((1-\alpha ){V}_{\mathrm{nominal}},(1+\alpha ){V}_{\mathrm{nominal}})$ for $\alpha \in (0,1)$, with $\alpha =0.1$ typically [36]. Thus, to find the max-flow allocation of power to the vehicles, we solve the optimization problem for fixed t:

$$\begin{eqnarray}&&\underset{V(t)}{\mathrm{maximize}}\quad U(t)=\displaystyle \sum _{l=1}^{N(t)}{P}_{l}(t)\end{eqnarray} \tag{ 4a }$$

$$\begin{eqnarray}&&{\rm{subject}}\;{\rm{to}}\quad (1-\alpha ){V}_{\mathrm{nominal}}\leqslant {V}_{i}(t)\leqslant (1+\alpha ){V}_{\mathrm{nominal}},i\in {\mathcal{V}}\end{eqnarray} \tag{ 4b }$$

$$\begin{eqnarray}&&{V}_{i}(t){V}_{j}(t)-{V}_{j}(t){V}_{j}(t)-{P}_{\pitchfork(j)}(t){R}_{{ij}}-{Q}_{\pitchfork(j)}(t){X}_{{ij}}=0,{e}_{{ij}}\in {\mathcal{E}}.\end{eqnarray} \tag{ 4c }$$

Constraint (4b) sets the limits on the node voltage. Equation (4c) is the physical law coupling voltage to power, generalized from equation (3) 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 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 corresponds to the 2 × 2 principal submatrix $W({e}_{{ij}},t)$ defined by [38, 39]:

$$\begin{eqnarray}&&W({e}_{{ij}},t)=(\begin{array}{c}{V}_{i}(t)\\ {V}_{j}(t)\end{array}){V}_{i}(t){V}_{j}(t)=(\begin{array}{cc}{V}_{i}^{2}(t) & {V}_{i}(t){V}_{j}(t)\\ {V}_{j}(t){V}_{i}(t) & {V}_{j}^{2}(t)\end{array})=(\begin{array}{cc}{W}_{{ii}}(t) & {W}_{{ij}}(t)\\ {W}_{{ji}}(t) & {W}_{{jj}}(t)\end{array}),\end{eqnarray} \tag{ 5 }$$

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 ${{vv}}^{{\rm{T}}}$. Hence, constraint (4c) 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 (4) is on the Pareto frontier ⁶ , since we maximize an increasing function in the objective. The rank one constraint is nonconvex, but it does not change the Pareto frontier or the optimum [25, 40], and we remove it to relax problem (4) to:

$$\begin{eqnarray}&&\underset{W(t)}{\mathrm{maximize}}\quad U(t)=\displaystyle \sum _{l=1}^{N(t)}{P}_{l}(t)\end{eqnarray} \tag{ 6a }$$

$$\begin{eqnarray}&&{\rm{subject}}\;{\rm{to}}\quad {((1-\alpha ){V}_{\mathrm{nominal}})}^{2}\leqslant {W}_{{ii}}(t)\leqslant {((1+\alpha ){V}_{\mathrm{nominal}})}^{2},\quad i\in {\mathcal{V}}\end{eqnarray} \tag{ 6b }$$

$$\begin{eqnarray}&&{W}_{{ij}}(t)-{W}_{{jj}}(t)-{P}_{\pitchfork(j)}(t){R}_{{ij}}-{Q}_{\pitchfork(j)}(t){X}_{{ij}}=0,\quad {e}_{{ij}}\in {\mathcal{E}}\end{eqnarray} \tag{ 6c }$$

$$\begin{eqnarray}&&W({e}_{{ij}},t)\;\succcurlyeq \;0,\quad {e}_{{ij}}\in {\mathcal{E}},\end{eqnarray} \tag{ 6d }$$

where the generalized inequality in constraint (6d) means the $W({e}_{{ij}},t)$ matrices are positive semidefinite [41].

The problem of allocating power to vehicles in a proportional fair way has the same constraints 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 maximize 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 ${{\rm{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:

$$\begin{eqnarray}&&{P}_{l}(t)=\displaystyle \frac{{{\rm{P}}}_{i}(t)}{{w}_{i}(t)},\end{eqnarray} \tag{ 7 }$$

where ${w}_{i}(t)={\displaystyle \sum }_{l=1}^{N(t)}{{\rm{\Delta }}}_{{il}}(t)$ is the number of electric vehicles charging on node i at time t, and ${{\rm{\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 C):

$$\begin{eqnarray}&&\underset{W(t)}{\mathrm{maximize}}\quad U(t)=\displaystyle \sum _{i\in {{\mathcal{V}}}^{+}}{w}_{i}(t){\mathrm{logP}}_{i}(t)\end{eqnarray} \tag{ 8a }$$

$$\begin{eqnarray}&&{\rm{subject}}\;{\rm{to}}\quad {((1-\alpha ){V}_{\mathrm{nominal}})}^{2}\leqslant {W}_{{ii}}(t)\leqslant {((1+\alpha ){V}_{\mathrm{nominal}})}^{2},\quad i\in {\mathcal{V}}\end{eqnarray} \tag{ 8b }$$

$$\begin{eqnarray}&&{W}_{{ij}}(t)-{W}_{{jj}}(t)-{P}_{\pitchfork(j)}(t){R}_{{ij}}-{Q}_{\pitchfork(j)}(t){X}_{{ij}}=0,\quad {e}_{{ij}}\in {\mathcal{E}}\end{eqnarray} \tag{ 8c }$$

$$\begin{eqnarray}&&W({e}_{{ij}},t)\;\succcurlyeq \;0,\quad {e}_{{ij}}\in {\mathcal{E}},\end{eqnarray} \tag{ 8d }$$

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 equation (7). The complexity of the problem (8) thus scales with the number $| {\mathcal{V}}|$ of nodes, which is typically smaller than the number N(t) of vehicles for large arrival rates λ. Similarly, we also aggregated vehicles in the implementation of problem (6), 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 λ) 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 (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. We simulated vehicles charging on the realistic SCE 47-bus and SCE 56-bus distribution networks [38], which are detailed in figure 2. To characterize 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 realizations 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 [42] on the ETHZ Brutus cluster ⁷ 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 λ. 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 h, but 4 min for $\lambda =0.05$ .

Zoom In Zoom Out Reset image size

We set the battery capacity B = 1 for all vehicles, and the nominal voltage ${V}_{\mathrm{nominal}}=1$. Scaling V_nominal by β, 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 (4) and (8) are invariant upon the scaling ${V}_{\mathrm{nominal}}^{\prime }=\beta {V}_{\mathrm{nominal}}$, ${V}_{i}^{\prime }(t)=\beta {V}_{i}(t)$ for all nodal voltages, ${P}_{l}^{\prime }(t)={\beta }^{2}{P}_{l}(t)$, ${P}_{\pitchfork}^{\prime }(t)={\beta }^{2}{P}_{\pitchfork}(t)$ and ${Q}_{\pitchfork}^{\prime }(t)={\beta }^{2}{Q}_{\pitchfork}(t)$, and ${B}^{\prime }={\beta }^{2}B$. Considering these scaling properties, our simulations can be extended to values of ${V}_{\mathrm{nominal}}\ne 1$, provided the vehicle capacity B is rescaled accordingly, and we use this property to rescale the problem when convenient.

We find critical behaviour that resembles results found in communication networks, in that both systems undergo a continuous phase transition [43]. 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 [43]:

$$\begin{eqnarray}&&\eta (\lambda )=\underset{t\to \infty }{\mathrm{lim}}\displaystyle \frac{1}{\lambda }\displaystyle \frac{\langle {\rm{\Delta }}N(t)\rangle }{{\rm{\Delta }}t},\end{eqnarray} \tag{ 9 }$$

where ${\rm{\Delta }}N(t)=N(t+{\rm{\Delta }}t)-N(t)$ and $\langle \ldots \rangle$ indicates an average over time windows of width ${\rm{\Delta }}t$. We calculate $\eta (\lambda )$ in the steady state, that is ${\mathrm{lim}}_{t\to \infty }{\rm{\Delta }}N(t)\propto {\rm{\Delta }}t$. For arrival rates $\lambda \lt {\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 \gt {\lambda }_{c}$ some vehicles have to wait for increasingly long times to fully charge (congested phase). The order parameter characterizes the phase transition: $\eta (\lambda )=0$ in the free-flow regime, and $\eta (\lambda )\gt 0$ in the congested phase, a higher order parameter meaning that queues of charging vehicles build up more rapidly.

Figure 3 is a plot of the order parameter for the 47- and 56-bus networks and the two congestion control methods, as a function of λ. Simulation results shown in figure 3(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 figures 3(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.

Zoom In Zoom Out Reset image size

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 figure 3. To overcome this limitation, we adopt the susceptibility-like function [43]:

$$\begin{eqnarray}&&\chi (\lambda )=\underset{{\rm{\Delta }}t\to \infty }{\mathrm{lim}}{\rm{\Delta }}t\ {\sigma }_{\eta }({\rm{\Delta }}t),\end{eqnarray} \tag{ 10 }$$

where ${\rm{\Delta }}t$ is the length of a time window, and ${\sigma }_{\eta }({\rm{\Delta }}t)$ is the standard deviation of the order parameter η. To compute $\chi (\lambda )$, we first consider a long time series and split it into windows with length ${\rm{\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 figure 4) , allowing us to study the dependencies of the critical arrival rates on the congestion control algorithm, as well as network topology and size.

Zoom In Zoom Out Reset image 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 figures 4(c) and (d) four representative instances of the time series of the number of vehicles charging on the 47-bus network at λ = 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 λ = 0.39 for this algorithm. In contrast, N(t) oscillates in proportional fairness, suggesting that the critical point is above λ = 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 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 (4b) 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 maximize the aggregate power allocated to vehicles that can take all instantaneous power they are allocated (elastic demand), on a network with bounded voltage drops (i.e. 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.

To characterize 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 [44]:

$$\begin{eqnarray}&&G=\displaystyle \frac{1}{2\mu }{\rm{E}}[| u-v| ]=\displaystyle \frac{1}{2\mu }{\displaystyle \int }_{0}^{\infty }{\displaystyle \int }_{0}^{\infty }| u-v| f(u)f(v)\ {\rm{d}}u\ {\rm{d}}v,\end{eqnarray} \tag{ 11 }$$

where u and v are independent identically distributed random variables with probability density f and mean μ. 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

$$\begin{eqnarray}&&\widehat{G}=\displaystyle \frac{{\displaystyle \sum }_{i=1}^{n}{\displaystyle \sum }_{j=1}^{n}| {x}_{i}-{x}_{j}| }{2{n}^{2}\mu }.\end{eqnarray} \tag{ 12 }$$

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 figure 5 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 λ. 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 [45]. 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.

Zoom In Zoom Out Reset image size

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 λ in proportional fairness than in max-flow. This result is surprising, because one would expect that, for a chosen arrival rate λ, the maximization of the aggregate instantaneous power would also lead to a maximization of the energy (the time integral of power), and hence to a maximization 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).

The angle θ between V_i and V_j is small in distribution networks [36] (see figure A1 ), 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 figure 1(b),

$$\begin{eqnarray}{\rm{\Delta }}{V}_{{ij}} & = & | {V}_{i}| -| {V}_{j}| \\ & \simeq & {V}_{i}-{V}_{j}\\ & = & {\mathfrak{R}}({I}_{{ij}}{Z}_{{ij}})={\mathfrak{R}}(\displaystyle \frac{{I}_{{ij}}{Z}_{{ij}}{V}_{j}^{*}}{{V}_{j}^{*}})={\mathfrak{R}}(\displaystyle \frac{{S}_{\pitchfork(j)}^{*}{Z}_{{ij}}}{{V}_{j}^{*}})\\ & = & {Re}(\displaystyle \frac{({P}_{\pitchfork(j)}-{\rm{i}}{Q}_{\pitchfork(j)})({R}_{{ij}}+{\rm{i}}{X}_{{ij}})}{{V}_{j}^{*}})\simeq \displaystyle \frac{{P}_{\pitchfork(j)}{R}_{{ij}}+{Q}_{\pitchfork(j)}{X}_{{ij}}}{{V}_{j}},\end{eqnarray} \tag{ A1 }$$

where the superscript asterisk denotes the complex conjugate transpose.

Zoom In Zoom Out Reset image size

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{eqnarray}&&{P}_{\pitchfork(k)}=\displaystyle \sum _{i\in {{\mathcal{V}}}_{\pitchfork(k)}}\displaystyle \sum _{l=1}^{N(t)}{{\rm{\Delta }}}_{{il}}(t){P}_{l}(t)+\displaystyle \sum _{i\in {{\mathcal{V}}}_{\pitchfork(k)}}\displaystyle \sum _{j:{e}_{{ij}}\in {{\mathcal{E}}}_{\pitchfork(k)}}{P}_{{ij}}(t),\end{eqnarray} \tag{ B1 }$$

and

$$\begin{eqnarray}&&{Q}_{\pitchfork(k)}=\displaystyle \sum _{i\in {{\mathcal{V}}}_{\pitchfork(k)}}\displaystyle \sum _{j:{e}_{{ij}}\in {{\mathcal{E}}}_{\pitchfork(k)}}{Q}_{{ij}}(t),\end{eqnarray} \tag{ B2 }$$

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{eqnarray}{S}_{{ij}}(t) & = & {P}_{{ij}}(t)+{\rm{i}}{Q}_{{ij}}(t)=({V}_{i}(t)-{V}_{j}(t)){I}^{*}\\ & = & ({V}_{i}(t)-{V}_{j}(t)){(\displaystyle \frac{{V}_{i}(t)-{V}_{j}(t)}{{R}_{{ij}}+{{iX}}_{{ij}}})}^{*}\\ & = & ({V}_{i}(t)-{V}_{j}(t)){({V}_{i}(t)-{V}_{j}(t))}^{*}{(\displaystyle \frac{{R}_{{ij}}-{\rm{i}}{X}_{{ij}}}{{R}_{{ij}}^{2}+{X}_{{ij}}^{2}})}^{*}\\ & = & ({W}_{{ii}}(t)-{W}_{{ij}}(t)-{W}_{{ji}}(t)+{W}_{{jj}}(t))\displaystyle \frac{{R}_{{ij}}+{\rm{i}}{X}_{{ij}}}{{R}_{{ij}}^{2}+{X}_{{ij}}^{2}}.\end{eqnarray} \tag{ B3 }$$

Since, the voltages are real, ${W}_{{ij}}(t)={W}_{{ji}}(t)$, and thus

$$\begin{eqnarray}&&{P}_{{ij}}(t)=({W}_{{ii}}(t)-2{W}_{{ij}}(t)+{W}_{{jj}}(t))\displaystyle \frac{{R}_{{ij}}}{{R}_{{ij}}^{2}+{X}_{{ij}}^{2}},\end{eqnarray} \tag{ B4 }$$

and

$$\begin{eqnarray}&&{Q}_{{ij}}(t)=({W}_{{ii}}(t)-2{W}_{{ij}}(t)+{W}_{{jj}}(t))\displaystyle \frac{{X}_{{ij}}}{{R}_{{ij}}^{2}+{X}_{{ij}}^{2}}.\end{eqnarray} \tag{ B5 }$$

In proportional fairness, we maximize the sum of the logarithm of the instantaneous power allocated to electric vehicles:

$$\begin{eqnarray}&&{\displaystyle \sum }_{l=1}^{N(t)}\mathrm{log}{P}_{l}(t)={\displaystyle \sum }_{i\in {{\mathcal{V}}}^{+}}{\displaystyle \sum }_{l=1}^{N(t)}{{\rm{\Delta }}}_{{il}}(t)\mathrm{log}\displaystyle \frac{{{\rm{P}}}_{i}(t)}{{\displaystyle \sum }_{l=1}^{N(t)}{{\rm{\Delta }}}_{{il}}(t)},\end{eqnarray} \tag{ C1 }$$

where P_l(t) is the instantaneous power allocated to electric vehicle l, and ${{\rm{P}}}_{i}$ the instantaneous power allocated to node i. To maximize equation (C1), 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 equation (C1), yielding:

$$\begin{eqnarray}&&U(t)=\displaystyle \sum _{i\in {{\mathcal{V}}}^{+}}\displaystyle \sum _{l=1}^{N(t)}{{\rm{\Delta }}}_{{il}}(t)\mathrm{log}\;{{\rm{P}}}_{i}(t)=\displaystyle \sum _{i\in {{\mathcal{V}}}^{+}}{w}_{i}(t)\mathrm{log}\;{{\rm{P}}}_{i}(t).\end{eqnarray} \tag{ C2 }$$