Congestion control in charging of electric vehicles

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 distribution networks. We analyse the inequality in the charging times as the vehicle arrival rate increases, and show that charging times are considerably more uneven in max-flow than in proportional fairness. We also analyse the onset of instability, and find that the critical arrival rate is indistinguishable between the two protocols.


I. INTRODUCTION
Electric vehicles may become competitive, in terms of total ownership costs, with internalcombustion 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 hours [2]. 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 [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 local distribution networks and requiring local infrastructure reinforcement [4][5][6][7].
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 synchronization [10] and stability [11,12] on the power grid. In parallel, recent advances in optimization and phase transitions [13,14] suggest that the tools of critical phenomena and optimization are converging, opening up new horizons for physicists. Here, we show how the mathematical toolbox of convex optimization [15] can provide important insights into the critical behaviour of a fleet of electric vehicles charging on a distribution network.
Local distribution networks are approximately trees, and thus are cycle-free. Hence, Kirchhoff's voltage law is never violated on network loops. 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 current law and keeping voltage drops bounded. Such congestion management requires a protocol so that vehicles charge in a coordinated way, and here we explore two such mechanisms: max-flow and proportional fairness. If the inter-arrival time becomes too short, the charging of electric vehicles takes too long and more cars arrive for charging than leave fully charged, and the system loses stability [16,17]. Hence, gaining insights into the onset of instability is of paramount importance for the design of future distribution networks. Further, controlling congestion on the network offers a natural way to coordinate charging strategies, and could pave the way for the design of decentralised congestion control algorithms [18].
The efficient solution to the problem of maximizing network flows is given by the non-unique allocations of max-flow. Such allocations maximize network throughput, however, they can also leave users with zero flow, which is considered unfair from the user point of view. It turns out that fairness and congestion control are two sides of the same coin [19][20][21][22][23][24][25]. 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 [25,Section 7.2]. Mathematically, the problem is to find the feasible allocation that maximizes the sum of the logarithm of user rates. The proportional fairness allocation is especial, because the users and the network operator simultaneously maximize their utility functions [25]. Furthermore, the problem is convex, and so can be solved in polynomial time [15].

II. MATERIALS AND METHODS
Consider a tree topology, such that electric vehicles can charge at each node. Vehicles arrive in continuous time (following a Poisson process with rate λ) with empty batteries, choose a node with uniform probability amongst all nodes (except for 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 in to a node, the network will allocate it an instantaneous power, which is a function of the network topology, electrical elements, as well as the state of charge and other electric vehicles' charging strategies.
We abstract the distribution network to a rooted directed tree ⋔ (r) with node (often called bus) set V, edge (also called branch) set E, and a root node r (feeder) that injects power into the tree [26]. Edge e i j ∈ E connects node i to node j, where i is closer to the root than j, and is characterised by the complex resistance Z i j = R i j + iX i j , where R i j is the edge resistance and X i j the edge reactance. The vector V(t) denotes the voltage allocated to the nodes. The power loss along edge e i j is given by S i j (t) = P i j (t)+iQ i j (t), where P i j (t) is the real power loss, and Q i j (t) the reactive power loss. Electric vehicle l receives active power [27] P l (t) until charged-see Fig. 1(a). The voltage drop ∆V i j down the edge e i j is the difference between the amplitude of the voltage phasors V i at node i and V j at node j [28]. Kirchhoff's voltage law applied to the circuit in Fig. 1 The power consumed by the subtree , where vehicles consume real power only, but network edges have both active (real) and reactive (imaginary) power losses. edge resistance and reactance is taken from [29]. Node 1 is the root node in both networks. Nodes 13, 17, 19, 23 and 24 of the SCE 47-bus network (in lighter colour) are photovoltaic generators, and we removed them from the network.

A. The optimization problem
Let N(t) be the number of electric vehicles in the network at time t. 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). Let ⋔ ( j) denote the subtree of the distribution network rooted in node j, with node set V ⋔( j) and edge set E ⋔( j) . Let P ⋔( j) denote the active power, and Q ⋔( j) the reactive power consumed by the subtree ⋔ ( j)-see Fig. 1. The maximum 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 Eq. (1), and node voltages are within (1 −α)100% of the nominal voltage (typically, α = 0.1) [28]. Thus, to find the max-flow allocation of power to the vehicles, we solve the optimization problem for fixed t: is quadratic, hence not convex in general, which implies that the problem is not accessible to polynomial time methods. To convexify problem (2), we define a weighted adjacency matrix W(t), such that edge e i j corresponds to the 2 × 2 principal submatrix W(e i j , t) defined by [29,30]: The matrices W(e i j , t) are positive semidefinite, because their eigenvalues (λ 1 = 0 and λ 2 = V 2 i +V 2 j ) are non-negative, and rank one, because they are of the form vv T . Hence, constraint (2c) can be replaced by three constraints: the first substitutes the quadratic terms in the voltages with linear terms in the W(e i j , t), and the second and third constraints guarantee that the W(e i j , t) are positive semidefinite and rank one.
The solution of problem (2) is on the Pareto frontier [31], since we maximize an increasing function in the objective. The rank one constraint is nonconvex, however, it does not change the Pareto frontier [32,33], and we remove it to relax the problem (2) to: where the generalized inequality in constraint (4d) means the W(e i j , t) matrices are positivesemidefinite [34].
The problem of allocating power to vehicles in a proportional fair way has the same constraints as problem (2), however, the objective function is the sum of the logarithm of the power allocated to the vehicles. 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 have equivalent instantaneous demand, and thus the power 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: where w i (t) = N(t) l=1 ∆ il (t) is the number of electric vehicles charging on node i at time t, and ∆ 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): 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. (5). The complexity of the problem (6) thus scales with the number V of nodes, which is typically smaller than the number N(t) of vehicles. Similarly, we also aggregated vehicles in the implementation of problem (4), but omit the proof.

B. The Gini Coefficient
The Gini coefficient is defined as [35] where u and v are i.i.d. 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, . . . , n), the empirical Gini coefficient, G, may be estimated by a sample mean The Gini coefficient is used as a measure of inequality, because a sample where the only non-zero value is x has µ = x/n and hence G = (n − 1)/n → 1 as n → ∞, whereas G = 0 when all data points have the same value.

III. RESULTS AND DISCUSSION
We simulated vehicles charging on the SCE 47-bus and SCE 56-bus networks [29], which are detailed in Fig. (2). 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 runs of vehicles charging, each simulation running for 15000 time units (we discretized each unit into 10 equal subintervals).
At each time step, we checked if the configuration of vehicles changed and if it did, we solved the max-flow problem (4) and the proportional fairness problem (6). We ran the simulations on the ETHZ Brutus cluster [36] due to the high computational requirements.
We set the vehicle battery capacity B = 1 for all vehicles, and V nominal = 1. Scaling the nominal voltage by β implies scaling the power delivered to vehicles by β 2 , for β ∈ (0, ∞). To see this, observe that problems (2) and (6)  runs and shaded areas represent 95% confidence intervals.
when the system is stable, vehicles will experience a slower increase in inequality in proportional fairness than in max-flow with the increase of the vehicle arrival rate λ. Figure 4 is a plot of the order parameter where ∆N(t) = N(t + ∆t) − N(t) and . . . indicates average over time windows of width ∆t. The order parameter represents the ratio between the number of uncharged vehicles and the number of vehicles that arrive at the network to be charged. The system is unstable when the order parameter is positive (for arrival rates λ > λ c ), and a higher order parameter means that queues are building up more rapidly. Numerical simulations suggest that λ c depends on several factors (the network In conclusion, we modeled the max-flow and proportional fairness protocols for the control of congestion caused by a fleet of vehicles charging on distribution networks. We analyzed the inequality in the charging times as the vehicle arrival rate increases, and showed that charging times are considerably more uneven in max-flow than in proportional fairness. We also analysed the onset of instability, and found that the critical arrival rate λ c is indistinguishable between the two protocols.
It is an open question whether λ c is the same for max-flow and proportional fairness.

ACKNOWLEDGMENTS
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 highperformance 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).

Appendix A: Voltage drop on one edge
The voltage drop ∆V i j down the edge e i j 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 [28]. The angle θ between V i and V j is small in distribution networks [28] (see Fig. 5), 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. 1(b), where the superscript asterisk denotes the complex conjugate transpose.
FIG. 5. The difference between the V i and V j phasors Z i j I i j , decomposed along the V j vector and its orthogonal direction. The phase angle θ difference between V i and V j is small, and hence the voltage drop can be approximated by ∆V i j ≃ ℜ I i j Z i j .

Appendix B: Active and reactive loads on a subtree
The active and reactive power consumed by the loads in the subtree rooted in node k can be computed as: and Q ⋔(k) = i∈V ⋔(k) j:e i j ∈E ⋔(k) where P i j (t) is the active and Q i j (t) the reactive power dissipated on a cable connecting nodes i and j. The complex power is given by: Since, the voltages are real, W i j (t) = W ji (t), and thus and In proportional fairness, we maximize the sum of the logarithm of the instantaneous power allocated to electric vehicles: where P l (t) is the instantaneous power allocated to electric vehicle l, and P i the instantaneous power allocated to node i. To maximize Eq. (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 Eq. (C1), yielding: