A Fluid Model of an Electric Vehicle Charging Network

We develop and analyze a measure-valued fluid model keeping track of parking and charging requirements of electric vehicles in a local distribution grid. We show how this model arises as an accumulation point of an appropriately scaled sequence of stochastic network models. The invariant point of the fluid model encodes the electrical characteristics of the network and the stochastic behavior of its users, and it is characterized, when it exists, by the solution of a so-called Alternating Current Optimal Power Flow (ACOPF) problem.


Introduction
To deal with the effects of climate change, many countries are in the process of implementing new policy measures to stimulate the usage of electricity generated by renewable sources such as solar and wind.This comes with many societal challenges as well as opportunities for research.The supply of energy is less predictable, which makes the task of keeping high-voltage transmission networks reliable more challenging.In the local distribution grids, new products and services that can be used to balance the grid emerge (such as smart devices), but also create more intermittency.
In particular, electric vehicles (EVs) can cause a substantial additional load on local distribution grids [16].This has led to a significant interest in the stochastic scheduling of electric vehicle networks.There are different ways to replenish the batteries of an EV.A stochastic network analysis of fastcharging networks has been performed in [35].The analysis of a stochastic network of battery swapping infrastructures is performed in [32].
The focus of the present paper is on analyzing congestion associated to slow charging, which happens when a car is parked while its owner is at home, at work, or shopping.In [10], it was suggested to model the evolution of slowly-charging electric vehicles in a local grid by bandwidth sharing networks, approximating the instantaneous allocation of electricity to vehicles by proportional fairness.The main constraint that needs to be satisfied is that the voltage drop in the network needs to remain bounded.The focus in [10] was solely on simulation, assuming a Markovian model and infinitely many parking spaces for EVs.Using simulations, the stability of proportional fairness and max-min fairness was examined.
In a recent paper [3], we proposed an extension of [10] by allowing for load limits, finitely many parking spaces, and deadlines (associated with parking times).The joint distribution of charging requirements and parking times was not restricted to Markovian or independence assumptions.Using heuristic arguments, [3] proposed a fluid model keeping track of the number of charged and uncharged cars in the system and an associated invariant point.This invariant point is shown to be computationally tractable in [3], as it is formulated in terms of an AC Optimal Power Flow problem with an exact convex relaxation.
The goal of the present paper is to put the analysis of [3] on rigorous footing using measure-valued fluid limits.As in [3], we allow the parking times and the charging requirements of electric vehicles to be dependent and generally distributed random variables.In addition, we consider general arrival processes with time-varying arrival rates and multiple electric vehicle types.The distribution grid is explicitly modeled and we allow for multiple parking lots, each with finitely many parking spaces.(The fluid approximation of [3] did not take the subtleties of dynamically rejecting vehicles at parking lots into consideration.) Our work is connected to the literature on bandwidth-sharing networks.Such networks have been successfully used to model communication networks where the set of feasible schedules is determined by the maximum amount of data a communication channel can carry per time unit [26].The stochastic analysis of bandwidth-sharing networks was initially restricted to specific networks [7,8].The application of fluid and diffusion approximations led to computationally tractable approximations of a large class of networks; see for example [9,20,29,30,33,34].
In the context of communication networks, proportional fairness is a non-trivial but justified approximation of the transmission control protocol (TCP).A similar justification in the context of EV charging is performed in [1,13].In these papers, by using arguments similar to the seminal work [22], it is shown how algorithms like proportional fairness emerge in decentralized EV charging.Our class of controls contain proportional fairness as a special case.
Our analysis is mostly related to [30].The main difference is that, in the setting of electric vehicles, an important constraint that needs to be satisfied is to keep the voltage drops bounded, making the bandwidth-sharing network proposed here different.This also causes new technical issues, as the capacity set can be non-polyhedral or even non-convex.In addition, arriving vehicles finding a full parking lot are discarded; we assume such cars park on a regular parking spot.This leads to the additional technical complication of a loss process in a measure-valued context.
We now describe our contributions in more detail.We develop a measure-valued fluid model for the vector process which describes the number of total and uncharged EVs in each parking lot, allowing the dynamics of the stochastic model to be approximated with a deterministic model.This model depends on the joint distribution of the charging requirements and the parking times.We show that our measurevalued fluid model arises as a weak limit of a vector of measure-valued processes under an appropriate scaling.Moreover, under an additional assumption on the network topology, we show that the invariant point of this dynamical system is unique and can be characterized in a computationally friendly manner by formulating a nontrivial AC optimal power flow problem (ACOPF), which is tractable as its convex relaxation is exact in many cases; as mentioned before, this characterization was observed and applied in [3].
In order to prove properties for the solutions of the fluid model, we investigate the properties of the bandwidth allocation function in our setting where the capacity set is convex.We establish similar continuity properties of the allocation function as in [29].While the structural properties of a linearized voltage model can be developed in full, for the AC power flow equations we were only able to show continuity of the allocation function.We conjecture that Lipschitz continuity and a monotonicity property hold as well, but leave these problems open; we refer to Sections 3 and 6 for more specific comments.
The rest of this paper is organized as follows.In Section 2, we provide a detailed model description.In particular, we introduce our stochastic model, the power flow models that we use, and we define the system dynamics.Next, in Section 3, we present a continuity property of the optimal power allocation.Then, we move to the stochastic model.A fluid model is presented in Section 4, where we also study its properties.Section 5 shows that the fluid model can arise as a weak limit of the fluid-scaled processes.In Section 6, we focus on the invariant analysis of the fluid model under an additional assumption on the network topology and conclude with a counterexample of monotonicity of general tree networks.All proofs are gathered in Sections 7-10.

Model description
In this section, we provide a detailed formulation of our model and explain various notational conventions that are used in the remainder of this work.The model description in this section is nearly identical to that in [3].We include all details on the network structure and physical characteristic for completeness; the main difference is that the measure-valued state descriptor is fully developed and analyzed in the present paper.To this end, we also require more sophisticated notation, which we introduce first.

Preliminaries
We introduce the notational conventions that will be used throughout the paper.All vectors and matrices are denoted by bold letters.Further, R is the set of real numbers, R + is the set of nonnegative real numbers, and N is the set of strictly positive integers.For two real numbers x and y, define x ∨ y := max{x, y}, x ∧ y := min{x, y}, and x + := x ∨ 0. For two vectors x, y ∈ R I , define the coordinate-wise product x • y := (x 1 y 1 , . . ., x I y I ) (i.e, the Hadamard product) and the maximum norm x := max Vector inequalities hold coordinate-wise, namely x > y implies that x i > y i for all i.Furthermore, I represents the identity matrix and e and e 0 are the vectors consisting of 1's and 0's, respectively, the dimensions of which are clear from the context.Also, e i is the vector whose i th element is 1 and the rest are all 0.
Let Y be a metric space.Let M(Y ) be the space of Randon measures (i.e., locally finite and inner regular measures) on Y , endowed with the Borel σ-algebra denoted by B(Y ).Further, M F (Y ) is the space of the finite nonnegative measures in M(Y ) equipped with the weak topology.We say that a sequence of measures µ n in M F (Y ) converges to µ in the weak topology and we write It is known that M F (Y ) k , d k is a separate and complete space [5]; i.e., a Polish space.When Y = R n + , we simplify the notation to M F .

Network and infrastructure
We consider the typical situation where a low-voltage distribution network has a tree structure.Thus, take a rooted tree G = (I, E), where I = {0, 1, . . ., I}, denotes its set of nodes (buses) and E is its set of directed edges, assuming that node 0 is the root node (known as "feeder ").Denote by ik ∈ E the edge that connects node i to node k, assuming that i is closer to the root node than k.Let I(k) and E(k) be the node and edge set of the subtree rooted in node k ∈ I.The active and reactive power consumed by the subtree (I(k), E(k)) are P I(k) and Q I(k) .The resistance, the reactance, and the active and reactive power losses along edge ik are denoted by r ik , x ik , L P ik , and L Q ik , respectively.Moreover, V i is the voltage at node i and V 0 is known.At any node, except for the root node, there is a charging station with K i > 0, i ∈ I \ {0}, parking spaces (each having an EV charger).Further, we assume that there are J = {1, . . ., J} different types of EVs indexed by j.

Stochastic model for EVs
Type-j EVs arrive at node i according to a counting process E ij (•) := {E ij (t), t ≥ 0}; i.e., E ij (t) is the number of EVs that arrive into the parking lot in the time interval (0, t].We assume that all E ij (•) are finite, nondecreasing processes with where λ ij (s) > 0 are integrable functions.Moreover, let ζ ijl denote the arrival time of the l th type-j EV at node i.If all spaces are occupied, a newly arriving EV does not enter the system, but is assumed to leave immediately.
We now turn to the model characteristics of the EVs.Each EV has a random charging requirement (counted in time) and a random parking time.These depend on the type of the EV and the location that it is parked (i.e., the node), but are independent between EVs.Our framework is general enough to distinguish between types.For example, we can classify types according to intervals of ratio of the charging requirement and parking time and/or according to the contract they have with the network provider.An EV leaves the system after its parking time expires.It may not be fully charged.If an EV finishes its charge, it remains at its parking space without consuming power until its parking time expires.EVs that have finished their charge are called "fully charged ".
Let B ijl and D ijl denote the charging requirement and the parking time of the l th EV of type-j at node i.In queueing terminology, these quantities are respectively called service requirements and deadlines.Moreover, we assume that the sequence {B ijl , D ijl , l ∈ N} is a sequence of i.i.d.copies of a random vector (B ij , D ij ) with distribution law F ij (A) = P ((B ij , D ij ) ∈ A) for any Borel set A ∈ B(R 2 + ).Further, for l = 1, . . ., Q ij (0) we denote by (B 0 ijl , D 0 ijl ) the residual charging requirement and the residual parking time of the initial population of type-j at node i.Moreover, we assume the probability density function (pdf) of the parking times f Dij (•) exists with f Dij (0) > 0 for any i, j ≥ 1.Note that it is possible to communicate an indication of the charging time B and the parking time D by the owner of an EV [2].The model is illustrated in Figure 1.

Charging control rule
An important part of our framework is the way we specify how the charging of EVs takes place.Let the number of uncharged vehicles (of all types and in all nodes) be given by the vector z ∈ [0, ∞) I×J ; i.e., z ij is the number of uncharged vehicles of type-j in node i.We assume the existence of a vector function p(z) = (p ij (z) : i ∈ I \ {0}, j ∈ J ) that specifies the instantaneous rate of power each uncharged vehicle receives.Moreover, we assume that this function is obtained by optimizing a "global" function.Specifically, for a type-j EV at node i, we associate a function u ij (•) which is strictly increasing and concave in R + , twice differentiable in (0, ∞) with lim x→0 u ij (x) = ∞.The charging rate p(z) is then determined by max p I i=1 J j=1 z ij u ij (p ij ) subject to a number of constraints that take into account physical limits on the charging of the batteries, load limits, and most importantly voltage drop constraints.An important example is the choice u ij (p ij ) = w ij log p ij , which is known as weighted proportional fairness.
We next introduce the physical constraints of the network.The maximum electric power that can be consumed in total by all cars is M i > 0 at node i.Each type-j EV can be charged at a rate that is at most equal to c max j .That is, (2.1) We refer to (2.1) as "load constraints".In addition, we impose voltage drop constraints.These constraints rely on the load flow model used.Two of these models that we consider are described next.

A simplified AC voltage model
We consider a simplification of the full AC power flow equations, based on the typical situation that voltage angle differences in distribution networks are negligible [23,Chapter 3].Under this assumption, Kirchhoff's law [24, Eq. 1] takes the form, for pk ∈ E, where p ∈ I denotes the unique parent of node k.The previous equations are non-linear.Applying the transformation Note that W ( pk ) are positive semidefinite matrices (denoted by W ( pk ) 0) of rank one.The active and reactive power consumed by the subtree (I(k), E(k)) are given by where by [10, Appendix B], Note that W kk are dependent on the vectors p and z.We sometimes write W kk (p, z) when we wish to emphasize the dependence.The function p(z) is given by max p,W for z ij > 0. If z ij = 0, then p ij = 0.In addition, 0 < υ k ≤ W 00 ≤ υ k are the voltage limits.Observe that the optimization problem (2.5) is non-convex and in general NP hard due to the rank-one constraints. (2.7) Note that the constraints W ( ik ) 0 are equivalent to W ii W kk − W 2 ik ≥ 0, since we consider W ii > 0 for any node i ≥ 1.In the sequel, we freely use both formulations.

Linearized Distflow model
Though the previous voltage model is tractable enough for a convex relaxation to be exact, it is rather complicated.Assuming that the active and reactive power losses on edges are small relative to the power flows, but now allowing the voltages to be complex numbers, we arrive at a linear approximation of the previous model, called the linearized (or simplified) Distflow model [4].In this case, the voltage magnitudes W lin kk := |V lin k | 2 have an analytic expression [24, Lemma 12]: where the P(k) is the unique path from the feeder to node k.
Remark 2.1.Note that W lin kk ≤ W 00 for all nodes k, as we assume that the nodes only consume power, and by [24,Lemma 12] we obtain W kk (p, z) ≤ W lin kk (p, z).That is, we can remove the constraints W kk (Λ) ≤ υ k from (2.5).
To derive the representation of the power allocation mechanism p(z) in this setting, one replaces the constraints in (2.5) by (2.1) and υ k ≤ W lin kk (p, z).

State descriptor
In this section, we introduce the dynamics that describe the evolution of the system.Specifically, we will now incorporate in the system dynamics all residual processes needed to obtain a Markovian system.Let Q ij (•) and Z ij (•) be non-negative discrete measures for i, j ≥ 1.The total number of type-j EVs at node i at time t > 0 and the number of uncharged EVs are given by gives the total number of EVs at node i ≥ 1. Recall that ζ ijl is the arrival time of the l th EV of type-j at node i.The residual parking time of the l th newly arriving EV of type-j can be written as D ijl (t) := (D ijl − (t − ζ ijl )) + , l = 1, . . ., E ij (t) and for the initial population D 0 ijl (t) := (D 0 ijl − t) + , l = 1, . . ., Q ij (0).In order to define the residual charging requirements, we first introduce the following operators: (2.9) For s ≤ t, S ij (Z, s, t) is the cumulative bandwidth allocated per type-j EV at node i during time interval [s, t].The residual charging requirement of the l th type-j EV at node i at time t ≥ 0 is given by for the newly arriving EVs and B 0 ijl (t) = (B 0 ijl −S ij (Z, 0, t)) + , l = 1, . . ., Z ij (0), for the initially uncharged EVs.Now, we define the measure-valued state descriptor for any t ≥ 0 and for any Borel set B ⊆ R, (2.10) The measure δ + • (B) is the Dirac measure restricted on (0, ∞); i.e., δ + x (B) := δ x (B∩(0, ∞)) and counts the total number of type-j EVs in node i whose residual parking time belongs to the Borel set B.
The number of uncharged EVs for which the minimum between their residual charging requirement and their residual parking time belongs to any Borel set B ⊆ R 2 is given by (2.11) represents the event that there is an idle EV charger right before the arrival of the l th type-j EV.As not all EVs enter the system, we naturally define the following stochastic processes.First, the number of accepted type-j EVs at node i until time t > 0 is given by (2.12) Next, the number of rejected EVs until time t > 0 is given by (2.13) Observe that the following relation holds: Having introduced the stochastic model, which is defined through equations (2.10)-(2.13),we move to the main results of this paper.We first study some properties of the bandwidth allocation function in Section 3. We then define an appropriate fluid model in Section 4 and derive some of its properties.

Continuity of the optimal allocation function
In this short section, we state some structural properties of the optimal allocation function, which may be of independent interest.In particular, we show that the optimal solution of (2.7) is continuous under the AC power flow model (2.3).This result is needed in Section 5 in order to show convergence of the fluid-scaled processes.Last, in power system analysis, rigorous proofs are typically difficult and require additional assumptions on the distribution system [12], even if one ignores the stochastic dynamics.In the rest of this section, we make the additional assumption that the ratio of resistance and reactance is constant for all the lines, i.e., r pl x pl remains constant for any pk ∈ E. We show that the optimal aggregated power allocation Λ(z), z ∈ (0, ∞) I×J , is a continuous function in z.In order to establish this property, we first present a preliminary result.
we have that Λ is also a feasible point of (2.7).
Observe that in case the feasible set of (2.7) is polyhedral, the conclusion is immediate.The proof of the previous proposition is given in Section 7. The main idea of the proof is to construct a new solution (W , Λ ).Then, using the feasibility of the point Λ(z) and induction starting from the leaf nodes, we show that the point (W , Λ ) lies in the feasible set of (2.7).In the sequel, we present the main result of this section, which says that Λ(•) is a continuous function.
The proof of Theorem 3.2 is given in Section 7 and it combines Proposition 3.1, the continuity property of the voltages as functions of loads, and arguments from [30, Lemma 1].
In Section 4, and more specifically when we prove that the fluid model solution is unique, we need the stronger property that the optimal solution of (2.7) is Lipschitz continuous.If we assume the linearized Distflow power model (see Section 2.4), then the feasible set of (2.7) is polyhedral, and hence Λ(•) is Lipschitz continuous by applying directly [29, Theorems 3.1 and 3.2].In the case of the AC power flow model, where the feasible set of (2.7) is convex, we need to make an additional assumption that the strict complimentary condition holds for some constraints before we can conclude Lipschitz continuity.While we have not been able to establish this property without the aforementioned assumption, we conjecture that Λ(•) is Lipschitz continuous, and leave this question open.
We now move to the original stochastic network and its fluid model.

Fluid model definition
In this section, we define and study the properties of a deterministic fluid model, associated with the stochastic model introduced in Section 2. All proofs of this section are gathered in Section 8. Define the following classes Further, for any A ∈ C and s ∈ R, define A + s := {s + y, y ∈ A} and for any A ∈ C and (s, where E ij (t) = t 0 λ ij (s)ds.We say that the vector Furthermore, for any t ≥ 0, A ∈ C, and A ∈ C the following relations hold Moreover, the functions and We call the vectors (Q(•), Z(•)) and (Q(•), Z(•)) the measure-valued fluid model solution and the numeric fluid model solution, respectively.
The fluid model equations, though still rather complicated, have an intuitive meaning.For instance, the term P B ij ≥ S ij (Z, s, t), D ij ≥ t − s resembles the fraction of EVs of type-j admitted to the system at time s at node i that are still in the system at time t.For this to happen, their deadline needs to exceed t − s and their service requirement needs to be bigger than the service allocated, which is S ij (Z, s, t).In addition, R ij (t) represents the lost fluid of type-j EVs at node i due to a full system until time t ≥ 0. We next show that the total number of EVs in the fluid model can be rewritten in a familiar form for queueing systems and the departure process in the fluid model can be written as a function of the total number of EVs.Proposition 4.1.We have that for any i ≥ 1 and j ≥ 1, where D ij (t) represents the amount of fluid that departs from the system in time interval [0, t), and The last proposition uses the assumption of existence of the density of the parking times in order to ensure that the limit in (4.4) exists, and this is the only point where we need this assumption.It follows from Proposition 4.1 that the total number of EVs can be written with the help of a one-dimensional reflection mapping.This result will be helpful, when we show uniqueness of the fluid model solution in Theorem 4.2.The novelty in our setting is (4.4),where an intuitive explanation is as follows., ∞)) represents the amount of fluid of type-j EVs at node i for which its residual parking time lies in the interval (0, ).It is natural now to expect that by dividing the last difference by > 0 and by allowing to be arbitrary small, the quantity lim the departure rate of an EV from the parking lot at time s > 0. Observe that (4.4) corresponds to [17, Equation 3.2].However, in the latter, the authors use different test functions to define the fluid model and they write the departure rate in terms of the hazard rate function.
Before we continue our analysis, we present an example in case of a Markovian model.We shall show that the departure process given in (4.4) can be written in the well-known form for a Markovian model [27], namely To show (4.5), use (4.2) and A = [ , ∞) in (4.1) to get Observing that lim An important question is when a solution of the fluid model equations (if it exists) is unique.The next theorem answer this question.
and consider the linearized Distflow power model 2.8.Suppose that Z(0) = 0 or that Z(0) ∈ (0, ∞) I×J and the first projection of Z(0) is Lipschitz continuous, i.e., there exists L > 0 such that for any i, j ≥ 1, x < x , and y > 0, Then there exists a unique solution of the fluid model equations.
The proof of Theorem 4.2 is given in Section 8 and the main steps of the proof are as follows.
1.The first step is to show that each pair Q i (•), R i (•) satisfies a one-dimensional reflection mapping and (each pair) is unique.
2. Second, we show that Z ij (t) > 0 for any i, j ≥ 1 and t > 0.
Having defined a fluid model and studied its properties, the next main step is to show that the fluid model arises as a weak limit of the original stochastic model under an appropriate scaling.This is the topic of the next section.

Fluid limit theorem
In this section, we study the asymptotic behavior of the stochastic network described in Section 2. Consider a family of systems indexed by n ∈ N, where n tends to infinity, with the same basic structure as that of the system described in Section 2. To indicate the position of the system in the sequence of systems, a superscript n will be appended to the system parameters and processes.
First, we introduce our asymptotic regime.We assume that the scaled capacity at node i is given by M n i = nM , the scaled number of EV chargers at node i is K n i = nK, and the scaled resistance and reactance on line pk are given by r n pk = r pk /n and x n pk = x pk /n.Note that in our setting we need to scale the physical parameters of the system in contrast to the typical scalings in stochastic networks that arise in communication networks.We summarize below the assumptions we make in this section. Assumptions: 1.The scaled parameters are given by K n i = nK, M n i = nM , r n pk = r pk /n, and x n pk = x pk /n.

The external arrival process satisfies
3. The limit of the external arrival process is Lipschitz continuous; i.e., there exists 4. The scaled initial configurations converge to random vectors of finite measures, are almost surely free of atoms.
Having introduced our scaling regime, we now move to the fluid-scaled state descriptor.The fluidscaled measure-valued processes are given by Q n and the fluid-scaled counting processes are given by Q . Moreover, our fluid scaling leads to the following relation p n (z) = p( z n ).To see the latter, observe that under our scaling the feasible set of (2.6) can be written as follows Furthermore, by (2.9), we have that The next theorem states that the fluid model arises as a limit of the fluid-scaled state descriptor under our assumptions.
Theorem 5.1 (Fluit limit).The sequence of the fluid-scaled measure-valued vector process Q is tight and every accumulation point When the allocation mechanism is given by the linearized Distflow power model 2.8, we can invoke Theorem 4.2 to strengthen this result to a convergence result.For the full AC case, the same can be concluded if the bandwidth allocation function is Lipschitz continuous.The proof of Theorem 5.1 is given in Section 9, which is organized as follows.
2. We then show tightness for the fluid-scaled stochastic process describing the number of rejected customers, i.e., R n (•).
3. The last step is to show that the limit of any convergent subsequence of Q n (•), Z n (•) satisfies the fluid model equations.
Remark 5.1.The fluid limit theorem holds even if the external arrival process is a process with a general mean E ij (•).In this case, we need to modify the definition of a fluid model solution such that However, it seems that the uniqueness of the fluid model solutions does not hold.
In this section, we have obtained a fluid limit that holds for a general tree network.In the next section, we investigate under what assumptions the fluid limit converges to an invariant point.

Invariant analysis
In this section, we study the behavior of the system as time goes to infinity.To do so, we assume that the arrival rate is constant, i.e., E ij (t) = λ ij t.
First, we prove a result equivalent to [17,Theorem 3.6].There, a characterization of the invariant point for a loss system is shown, which we now prove in our setting.However, the difference with [17] is that we consider different test functions to define the fluid model and second, we consider multiple types of customers.All proofs are gathered in Section 10.
Define the traffic intensity at node i of type-j EVs by ] and the total traffic intensity an node i by ρ i := J j=1 ρ ij .The following result characterizes the invariant states of a loss system with multiple types of EVs.Proposition 6.1.Let λ ij ∈ (0, ∞).We have that (Q * , q * ) is invariant if and only if for any Borel set A ∈ B(R + ) and i, j ≥ 1, and In the sequel, we examine the asymptotic behavior of (Z(•), Z(•)).We make an additional assumption that the network is monotone as it is stated in the following definition.Definition 6.1.An allocation mechanism is called "monotone" if 0 < y ≤ z implies that p ij (y) ≥ p ij (z).
For instance, this property holds when the network has a line topology under the linearized Distflow model described in Section 2.4.2.In this case, [9, Proposition 5] can be applied directly in order to show the desired monotonicity property.We conjecture that this monotonicity holds true for a line network under the AC power flow model as well, but have not been able to prove this, apart from the case of two nodes.Proposition 6.2.If the network is monotone, then we have that (Z(t), Z(t)) → (Z * , z * ) as t → ∞.Furthermore, the vector (Z * , z * ) satisfies the following relation.For any Borel set A ∈ B(R 2 + ) and i, j ≥ 1, and z * is given by the solution of the fixed-point equation The proof of the last proposition combines Proposition 6.1 and arguments from [30, Theorem 2].Moreover, using similar arguments from [30,Theorem 6] and [19,Theorem 3.3], it can be shown that the fluid and steady-state limits can be interchanged and the sequence of fluid-scaled stationary distributions (Z n (∞), Z n (∞)) converges weakly to the invariant point as n → ∞, provided that the invariant point is unique.
The invariant point z * , when it is unique, can be computed by solving a single ACOPF problem, which is convex in our case since it admits an exact convex relaxation.Define the functions and recall that the aggregated allocation (the total power which type-j EVs consume) at node i is Λ ij (z) := z ij p ij (z).Also, for a random variable Y , denote by inf(Y ) the leftmost point of its support.
Proposition 6.3 (Characterization of the invariant point).Let inf (D j /B j ) ≤ 1/c max j .The solution z * of (6.1) is unique and is given by z , where Λ * is the unique solution of the optimization problem By (2.4), observe that W kk depends on z through the products z ij p ij (z).By the definition of Λ, we have that W kk depends only on Λ.That is, the previous optimization problem is indeed independent of the fixed point z * .Note that when the assumption inf (D j /B j ) ≤ 1/c max j is violated, there can be a continuum of invariant fluid model solutions [30].To arrive at (6.3), the essential idea is to add Little's law (6.1) to the set of Karush-Kuhn-Tucker (KKT) conditions that characterize p(z) and rewrite all equations in such a way that they form the KKT conditions for the problem (6.3).The proof of Proposition 6.3 is given in [3, Theorem 1].
A natural question is if Proposition 6.2 holds for non-monotone networks.However, this is an open problem even in the area of communication networks; see [30] and [9].In bandwidth-sharing networks, the feasible set is polyhedral.It is further proved that if the network topology is radial, then the monotonicity property holds [9, Proposition 5].Unfortunately, this is not the case in our setting.Distribution networks are not in general monotone.More surprisingly the monotonicity property for the tree networks does not hold even if we consider a polyhedral feasible set, i.e., linearized Distflow model.In the next section, we explain the reason why the monotonicity property fails.

A counterexample of monotonicity for a general tree network
A line network is monotone under the linearized Distflow model as we have already discussed and we conjecture this property is true for the AC power flow model as well.However, if we extend the line network to a tree network, the monotonicity property may fail to hold for both power flow models.Below we present a counterexample for both of power flow models.
Assume a tree network with four nodes, i.e., one feeder and three load nodes.The feeder (node 0) is connected to node 1, which has two children, nodes 2 and 3.Moreover, assume that the network is constant in the sense that all the resistances and reactances are the same for all the lines and take r pk = x pk = 0.1 for any edge pk .Further, the voltage magnitude at the feeder is fixed and taken to be W 00 = 1 and the lower bound for the voltage magnitude is given by υ = (0.9) 2 = 0.81.We shall show numerically that monotonicity does not hold for this simple tree network using the proportional fairness allocation mechanism.Indeed, we solve the optimization problem for the vectors z = (1, 1, 1) and y = (1, 2, 1).The allocated power to cars is given in Table 1, where we observe that p 3 (y) > p 3 (z)  2. The intuition behind this counterexample is as follows.The voltage constraints at the leaf nodes are both active for vector z as the network is constant.The voltage magnitude constraint remains active in the node where we increase the number of EVs, i.e., node 2 in this example.The total allocated .However, the feeder can not allocate more power (than Λ 2 (y)) to node 2 because of the voltage drop constraints.As a result, the feeder allocates more power to node 3 and hence p 3 (y) > p 3 (z), even though the number of EVs at node 3 does not increase.To see this, we remove the voltage drop constraints from the model and solve again the optimization problem; see Tables 3 and 4. Observe now that p(z) ≥ p(y) and the voltage magnitude at node 2 decreases.Even more surprisingly, the same behavior holds even if we use the linearized Distflow model.To see that, observe that by (2.8), the voltages at leaf nodes have an explicit solution, namely In other words, the voltages at leaf nodes depend on the power which is allocated to all three nodes.This is not the case for the constraints in bandwidth-sharing networks studied in [9].Constraints like the above correspond to non-tree networks in their setting.Hence, the desired monotonicity property does not hold in general in our model.Our intuition agrees with the numerical results in Table 5.
As the number of nodes I + 1 is finite, there exists I ≤ I + 1 such that L I = {0}, i.e., L I contains only the feeder node.Note that L 0 is the set of leaf nodes and the family L := {L m , 0 ≤ m ≤ I } is a partition of the set I. Indeed, we have that ∅ / ∈ L, I m=0 L m = I, and L i ∩ L k = ∅ for i = k.In Figure 2, we depict an example of a partition with five sets.

Feeder
Figure 2: The sets L i in a tree network.In this case I = 4.The red nodes are in L 0 , the blue nodes are in L 1 , the yellow node is in L 2 , the green node is in L 3 , and L 4 includes only the feeder.
Without loss of generality, we consider a single type of EVs; otherwise set Λ k := J j=1 Λ kj .To simplify the notation, in the rest of the proof we write Λ instead of Λ(z) and W kk instead of W kk (Λ).Recalling that Λ is a feasible point of (2.7), we have that Clearly, Λ satisfies the linear constraints of (2.7 To show that Λ is a feasible point of (2.7), we need to construct W il , i, l ≥ 0 such that the additional constraints of (2.7) are satisfied if we replace Λ by Λ .To this end, set W 00 = W 00 , W pk = W pk , for pk ∈ E. Further, W kk for k ≥ 1, are given by the solution of We shall show that W kk ≤ W kk for k ∈ I.The proof is then concluded by observing that by the inequality W kk ≤ W kk , we have that υ k ≤ W kk (Λ ) for k ≥ 1.Furthermore, by the third equation of (7.1), we get for pk ∈ E, Thus, Λ satisfies all the constraints of (2.7), and hence it is a feasible point.We now proceed to the proof of the claim that W kk ≤ W kk for k ∈ I. Define a pkls := r pk r ls +x pk x ls r 2 ls +x 2 ls and I(k) − := I(k) \ {k}.For k ∈ L m for some 0 ≤ m < I we have that The last equation can be rewritten as follows a pkks (W ss − W ss ). (7.3) We now show the inequality W kk ≤ W kk for each k by induction.Let k ∈ L 0 .By (7.3), we have that where p is the unique parent of node k.If m = 1 (i.e., k ∈ L 1 ), then we have that I(k) − = I(k) \ {k} = L 0 ∩ I(k) \ {k} and { ls ∈ E(k) : l ∈ I(k) \ {k}} = ∅.Further, {s : ks ∈ E(k)} = L 0 ∩ I(k) \ {k}.By (7.3) and (7.4), we obtain Now, observe that where the last equation holds by the assumption that r pk x pk is constant.That is, W kk ≤ W kk , for k ∈ L 1 .Suppose now that k ∈ L 2 .By (7.3), we have that The last equation can be equivalently rewritten as follows Applying (7.5) in the last equation, we obtain the following relation Now, observe that using the assumption that r pk x pk is the same for all edges, we have that r pk −r ks a pkks = 0. Further, we have that Thus, recalling that Λ k − Λ k ≤ 0 for any k ∈ I, we derive that W kk − W kk ≤ 0 for k ∈ L 2 .Suppose now that for all k ∈ L j , j = 0, . . ., m, We shall show that the same holds for k ∈ L m+1 .To this end, by (7.3) and (7.6), we have that Using again the assumption that r pk x pk is the same for all lines, we obtain This concludes the proof.
Proof of Theorem 3.2.We follow the argument in [29,Lemma 7.1].Take a sequence z k ∈ (0, ∞) I×J such that z k → z as k → ∞.We proceed by contradiction.Let us assume that Λ(•) is not continuous at point z.That is Λ(z k ) → Λ and Λ = Λ(z).Note the limit Λ exists as the sequence Λ(z k ) lives in a subset of the compact set {Λ ∈ [0, ∞) I×J : Λ ≤ M }.First, we show that Λ is a feasible point of (2.7).As Λ(z k ) is the optimal solution of (2.7), replacing z by z k we have that Taking the limit as k → ∞, we derive )) ≥ 0 (as we assume υ i > 0).Now, by continuity of the voltage magnitudes [12,Theorem 3] we obtain W ii (Λ ) ≥ υ i and W ii (Λ )W ll (Λ ) − W il (Λ ) 2 ≥ 0. That is, Λ is a feasible point of (2.7).Recalling that Λ(z) is the optimal solution of (2.7), we have that To derive the contradiction we construct a point Λ k which is feasible for (2.7) if we replace z by z k .To this end, define for any k ≥ 1, , by Proposition 3.1, we have that Λ k is a feasible point of (2.7) by replacing z by z k for k ≥ 1.It follows that as k → ∞, and That is, by (7.7) there exists a sufficiently large k such that The last inequality yields a contradiction as Λ(z k ) is the optimal solution of (2.7) by replacing z by z k .

Proofs for Section 4
Proof of Proposition 4.1.Using the identity P (D ij < t) + P (D ij ≥ t) = 1, (4.2) can be written as where In the sequel, we show that D ij (t) can be written as in (4.4).By the definition of the fluid model, we have that Observing that we have that By the assumption of existence of the pdf f Dij (•), we have that Dividing the last equation by and letting go to zero, we have that In other words, the limit of the left-hand side of (8.2) exists.Integrating (8.2) from 0 to t and interchanging the integrals by using Tonelli's theorem [31], we derive Furthermore, the following inequality holds for any t ≥ 0, That is, D ij (t) represents the departure process which proves (4.3) and (4.4).
The first step to prove Theorem 4.2 is to show that the fluid model solutions are bounded away from zero.This is stated in the following proposition.Proposition 8.1.Under the assumptions of Theorem 4.2, we have that for any > 0, inf Further, by our assumptions there exists the probability density function of parking times and f Dij (0) > 0 for any i, j ≥ 1.It is enough to show that Z ij (•) remains positive when the system is not full.Assume that Z(0) = 0 and define τ =: {s > 0 : and this covers also the case that τ = ∞.Note that if the arrival rate is constant then the last bound coincides with the one in [30, Lemma 3].Now, for t > τ , we have that Q i (t) = K i and by the continuity of the fluid model solutions, we have that . Further, by (4.3), we have that and using (4.4), we obtain where we define δ ij (s) := lim , and the fact that R ij (s) = 0 for s ≤ τ we have that δ ij (τ ) > 0. Hence, Proof of Theorem 4.2.We first show that each pair ) is unique for any i ≥ 1.Note that by Remark 4.1, fluid model solutions are invariant with respect to time shifts, and hence it suffices to show that ( By Proposition 4.1, we have that where . Now, by the one-dimensional reflection mapping [11, Chapter 6], we have that where It is known that the reflection mapping Ψ(•) is Lipschitz continuous [11].Now, for each i ≥ 1, define the mapping B i for each function a(•) on [0, ∞), where Observing that By [17,Lemma 3], the following functional equation for any i ≥ 1 has a unique solution on [0, T ]: The main idea now is to show that each function R i (•) satisfies (8.5), and hence it is unique.To this end, by the proof of Proposition 4.1, the relation J h=1 λ ih (s) dR i (s), and the properties of the Riemann-Stieltjes integral, we obtain Using the last equation and replacing D ij (t) in (8.3), we have that Using again the reflection mapping, we obtain The last equation and (8.4) yield Combining (8.3) and (8.4), we derive Now, replacing Φ i (•) in the last equation by the right hand side of (8.6) leads to Thus, R i (•) is a solution of (8.5), and hence unique.This implies that R ij (•) is unique for any i, j ≥ 1 and hence, ( We now proceed to show the uniqueness of the Z ij (•).First, we show that Z ij (•) has a Lipschitz continuous first projection.Indeed, let x < x and y ≥ 0. For any i, j ≥ 0, we have that By the Lipschitz continuity of E ij (•), the previous bound becomes By the assumption of the Lipschitz continuity of the initial condition, the change of variable v = Θ(s) = S ij (Z, s, t), and [30, Lemma 5], we have that That is, the first projection of Z ij (•) is Lipschitz continuous with constant where the last inequality follows by Theorem 3.2.Note now that point 0 is a feasible point of (2.5).Further, for a vector z such that z ij is small enough the power flow constraints are satisfied and hence p ij (z) = c max j .Moreover, the power allocation function is Lipschitz continuous since we consider the linearized Distflow power flow model as we discussed in Section 3. Now, by the Lipschitz continuity of E ij (•) and by applying [30, Theorem 1], we obtain that the fluid model solution (Z(•), Z(•)) is unique.9 Proof of fluid limit Theorem 5.1

Establishing tightness
The first step of the proof of Theorem 5.1 is to show that Q n (•), Z n (•) is C-tight, i.e., tight with continuous weak limits.To do so, we follow the idea of proof of [30,Theorem 5].First, we show that both processes satisfy the compact containment property.To this end, note that the following bounds hold almost surely and Moreover, by our assumptions, Hence, all the bounds in [30, Lemma 9] hold true for the measure-valued processes Q ij (•) and Z ij (•).That is, for any T > 0 and > 0, there exist compact sets and lim inf n→∞ Next, we shall show the oscillation control.To do so, we first show a preliminary result.Define Proposition 9.1.For any T > 0, δ > 0, and > 0, there exist α > 0 and b > 0 such that We shall show that lim inf Then, the result follows.Indeed, by (9.1), (9.2), we have that Now, Proposition 9.1 follows by using the last inequalities, (9.5)-(9.8),and [30, Lemma 12].

Fluid limits satisfy the fluid model solutions
Note that the total number of EVs can be written as follows and by [28] the latter is a weakly convergent sequence in (D[0, ∞), J 1 ) and hence it is tight.By Prokhorov's Theorem it is also relatively compact.That is, In other words, D n ij (•) is stochastically bounded and hence satisfies the first property of Kurtz's criteria.To show that it also satisfies the second property, we write Note that for any t ≥ 0 and n ≥ 1, Further, by continuity of the random variables ζ ijl and D ijl , we have that as δ → 0, Putting all the pieces together, Next, we show that the fluid limits are bounded away from zero.
For any δ > 0, there exist C δ > 0 and C δ > 0 such that almost surely inf Proof.First, we shall show that Q ij (•) is strictly positive.It is enough to show this inequality when the system is not full.Fix ∆ > δ.It is enough to show the result for t ∈ [δ, ∆].Define By our assumptions for the external arrival process, we have that for any m, where b is a continuity point for the distribution F Dij (•) with P (D ij > b) > 0, and the last inequality fol- Then, for large enough n, we have that for any i, j ≥ 1, Further, note that by continuity of the limit we have . Finally, we have that there exists C δ > 0 such that, for any ∆ > δ, as n → ∞.For any compact set C ⊆ R + , define the mapping φ C : D(R + , R I×J ) → R, given by φ C (y) := inf t∈C min i,j y ij (t).Note that φ C (y) is continuous at continuous y(•), which implies that By the Portmanteau theorem [6, Theorem 2.1], we have that We now move to the proof of Z ij (t) > 0 for t > 0. We first note that Q(•) is independent of Z(•), and hence we can assume that the fluid limit (Q(•), Q(•)) satisfies the fluid model equations as we shall show later.That is, (Q(•), Q(•)) satisfies the equations in Proposition 4.1.By Proposition 9.3, we have that the fluid-scaled process that describes the number of accepted EVs given in (2.12) converges weakly to A(t) := E(t) − R(t).First, we show that A ij (t) is strictly increasing for any i, j ≥ 1.Let t 1 , t 2 ≥ 0 with 0 ≤ t 1 < t 2 .Assume that there exists a subinterval in [t 1 , t 2 ] such that the total queue length at node i is full.Without loss of generality, assume that there exists τ ∈ [t 1 , t 2 ] such that Q i (s) = K i for any s ∈ [τ, t 2 ].First, assume that τ > t 1 , then we have that If τ = t 1 , then by (4.3), (4.4), and the fact that as n → ∞.Now, using again the Portmanteau theorem, we have that The total number of type-j EVs at node i can be written as follows Further, the above expression can be rewritten as where ξ ijl represents the time of the l th accepted EV and A n ij (•) represents the number of accepted type-j EVs at node i.In the same way, the number of uncharged type-j EVs at node i is given by In the sequel, we show that the fluid limit also satisfies the additional relations in Definition 4. We have proved that any subsequential limit (Q(•), Q(•), Z(•), Z(•), R(•)) satisfies the fluid model equations given in Definition 4.1, and hence the proof of Theorem 5.1 is completed.
10 Proofs for Section 6 Proof of Proposition 6.1.First assume that (Q * , q * ) is invariant, i.e., Q ij (t) = q * ij for any i, j ≥ 1 and t ≥ 0. We distinguish two cases i) J j=1 q * ij = K i and ii) J j=1 q * ij < K i .First, assume that J j=1 q * ij = K i .By (4.3) and (4.4), we have that Replacing (10.1) in (4.2) and taking the limit as time goes to infinity, we obtain Replacing the last equation in (4.1) and taking t → ∞, we have that for any Borel set A ∈ B(R + ), Last, by the nonnegativity of R ij (•), we obtain that ρ i > K i in this case.
In the second case, J j=1 q * ij < K i and by the fluid model equations R ij (t) = 0 for t ≥ 0 and for any i, j ≥ 1. Taking the limit as t → ∞ in (4.1), we have that for any Borel set A ∈ B(R + ), The "only if" part of the proposition follows using the same arguments as in the last part of proof in [17,Theorem 3.6].
Proof of Proposition 6.2.Uniqueness of the solution of the fixed-point equation follows by [30, Theorem 2] and [3].In order to study the asymptotic behavior of Z(•), let Q(0) = q * .By Proposition 6.1, we derive that . In We denote by C(Y, Y ) the space of continuous functions f : Y → Y and by C b (Y, Y ) the space of continuous and bounded functions f : Y → Y .By D(Y, Y ) denote the space of functions f : Y → Y that are right continuous with left limits endowed with the J 1 topology; i.e., the Skorokhod space.Further, we write X(•) := {X(t), t ≥ 0} to represent a stochastic process and X(∞) to represent a stochastic process in steady-state.Moreover, d = and d → denote equality and convergence in distribution (weak convergence).

Figure 1 :
Figure 1: A network with J types of EVs and constant arrival rates.

Definition 4 . 1 (
Fluid model).Let the initial data for the fluid model be given by

Remark 4 . 1 .
Note that the sets C and C generate the Borel σ−algebra of R and R 2 , respectively.Then, by Dynkin's π-λ theorem, the fluid model solutions hold for any Borel set.See Section 2.3 in[14] for more details.Moreover, by[30, Remark 3.2], fluid model solutions are invariant with respect to time shifts.

Example 4 . 1 (
Markovian model).Consider a Markovian model (i.e., Poisson arrival process with constant arrival rate and exponential parking times), and take J = 1 and Q i (0) = 0 for convenience.

7 Proofs for Section 3 Proof of Proposition 3 . 1 .
First, note that the point 0 lies in the feasible set by choosing W pk = W 00 .We now define a partition of the set I. Recall that I(k) denotes the subtree rooted in node k ∈ I (including node k).Let us define the following sets L 0 := {k ∈ I : I(k) = {k}} and for any m ≥ 1,

t+δ) l=1 1
{t<ζ ijl +D ijl ≤t+δ} .By (9.14), the fact that D n,∞ ij (•) is relatively compact and using the same arguments as in [21, Lemma 5.10], we conclude that D n ij (•) satisfies the second property of Kurtz's criteria.That is, D n ij (•) is relatively compact and hence tight.The tightness of R n (•) follows by (9.13) and by the tightness of D n (•) and

1 .λλ
Observe that by(2.13)  and the definition of the Riemann-Stieltjes integral, we have that ij (s)1 {Q n i (s − )=Ki} ds, ij (u)du .By our assumptions for the arrival process, we obtain thatR n i (•) − H n i (•) d → 0 as n → ∞, and hence H n i (•) d → R i (•).Now, by(2.13), the number of rejected type-j EVs at node i can be written as follows of the external arrival process and the fact that H

0 P
λij ρi (ρ i ∧ K i )t.Now, by the definition of the fluid model (Definition 4.1) we have thatZ ij (t) = Z ij (0)P B 0 ij ≥ S ij (z, 0, t), D 0 ij ≥ t + λ ij ρ i (ρ i ∧ K i ) t (B ij ≥ S ij (z, s, t), D ij ≥ t − s) ds.By [30, Theorem 3], there exist b u , b l ∈ (0, ∞) I×J such that0 < b l ij ≤ lim inf t→∞ Z ij (t) ≤ lim sup t→∞ Z ij (t) ≤ b u ij .(10.3) Furthermore, b u , b l satisfy the following relations b Now, we have assumed that the network is monotone, and hence sup b l ≤z≤b u p ij (z) = p ij (b l ) and inf b l ≤z≤b u p ij (z) = p ij (b u ).Applying the last relation in (10.3) and (10.4), we have that b other words, b u , b l satisfy the fixed-point equation (6.1) and by the uniqueness of the solution of the fixed-point equation, we obtain b u = b l = z * and hence lim t→∞ Z ij (t) = z * ij .

Table 5 :
Allocated power to EVs for the linearized Distflow model [30,her, by the proof of Proposition 8.1, we have that δ ij (s) =δ ij (t 1 ) > 0 for s ∈ [t 1 , t 2 ], and hence A ij (t 2 ) − A ij (t 1 ) > 0. Now, consider a type-j EV l at node i.Observe that by the constraints p ij (•) ≤ c max ∧ D ijl after its arrival.Hence, the stochastic process Z ij (•) is bounded from below by the queue length Q inf ij (•) of the infinite-server queue with arrival process A ij (•), Q inf ij (0) = 0, and i.i.d.service requirements { ∧ D ijl , l ∈ N}.Recalling that A ij (•) is strictly increasing by[30,   Lemma 3.14], there exists C δ > 0 such that, for any ∆ > δ, n ij (t) ≥ C δ ≥ P n inf δ≤t≤∆ min i,j 0) (A + (S ij (Z n , 0, t), t)) have that S ij (Z n , s, t) → S ij (Z, s, t), for s ∈ [t 0 , t] as n → ∞.Moreover, the function S ij (Z,s, t) is continuous and S ij (Z n , s, t) is monotone in s.Hence, we have that sup t0≤s≤t S ij (Z n , s, t) − S ij (Z, s, t) → 0. Now, the convergence follows by adapting the conclusion of the proof of [30, Theorem 5, Section 7.6].