Bounds and Limit Theorems for a Layered Queueing Model in Electric Vehicle Charging

The rise of electric vehicles (EVs) is unstoppable due to factors such as the decreasing cost of batteries and various policy decisions. These vehicles need to be charged and will therefore cause congestion in local distribution grids in the future. Motivated by this, we consider a charging station with finitely many parking spaces, in which electric vehicles arrive in order to get charged. An EV has a random parking time and a random charging time. Both the charging rate per vehicle and the charging rate possible for the station are assumed to be limited. Thus, the charging rate of uncharged EVs depends on the number of cars charging simultaneously. This model leads to a layered queueing network in which parking spaces with EV chargers have a dual role, of a server (to cars) and customers (to the grid). We are interested in the performance of the aforementioned model, focusing on the fraction of vehicles that get fully charged. To do so, we develop several bounds and asymptotic (fluid and diffusion) approximations for the vector process, which describes the total number of EVs and the number of not fully charged EVs in the charging station, and we compare these bounds and approximations with numerical outcomes.


Introduction
The rise of electric vehicles (EVs) is unstoppable due to factors such as the decreasing cost of batteries and various policy decisions [23]. Currently, the bottlenecks are the ability to charge a battery at a fast rate and the number of charging stations, but this bottleneck is expected to move towards the current grid infrastructure. This is illustrated in [22], the authors evaluate the impact of the energy transition on a real distribution grid in a field study, based on a scenario for the year 2025. The authors confront a local low-voltage grid with electrical vehicles and ovens and show that charging a small number of EVs is enough to burn a fuse. Additional evidence of congestion is reported in [11]. This paper proposes to model and analyze such congestion by the use of the so-called layered queueing networks. Layered networks are specific queueing networks where some entities in the system have a dual role; e.g., servers (in our context: parking spaces with EV chargers) become customers to a higher-layer (here: the power grid). The use of layered queueing networks allows us to analyze the interaction of two sources of congestion; first, the number of available spaces with charging stations (as not all cars find a space) and second, the amount of available power that the power grid is able to feed to the charging station. [22] We consider a charging station (or parking lot) with finitely many parking spaces. Each space has an EV charger connecting with the power grid. EVs arrive at the charging station randomly in order to get charged. If an EV finds an available space, it enters the parking lot and charging starts immediately.
An EV has a random parking time and a random charging time. It leaves the parking lot only when its parking time expires; i.e., it remains at its space without consuming power until its parking time expires if finishing its charge withing the parking time. Both the charging rate per vehicle and the charging rate possible for the complete charging station are assumed to be limited. Thus, the charging rate of uncharged EVs depends on the number of cars charging simultaneously. Last, we assume that all available power is shared at the same rate to all cars that need charging. The available power that can delivered by the grid is assumed to be constant.
Using queueing terminology, our model can be described as a two-layered queueing network. An EV enters the charging station and connects its battery to an EV charger. In our context, EVs play the role of customers, while EV chargers are the servers. Thus, the system of EVs and EV chargers can be viewed as the first layer. Moreover, EV chargers are connected to the power grid. Thus, at the second layer, active EV chargers act as jobs that are served simultaneously by the power grid, which plays the role of a single server.
This paper focuses on the performance analysis of this system under Markovian assumptions. Specifically, we are interested in finding the fraction of fully charged EVs in the charging station, which is equivalent to the probability that an EV leaves the charging station with a fully charged battery. A mostly heuristic description of some partial results in this paper has appeared in [2]. We first start with the steady-state analysis of the original system, for which we can find explicit bounds for the fraction of fully charged EVs. To do so, we study three special cases of the original system: i) there is enough power for all EVs, ii) there are enough parking spaces for all EVs, and iii) the parking lot is full. In these cases, we are able to find the explicit joint distribution in steady-state of the total number of EVs and the number of not fully charged EVs in the charging station, which we call the vector process.
In order to improve the bounds for the fraction of fully charged EVs, we next develop a fluid approximation for the number of uncharged EVs in the parking lot. The mathematical results here are closely related to results on processor-sharing queues with impatience [21]. However, the model here is more complicated as there is a limited number of spaces in the system and fully charged cars may not leave immediately as they are still parked.
We then move to diffusion approximations, working in three asymptotic regimes. First, we consider the Halfin-Whitt regime, in which we prove a limit theorem for the vector process, showing that it converges to a two-dimensional reflected Ornstein-Uhlenbeck (OU) process with piecewise linear drift. Then, we consider an overloaded regime for the process describing the number of total EVs in the system. In this case, the limit reduces to a one-dimensional OU process with piecewise linear drift. Last, we approximate the vector process by a two-dimensional OU process when the parking times are sufficiently large. The mathematical results here are based on martingales arguments [30].
EVs can be charged in several ways. Our setup can be seen as an example of slow charging, in which drivers typically park their EV and are not physically present during charging (but are busy shopping, working, sleeping, etc). For queueing models focusing on fast charging, we refer to [5,46]. Both papers consider a gradient scheduler to control delays. Next, [45] presents a queueing model for battery swapping while [38] is an early paper on a queueing analysis of EV charging, focusing on designing safe control rules (in term of voltage drops) with minimal communication overhead.
Despite being a relatively new topic, the engineering literature on EV charging is huge. We can only provide a small sample of the already vast but still emerging literature on EV charging. The focus of [37] is on a specific parking lot and presents an algorithm for optimally managing a large number of plug-in EVs. Algorithms to minimize the impact of plug-in EV charging on the distribution grid are proposed in [36]. In [29], the overall charging demand of plug-in EVs is considered. Mathematical models where vehicles communicate beforehand with the grid to convey information about their charging status are studied in [35]. In [28], EVs are the central object and a dynamic program is formulated that prescribes how EVs should charge their battery using price signals.
In addition, layered queueing network have been successfully applied in analyzing interactive networks in communications networks and manufacturing systems. These are queueing networks where some entities in the system have a dual role. In such systems, the dynamics in layers are correlated and the service speeds vary over time. Layered queueing networks can be characterized by separate layers (see [34] and [43]) or simultaneous layers (such as our model); [3]. In the first case, customers receive service with some delay. An application where layered networks with separate layers appear is the manufacturing systems e.g., [16] and [17]. On the other hand, in layered networks with simultaneous layers, customers receive service from the different layers simultaneously. Layered networks with simultaneous layers have applications in communications networks, e.g. in web-based multi-tiered system architectures. In such environments, different applications compete for access to shared infrastructure resources, both at the software level and at the hardware level. For background, see [39], and [40].
The paper is organized as follows. In Section 2, we provide a detailed model description -in particular we introduce our stochastic model and we define the system dynamics. Next, in Section 3 we present some explicit bounds in steady-state for the fraction of fully charged EVs. Section 4 contains several asymptotic approximations. First, a fluid approximation is presented; we then derive diffusion limits and approximations in three asymptotic regimes. Numerical validations are presented in Section 5. Last, all proofs are gathered in Section 6.

Model
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.

Preliminaries
We use the following notational conventions. 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 real numbers x and y, we define x ∨ y := max{x, y} and x ∧ y := min{x, y}. Furthermore, I represents the identity matrix and e and e 0 are 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 (Ω, F, P) be a probability space. For T ] be the twodimensional Skorokhod space; i.e., the space of two-dimensional real-valued functions on [0, T ] that are right continuous with left limits endowed with the J 1 topology; cf. [10]. Observe that as all candidate limit objects we consider are continuous, we only need to work with the uniform topology. It is wellknown that the space (D[0, T ] 2 , J 1 ) is a complete and separate metric space (i.e., a Polish metric space); [6]. We denote by B(D[0, T ] 2 ) the Borel σ−algebra of D[0, T ] 2 . We assume that all the processes are defined from (Ω, F, P) to (B(D[0, T ] 2 ), D[0, T ] 2 ). 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). For two random variables X, Y , we write X ≤ st Y (stochastic ordering) if P (X > a) ≤ P (Y > a) for any a ∈ R. Further, Φ(·) and φ(·) represent the cumulative probability function and the probability density function (pdf) of the standard Normal distribution, respectively. Last, let C 2 b (G) denote the space of twice continuously differentiable functions on G such that their first-and second-order derivatives are bounded.

Model description
We consider a charging station with K > 0 parking spaces. Each space has an EV charger which is connected to the power grid. EVs arrive independently at the charging station according to a Poisson process with rate λ. They have a random charging requirement and a random parking time denoted by B and D, respectively. The random variables B and D are assumed to be mutually independent and exponentially distributed with rates µ and ν, respectively. If an EV finishes its charge, it remains at its space without consuming power until its parking time expires. We call these EVs fully charged EVs. Thus, EVs leave the system only after their parking time expires, which implies that an EV may leave the system without its battery being fully charged. Furthermore, if all spaces are occupied, a newly arriving EV does not enter the system but leaves immediately. As such, the total number of vehicles in the system can be modeled by an Erlang loss system, though we need a more detailed description of the state space.
We denote by Q(t) ∈ {0, 1, . . . , K} the total number of EVs (charged and uncharged) in the system at time t ≥ 0, where Q(0) is the initial number of EVs. Further, we denote by Z(t) ∈ {0, 1, . . . , Q(t)} the number of EVs without a fully charged battery at time t and by Z(0) the number of such vehicles initially in the system. Thus, C(t) = Q(t) − Z(t) represents the number of EVs with a fully charged battery at time t.
The power consumed by the parking lot is limited and depends on the number of uncharged EVs at time t. We let it be given by the power allocation function L : We assume that the parameter M is given and that 0 < M ≤ K. For example, the parameter M can depend on the contract between the power grid and the charging station. Alternatively, M can be thought as the maximum number of EVs the charging station can charge at a maximum rate, where without loss of generality we can assume that the maximum rate is one. The model is illustrated in Figure 1.
Last, note that the processes Q(·), Z(·), and C(·) depend on K and M . We write Q K M (·), Z K M (·), and C K M (·), when we wish to emphasize this. It is clear from our context that the two-dimensional process {(Q(t), Z(t)) : t ≥ 0}, is Markov. The transitions rates in the interior and on the boundary are shown in Figure 2. Figure 2: Transition rates in the interior (left) and on the boundary (right) of the process {(Q(t), Z(t)) : t ≥ 0}.

Alternative model description in case of infinitely many parking spaces
Here, we give an alternative description of our model in case there are infinitely many parking spaces; i.e., K = ∞. In this case, the model can be described as a tandem queue with impatient customers; see Figure 3. EVs arrive at the charging station, which has M servers, and charging starts immediately. There are two possible scenarios. First, an EV gets fully charged during D and moves to the second queue, which has and infinite number of servers. This happens with rate µ(Z(t) ∧ M ). In the second queue, EVs get served with rate νC(t). In the second scenario, an EV abandons its charging because its parking time expired (and thus leaves the first queue impatiently); this happens with rate νZ(t). Note, that the total "rate in" in the system is λ and the total "rate out" is ν(Z(t) + C(t)) = νQ(t). In other words, Q(t) describes the number of customers is an M/M/∞ queue, i.e., its steady-state distribution is a Poisson distribution with rate λ/ν. As we will see in Proposition 3.3, the process describing the number of uncharged EVs in the system (i.e., Z(·)) behaves as a modified Erlang-A queue. The transition rates are shown in Figure 4. Figure 4: Transition rates of the process Z(·) (Erlang-A).

System Dynamics
In this section, we introduce the dynamics that describe the evolution of the system. We avoid a rigorous sample-path construction of the stochastic processes and we refer to [10,30] for background. For a constant r, let N r (·) be a Poisson process with rate r. The total number of EVs in the system at time t ≥ 0, Q(t), is given by where N λ (·), N ν,1 (·), and N ν,2 (·) are independent Poisson processes. Here, the number of EVs that arrive at the charging station during the time interval [0, t], is given by the process N λ and hence (2.1) describes the population in a well-known Erlang loss queue [27]. Another important process is the number of uncharged EVs in the system, Z(t), which can be written in the following form: where N µ t 0 L(Z(s))ds is the number of EVs that get fully charged during [0, t] and is independent of the aforementioned Poisson processes.
Last, the process which describes the number of fully charged EVs is given by Observe that in case K = ∞, (2.3) is reduced to the Erlang-A system; [20,47]. All the previous equations hold almost surely and are defined on the same probability space. It is clear that the vector process (Q(·), Z(·)) constitutes a two-dimensional Markov process. In the sequel, we are interested in finding the joint stationary distribution of (Q(·), Z(·)) and in deriving the fraction of fully charged EVs. Although the computation of the exact joint distribution does not seem promising, we are able to obtain exact bounds for the fraction of fully charged EVs in the next section.

Explicit bounds
The goal of this section is to give explicit results on some performance measures. In an EV charging setting, one may be interested in finding the fraction of EVs that get fully charged. This is an important performance measure from the point of view of both drivers and of the manager of the charging station. Using arguments from queueing theory, it can be shown that the fraction of EVs that get fully charged . Note that P s gives the probability that a vehicle leaves the charging station with fully charged battery. Thus, it is clear that the computation of P s requires the computation of the (joint) distribution of the process (Q(·), Z(·)). For the general model (i.e., for any K < ∞ and M < ∞) given in Section 2, define the steadystate probabilities p(q, z) := lim t→∞ P Q K M (t) = q, Z K M (t) = z . For simplicity, we use p(q, z) instead of p K M (q, z). These steady-state probabilities are characterized by the following balance equations: for (q, z) ∈ {R 2 + : z ≤ q}, we have that +µL(z + 1)1 {q =z} p(q, z + 1) + (q − z + 1)ν1 {q<K} p(q + 1, z). A closed form solution of the balance equations for any K and any M does not seem possible. However, we are able to obtain explicit solutions in some special cases. Below we derive bounds for P s based on three different cases: i) there is enough power for everyone (M = ∞), ii) there are enough parking spaces for everyone (K = ∞), iii) a full parking lot (Q(t) ≡ K). In the next proposition, we give upper and lower bounds for the fraction of EVs that get fully charged.

A full parking lot
Last, we consider the case where the parking lot is always full, i.e., the total number of EVs (uncharged and charged) is equal to the number of parking spaces. Roughly speaking we assume that the arrival rate is infinite, and that we replace (immediately) each departing EV by a newly arriving EV, which we assume to be uncharged. Hence, the total number of EVs always remains constant and it is equal to K. In other words, the original two-dimensional stochastic model reduces to a one-dimensional model. For this model, we find its steady-state distribution below. This result yields an upper bound for the number of uncharged EVs in the original system and hence a lower bound for the fraction of EVs that get fully charged. As we shall see later, the result in this section plays a crucial role in the study of the diffusion limit in the overloaded regime. Also, in the numerics, we see that a modification of the full parking lot case gives a very good approximation for the fraction of fully charged EVs.
On these assumptions, all newly arriving EVs are uncharged and so it turns out that the process describing the number of uncharged EVs in the system, {Z f (t), t ≥ 0}, is a birth-death process. In particular, the birth rate is ν(K − Z f (t)) and the death rate is equal to µ(Z f (t) ∧ M ). The state-state distribution of the aforementioned birth-death process is given in the following proposition.
Proposition 3.4. The steady-state distribution of the Markov process {Z f (t), t ≥ 0} is given by In Section 5, we validate these bounds in the three regimes; moderately, critically, and over-loaded. As we will see, the bounds are not very close in general. For this reason, we move to asymptotic approximations.

Asymptotic approximations
In this section, we present asymptotic approximations. First, we focus on the fluid approximation and then we move to three diffusion approximations. 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 in the sequence of systems, a superscript n will be appended to the system parameters and processes. In the remainder of this section, we assume that E [Q n (0)] and E [Z n (0)] are finite. Last, the proofs of the limit theorems are based on martingale arguments and are given in Sections 6.3-6.5. We give a rigorous proof for Theorem 4.5 in Section 6.3 and we omit the full details for the other proofs.

Fluid approximation
Here, we study a fluid model, which is a deterministic model that can be thought of as a formal law of large numbers approximation under appropriate scaling. We develop a fluid approximation for finite K, following a similar approach as in [21]. The main differences here are the finitely many servers in the system and that the state space consists of two regions: To obtain a non-trivial fluid limit, we assume that the capacity of power in the n th system is given by nM , the arrival rate by nλ, the number of parking spaces by nK, and we do not scale the time. The fluid scaling of the process describing the number of uncharged EVs in the charging station is given by n . This scaling gives rise to the following definition of a fluid model.
Note that (4.1) can be written as Further, the operator R(·) is Lipschitz-continuous in R + , which guarantees that (4.1) has a unique solution. In the proof of Proposition 4.2 below, we shall see that if the initial state of the fluid model solution z(0) ∈ [0, K], then z(t) ≤ K for any t ≥ 0. The last statement ensures that the definition of our fluid model is well-defined.
Next, we see that the fluid model solution can arise as a limit of the fluid scaled process Z n (·) n . The proof of the following proposition is based on martingale arguments and is given in Section 6.2.
Moreover, the next proposition states that the fluid model solution converges to the unique invariant point as time goes to infinity. Proposition 4.2. Let B and D be exponential random variables with rates µ and ν. We have that for any z(0) ∈ [0, K], z(t) → z * exponentially fast as t → ∞. In addition, z * is given by the unique positive solution to the following fixed-point equation In the proof of Proposition 4.2, we shall see that if z(0) = z * then z(t) = z * for any t ≥ 0, i.e., z * is the unique invariant point of (4.1). The point z * can be view as an approximation of the expected number of uncharged EVs in the system for the original (stochastic) model. Observing that E min{D, B max{1, z * M }} is the actual sojourn time of an EV in the system and that the quantity (λ ∧ νK) plays the role of the arrival rate, (4.2) can be seen as a version of Little's law. Further, if we allow a processor sharing discipline and infinity many servers (i.e., L(·) ≡ 1 and K = ∞), then (4.2) reduces to [21,Equation 4.1].
Remark 4.1. We shall see in the proof of Proposition 4.2 that the invariant point z * has a simpler form than (4.2) but the latter holds much more generally. If the random variables B and D are generally distributed and possibly dependent with E [B ∧ D] < ∞, then (4.2) still holds. The mathematical analysis then requires the use of measure-valued processes, which is beyond of the scope of the current work; for a heuristic approach see [4]. Thus, we present the proofs only under Markovian assumptions.
To ensure that z * is indeed a fluid approximation, we show that we can interchange the fluid and the steady-state limits. First, note that Z(·) has a limiting distribution. To see this, observe that Z(·) is bounded almost surely from above by the queue length of an Erlang A queue with M servers and infinite buffer. Alternatively, we can bound it by the queue length of an M/G/∞ queue. Now, using the same arguments as in [32], we conclude that Z(·) is a regenerative process and that there exists a stationary limit, Z(∞). The next proposition says that the stationary scaled sequence of random variables converges to the unique invariant point z * . Note that the arrival rate in an Erlang loss queue is known and it is equal to λ(1 − B(λ/ν, K)), where B(λ/ν, K) is the blocking probability in a loss system with K servers and traffic intensity λ/ν. Furthermore, λ(1 − B(λ/ν, K)) is asymptotically exact for our fluid approximation in the sense that λ(1 − B(nλ/ν, nK)) → (λ ∧ νK), as n → ∞. To improve the fluid approximation, we replace (λ ∧ νK) by λ(1 − B(λ/ν, K)), leading to Heuristically, we assume that an EV sees the system in stationarity throughout its sojourn and we use Little's law and a version of the snapshot principle; [31].
Having found the fluid approximation for the number of uncharged EVs in the charging station, we derive the fluid approximation for the fraction of EVs that get successfully charged. Let P s denote the probability that an EV leaves the parking lot with fully charged battery in the fluid model. It is given by P s = P (D > B max{1, z * /M }), where z * is the unique solution of (4.3). Under our assumptions, the explicit expression for this probability can be found. That is, Note that in region {z * ≤ M } the fraction of fully charged EVs is nothing else that the probability of minimum of two exponential random variables. We now focus on the fluid approximation for the number of uncharged EVs when the parking lot is full; Section 3.3. Analogous to Definition 4.1, we can define a fluid model, which we call z f (·). Proposition 4.4. Assume that the scaled parking spaces and the scaled power capacity are given by K n = Kn and M n = M n, respectively. If , and z f (t) → z * f as n and t go to infinity. Further, the limits can be interchanged and z * f is given by the following formula We give a heuristic approach how we can derive (4.4), skipping the proof, which can be done by using a similar procedure as in the general case. The intuition behind (4.4) is as follows. The actual sojourn time (in steady-state) of an EV in the system is given by By Little's law, we have that Dividing the last equation by n, yields Now, taking the limit as n goes to infinity, leads to Finally, z * f is given by the following fixed-point equation and solving the last equation leads to (4.4).
We shall see in the numerical examples in Section 5 that the fluid approximation is a good approximation of the fraction of fully charged EVs in most cases. However, especially in the underloaded regime and for small number of EV chargers, the error becomes larger. In the next section, we move to diffusion approximations.

Diffusion Approximations
In this section, we show diffusion limit theorems and diffusion approximations for the process describing the number of uncharged and the total number of EVs in the parking lot (vector process). To do it, we follow the strategy set up in [30] using the martingale representation.
First, we work in the Halfin-Whitt regime; Section 4.2.1. Using the "square-root staffing rule" to scale the system parameters, we extend [30, Theorem 7.1] and we obtain a limit which is a reflected twodimensional OU process with piecewise linear drift. Then, we derive an equation which characterizes its steady-state distribution, the so-called Basic Adjoint Relation (BAR). However, it turns out that the computation of the steady-state distribution is a hard problem, which is beyond the scope of this paper and it remains an open problem. To overcome this difficulty we consider more tractable asymptotic regimes.
The second asymptotic regime we consider here is an overloaded regime; Section 4.2.2. Assuming that the process describing the total number of EVs is in an overloaded regime and using the "square-root staffing rule" to scale the total power capacity in the system, we can show that the scaled vector process converges weakly to a one-dimensional limit. Thus, we can compute its steady-state distribution.
Last, in Section 4.2.3, we focus on the case that the parking times of the EVs are sufficiently large. We give a heavy traffic limit and a two-dimensional approximation for the vector process.

Diffusion approximation in Halfin-Whitt regime
The main goal in this section is to prove a two-dimensional diffusion limit for the vector process. For −∞ < β, κ < ∞, consider the following scaling: Define a sequence of diffusion-scaled processesQ n (·) := The allocation function in the n th system is given by L n (Z n (·)) := Z n (·) ∧ M n . We can then prove the following theorem.
,Q(·)). The limit satisfies the following two-dimensional stochastic differential equation Further, WẐ(·) and WQ(·) are driftless, univariate Brownian motions such that 2(ν + µ)E WẐ(t)WQ(t) = (2ν + µ)t. In addition,Ŷ (·) is the unique nondecreasing nonnegative process such that (4.5) holds and Adapting [30, Section 7.3], we can show that the last theorem also holds if we allow the arrival process to be a general stochastic process under the assumption that it is satisfies the functional central limit theorem.
The proof of Theorem 4.5 is given in Section 6.3.1 and is organized as follows.
1. We first establish a continuity result and show the existence and uniqueness of the candidate limit. (Proposition 6.1.) 2. We then rewrite the system dynamics using appropriate martingales and filtrations; see equations (6.16), (6.17), and Proposition 6.2.
3. Next, we show in Proposition 6.3 that the corresponding fluid-scaled processes converge weakly to deterministic functions.
4. Last, the proof of Theorem 4.5 is done by applying the martingale central limit theorem in [19] and Proposition 6.1.
Next, we focus on characterizing the joint steady-state distribution of the limit given by (4.5). Our approach is to find a functional equation which describes the joint steady-state distribution, the so-called Basic Adjoint Relation. The next step to use the BAR in order to obtain a key relation for the moment generating function of the vector proces. The piecewise linear drift and the existence of the reflection in (4.5) makes the key relation complicated and its analysis is beyond of the scope of this paper.
For any t ≥ 0, we know thatẐ(t) ∈ R andQ(t) ≤ κ. It is more convenient to transform the previous processes such thatẐ(t) ∈ R and κ −Q(t) ≥ 0. To do so, we recall that b 2 (x) = −νx. Thus, the diffusion limit can be written in the following integral form -see (4.5): whereŶ (·) is defined in Theorem 4.5. Multiplying by (−1), adding and subtracting the terms κ, νκt in the last equation, we obtain The processQ κ (t) represents the number of available spots in the parking lot at time t ≥ 0 (after scaling and after taking the limit as n goes to infinity). Furthermore,Ŷ (t) increases if and only ifQ κ (t) = 0. Define X(·) := (Ẑ(·),Q κ (·)) and note that each component of The boundary and the closure of G are given by ∂G = {x ∈ R 2 : x 2 = 0} andḠ = G ∂G, respectively. Now, observe that X(·) ∈Ḡ for any t ≥ 0. A geometrical representation of the space G and its boundary is shown in the next figure. Figure 5: The spaces G and its boundary for β > 0.
Before we continue the analysis of deriving the BAR, we note some properties for the processŶ (·), which is known as regulator. It is known that (Q(·),Ŷ (·)) satisfies a one-dimensional reflection mapping (or one-dimensional Skorokhod problem). The regulatorŶ (·) is continuous, nondecreasing and has the property   For more details we refer to [7, Lemma 3.1] and [26].
Equation (4.9) is rather complicated due to the piecewise linear term and the existence of the boundary measure. Although the analysis of (4.9) is beyond of the scope of the current paper, we conjecture that the Wiener-Hopf method [12] and boundary value techniques [13] may be applied.
It turns out that (4.9) remains quite complicated even if we assume K = ∞, i.e., no boundary measure. Contrary to the one-dimensional case [8], the steady-state distribution of (Z(·), Q(·)) cannot be written as a linear combination of two distributions. To see this, define π − (x) to be a bivariate Normal distribution with mean vector µ − = (0, 0) and covariance matrix Σ − = 1 1 1 ν+µ ν . In addition, let π + (x) be a bivariate Normal distribution with mean vector µ + = (−µβ/ν, 0) and covariance matrix where Φ(·) represents the cumulative probability function of the standard Normal distribution. It can be easily verified that π ∞ is indeed a probability distribution but it does not satisfy the correct marginal distribution ofQ(∞), as it is shown in Figure 6. It is well-known thatQ(∞) follows a Normal distribution with zero mean and variance ν+µ ν . For a discussion on this topic see [14]. In the sequel, we move in different asymptotic regimes in order to overcome this difficulty.

Diffusion approximation in an overloaded parking lot
In this section, we study an overloaded parking lot. First, we show a diffusion limit in the case that the parking lot is always full (see Section 3.3). We then show that the diffusion scaled vector process for the original system collapses to the one-dimensional limit in this case. Our motivation in this section comes from results in [1]. Specifically, the authors there show that under an appropriate scaling (including parameter and time scaling), the number of empty spaces in an overloaded parking lot behaves like an M/M/1 queue. However, here we need a modification of this result by dropping the time scaling. First, we define the dynamical equation that describes the evolution of the process of the number of uncharged EVs when the parking lot is always full. Let N f ν (·) and N f µ (·) denote two independent Poisson processed with rates ν and µ, respectively. For any t ≥ 0, we have that where Z f (0) ≤ K almost surely. Next, introduce our asymptotic regime. Take K n and M n such that K n = nK and M n = ν ν+µ K n + √ nβ, where K, β ≥ 0. The following proposition gives a diffusion limit for the scaled process describing the number of uncharged EVs, i.e., Z f (·).
, where the limit satisfies the following stochastic differential equation Moreover, the drift function is given by and W (·) is a standard Brownian motion.
Next, we give the steady-state distribution of the processẐ f (·). This can be done following [8,Equation 3]. Define the following truncated Normal probability density functions , for x ≤ β, and π + , for x > β, where σ 2 1 = 1 ν+µ νµK ν+µ and σ 2 2 = µK ν+µ . Now, the pdf ofẐ f (·) is given by for x ∈ R. Moreover, the constants are d 1 = 1 1+r and . Having studied the system when it is always full, we now move to the original stochastic model. The first step is to find a relation between the process that gives the number of uncharged EVs and the process that gives the empty parking spaces. Recall that the total number of EVs in the system is given by the following equation and the number of uncharged EVs is Define the stochastic process that describes the number of empty parking spaces in the system, E(t) := K − Q(t), t ≥ 0. Using the definition of the process E(·), we have that the system dynamics can be rewritten as follows (4.13) By (4.12), it follows that Applying the last equation in (4.13), yields (4.14) The last relation and an asymptotic bound for the process E n (·) (see Proposition 6.4) are the core elements we use to prove the main result in this section.
A diffusion approximation for the expected number of the original system in an overloaded regime is now given by and by using (4.11), we have that The asymptotic regime for an overloaded system leads to a one-dimensional approximation. In the next section, motivated by [41], we consider an asymptotic regime where we scale the parking times, which leads to a two-dimensional diffusion approximation.

Diffusion approximation for small parking rates
In this section, we study a diffusion approximation in case the parking rate ν is "small". First, we focus on the system with infinitely many parking spaces and we show a heavy traffic limit theorem; see Section 2.2.1 for an alternative model description when K = ∞. In this case, the limit is a twodimensional OU process with reflection. Then, making an overloaded assumption (for the uncharged EVs), we derive a two-dimensional OU limit process and we obtain the same limit if we assume a sufficiently large number of parking spaces.
Assume that K = ∞. Define the traffic intensity for this model as ρ := λ µM . Let µ, M be fixed. Further, define ν n = 1 n and λ n = µM (1 − c √ n ) for some constant c. Note, that 1−ρ n √ n → c as n → ∞, which is our heavy traffic assumption. Moreover, define the diffusion scaled process as follows The next proposition states a heavy traffic result for the two-dimensional scaled process.
Observe that the limit process in the last proposition depends on the reflection at zero. To overcome this difficulty, we consider an overloaded regime for the number of uncharged EVs. In this regime, the fraction of time that the processZ(·) spends at state zero is negligible; [41]. To this end, let λ, µ, M be fixed with λ > µM and ν n = 1/n. Modifying slightly the scaled processes, i.e.,  ,Q o (·)). The diffusion limit satisfies the following two-dimensional stochastic differential equation where W (·) = (W 1 (·), W 2 (·), W 3 (·), W 4 (·)) T , with W i (·) independent standard Brownian motions.
Note that we derive the same limit if we assume that K < ∞ and we scale the number of parking spaces in the n th system K n such that K n −λn √ n → ∞, as n → ∞. In this case, the fraction of time that the scaled processQ n o (·) spends on the boundary is negligible. This is made rigorous in the following lemma.  . Then Lemma 4.10 follows and holds true for any deviating sequence R n instead of K n −λn √ n . Last, note that we only need the weak convergence of the processQ n o (·) and the fact that the quantity K n −λn √ n goes to infinity. For the last convergence, it is enough to choose K n > λn.
The joint steady-state distribution of (Z o (·),Q o (·)), say π o (·, ·), is given by a bivariate Normal distribution with mean µ = (0, 0) and covariance matrix Σ = λ ν 2λ−µM 2ν 2λ−µM 2ν λ ν . Note, that Σ is indeed a variance matrix as it is positive definite. To see this, observe that where the first inequality holds by the assumption λ > µM . Now, for parameters of the original system such that λ > µM , sufficiently "small" ν, and K > λ/ν, we suggest the following diffusion approximation

Numerical evaluation
In this section, we validate numerically the previous bounds and approximations for three cases: the moderately (λ < νK), critically (λ = νK), and overloaded (λ > νK) systems. We focus on the expected number of uncharged EVs in the system and the probability that an EV leaves the charging station with fully charged battery (success probability). In First, we evaluate the fluid approximation. Table 1 gives the relative error between the expected number of uncharged EVs for the original system and the fluid approximation given in (4.2) for different values of the number of parking spaces K. For a given K, we give only the maximum relative error for 0 < M ≤ K. As expected, the relative error decreases as λ and K increase. In table 2, we present the relative error between the expected number of uncharged EVs for the original system and the modified fluid approximation given in equation (4.3). Not surprisingly, the relative error is much smaller in this case, as we can see in Table 2. For high values of λ, K the relative error is approximately 2% rather than 10-20%. In addition, the modified fluid approximation seems to be reasonable also in the moderate regime. Next, we evaluate the approximation in case of a full parking lot (see Section (3.3)) and the diffusion approximation in an overloaded regime given by (4.15). To improve the approximations, we directly modify them by replacing the parameter K by the expected number of the total EVs in the original system, i.e., λ(1 − B(λ/ν, K)). Table 3 gives the relative error for E [Z f (∞)]. As we expect, it decreases as λ and K increase. Furthermore, this approximation results to small relative errors in all regimes   (< 3%). The (prelimit) approximation E [Z f (∞)] is better than the modified diffusion approximation in (4.15), as we see in Table 4.
In the sequel, we depict the bounds in (3.2), the modified bound in (3.3) -where we replace K by λ(1 − B(λ/ν, K)) -(dotted line), the modified fluid approximation (4.4) (dashed line), and the modified diffusion approximation in (4.15) (dash-dot line). In Figures 7-18, the vertical axes give the probability that an EV leaves the parking lot with fully charged battery (success probability) and the horizontal axes give the ratio M/K. For each regime, we plot the success probability for K = 10, 20, 30, 50. In the moderated regime, the lower bound (K = ∞) is very close for high values of parking spaces. This is not surprising because the time that the process spends on the boundary is negligible in this case. The fluid approximation seems to be quite good in most of the cases and the diffusion approximation does not improve the fluid one. Last, note that the modified bound (3.3) does not give a lower bound of the original system. However, it seems to be the best approximation for all the cases; even in the moderated regime.    Note that the probability (in stationarity) that an EV leaves with a full battery can also be given by Now, using the assumptions that Z K M (∞) ≤ K and M ≤ K, the following inequality holds almost surely By the last inequality, we have that P D > B max 1, which proves the upper bound. Note that the upper bound is nothing else but the minimum of two exponential random variables.
We move now to the proof of the lower bound in (3.2). We note that it is equivalent to the following inequality by using coupling arguments. Then, (6.1) follows. Fix a sample path ω ∈ Ω. Assume that Z ∞ M (0) = Z K M (0) and take identical arrival, charging, and parking times for both systems. Define T * = inf{t > 0 : Q K M (t) = K}. It follows that Z ∞ M (t) = Z K M (t), for t ≤ T * . After time T * , the blocked arrivals in the loss queue will enter in the queue with infinite many parking spaces. That is, Z ∞ M (t) ≥ Z K M (t) for all t ≥ 0. Removing now the conditioning on the sample path ω, we derive Z ∞ M (t) ≥ st Z K M (t) for all t ≥ 0, and by the existence of the stationary distribution we have that . It remains to show (3.3). Let (Q λ (·), Z λ (·)) denote the total number of EVs and the number of uncharged EVs if the arrival rate is λ. First, using coupling arguments, we prove that if λ 1 ≤ λ 2 then Q λ1 (t) ≤ st Q λ2 (t) and Z λ1 (t) ≤ st Z λ2 (t) for any t ≥ 0. Assume the following coupling: if an arrival occurs to the system with arrival rate λ 1 , it also occurs to system with arrival rate λ 2 . Hence, as λ 1 ≤ λ 2 , there are more arrivals in the second system. Further, we assume that all other parameters, i.e. µ, ν, M , K, are equal in both systems. Assume that both systems start empty and define T * * = inf{t > 0 : Q λ2 (t) = K}. As in the second system there more arrivals, we have that Q λ1 (t) ≤ Q λ2 (t) and Z λ1 (t) ≤ Z λ2 (t), for t ≤ T * * . By the Markovian assumptions, we have that the residual charging and parking times are exponential with rare µ and ν. That is, at any new event after time T * * , we can resample the charging and parking times and hence the probability of a departure in the system with arrival rate λ 1 is higher or equal to the probability of a departure in the system with arrival rate λ 2 . In other words, for t ≥ 0 and for In the sequel, we see that Z f (·) can be arise as the limit of Z λ (·) as λ → ∞, assuming that Q λ (0) . To see this, first observe that Q λ (·) d → K as λ → ∞. Now, combining (2.1) and (2.3), we have that Taking λ → ∞ and using the continuous mapping theorem, we have that and where the last equality follows by (4.10). Furthermore, Z λ (·) is non-decreasing. That is, Z λ (t) ≤ st Z f (t) for any t ≥ 0 and by the existence of the stationary distributions we obtain Proof of Proposition 3.2. Note that the distribution p Q (q) corresponds to the stationary distribution of a one-dimensional Erlang loss system. Furthermore, by [27, Section 1.3] we know that Thus, the probability of an empty system is p e (0, 0) = p Q (0). As it is well known that a solution of the balance equations of a Markon process is unique, we shall show that p e (q, z), for z ≤ q satisfies the flow balance equations (3.1). Then, the proof of the proposition is completed. First, we note the relations between p e (q + a, z + b) and p e (q, z) for a, b ∈ {−1, 0, 1}. By (3.5), we obtain that Now, applying the previous relation in (3.4) we have that Working analogously, we derive the following relations p e (q − 1, z − 1) =z ν + µ λ p e (q, z), (6.3) p e (q + 1, z + 1) = 1 z + 1 λ ν + µ p e (q, z), (6.4) p e (q, z + 1) = q − z z + 1 ν µ p e (q, z), (6.5) z). (6.6) Using the above equations and recalling that L(z) = z when M = K, the right-hand side of (3.1) for 0 < z, q < K and z = q can be written as follows That is, p e (q, z) satisfies (3.1) for 0 < z, q < K and z = q. To show that p e (q, z) satisfies (3.1) for 0 < q < K and z = 0, we apply (3.4) and (6.2) in the right-hand side of (3.1). This leads to (q + 1)νp e (q + 1) In the same way, we show that the right-hand side of (3.1) for 0 < z < K and q = K becomes λK ν λ The last quantity is equal to Using again the relations (6.3), (6.4) and (6.6), the right-hand side of (3.1) is written for q < K and q = z as follows Using again relations (6.3)-(6.6), it follows immediately that p e (q, z) satisfies Equation (3.1) also for the remaining cases, i.e., (q, z) = (0, 0), (q, z) = (K, K), and last (q, z) = (K, 0).
Proof of Proposition 3.3. We show that Z ∞ M (·) behaves as modified Erlang-A queue. Although we adapt the proof in [47, Section 6.6.1], we briefly describe it here for completeness. First, we write the flow balance equations for the Markov process Z ∞ M (·) and then we solve them. The balance equations for the Markov process Z ∞ M (·) are given by and for z = 0 we have that λp Z (0) = (ν + µ)p Z (1).
Using the last equation and (6.7), we derive inductively the following relations The balance equations now can be simplified as follows Note that the balance equations can be simplified to Applying z = M − 1 in the first equation of (6.10), we obtain and recursively we have that Change the variable z = M − i in the last equation yields Working analogously, by the second equation of (6.10) we derive which leads to Last, π f (M ) is determined by the normalization equation K z=0 π f (z) = 1.

Proofs for Section 4.1
Proof of Proposition 4.1. In this proof, we use martingales arguments. Define the following filtration F n t := σ (Z n (0), Q n (0), Z n (s), Q n (s) : 0 ≤ s ≤ t) augmented by including all the null sets for t ≥ 0 and n ≥ 1. Applying the fluid scaling to the dynamical equation (2.3), we have that Defining the operatorM n r = 1 n (N r (·) − r·) and following [30], we can write The term 1 n t 0 1 { Q n (s) n =K} dN (λns) denotes the number of EVs that are lost due finding the system full under the fluid scaling. By the discussion in [33, Section 6.7], [24,Equation 3.44], and by the assumption Further, observing thatM n r (·) are zero mean martingales with respect to the filtration F n t for any n ∈ N (cf. [30]) and taking n → ∞, we derive that Z n (·) n d → z(·). Moreover, the limit function is characterized by the following functional equation which is equivalent to (4.1).
Proof of Proposition 4.2. It is not hard to solve the ODE (4.1) explicitly in the two regions, namely So, for any t ≥ 0 it follows z(t) ≤ K and hence Definition 4.1 is well defined.
In the sequel, we show that z(t) converges as t goes to infinity to the following point Then, we show that z is unique and using the Markovian assumptions, we see that it is equivalent to (4.2) and hence z = z * . Also, observe that z * is an invariant point. Indeed, if we set z(0) = z * , then by (6.11) we have that z(t) = z * for t ≥ 0. First, assume that z 1 ≤ M . If z(0) ≤ M , then z(t) ≤ M for any t ≥ 0. To see this, note that if z(0) − z 1 ≤ 0, then we have that On the other hand, if the quantity z(0) − z 1 is positive, we show that there does not exist t * > 0 such that z(t * ) > M . Suppose that there exists t * such that z(t * ) > M , we have that by assumption z(0) − z 1 > 0, the last inequality leads to Now, by the assumption z(0) ≤ M , we obtain M −z1 z(0)−z1 ≥ 1 and hence t * ≤ 0 which yields a contradiction. Next, we assume that z(0) > M . In this case, we show that there exists t * 2 such that z(t * 2 ) ≤ M . First, observe that if z 1 ≤ M then z 2 ≤ M . Hence, z(0) − z 2 > M − z 2 ≥ 0. Now, we note that Let t * 2 be the first time that z(t) ≤ M starting from a point z(0) > M . Setting now as initial point z(t * 2 ), we conclude that if z 1 ≤ M then z(t) ≤ M for t ≥ t * 2 and hence z(t) = z 1 + (z(t * 2 ) − z 1 ) e −(ν+µ)t , for t ≥ t * 2 , which yields This concludes the proof for the case z 1 ≤ M . The case z 2 > M follows the same logic. Now, we prove that (4.2) is equal to z. To this end, observe that where B is an exponential random variable with E [B ] = 1 µ max 1, z * M . We note that (4.2) can be written as Solving the last equation yields z * = z.
To conclude the proof of Proposition 4.2, it remains to show the uniqueness of the invariant point z * . In other worlds, we show that (6.12) (and hence (4.2)) has a unique solution. It is not hard to see that if z 1 < M then z 2 < M . So, z 2 cannot be the solution of (6.12). That is, z 1 is the unique solution of (6.12). On the other hand, if z 1 > M (i.e., it is not solution of (6.12)), then z 2 > M . That is, z 2 is the unique solution of (6.12). Last, if z 1 = M , then we have that z 2 = z 1 = M . In any case, (6.12) has a unique solution.
Proof of Proposition 4.3. Let Z n (∞) n be the stationary fluid scaled number of uncharged EVs. We know that 0 ≤ Z n (∞) ≤ K n , which yields Z n (∞) n ≤ K almost surely. In other worlds, the sequence of random variables Z n (∞) n is stochastically bounded in R and hence it is tight. Now, we consider the process {Z n (t), t ≥ 0} starting at point Z n (∞). That is, Z n (t) We start the analysis by establishing a continuity result, which can be proved by using results in [30] Proposition 6.1. Let t ≥ 0 and −∞ < κ < ∞. Consider the following system 13) where b i are positive constants, h i : R → R satisfy h i (0) = 0 and are Lipschitz continuous functions for i = 1, 2, and x 2 (t) ≤ κ. In addition, y(·) is a nondecreasing nonnegative function in D[0, ∞) such that (6.13) holds and ∞ 0 1 {x2(t)<κ} dy(t) = 0. Given b i ∈ R and g i (·) ∈ D[0, ∞), we have that the system (6.13) has a unique solution (x 1 (·), x 2 (·), y(·)). Moreover, the functions (x 1 (·), x 2 (·), y(·)) are continuous in D[0, ∞) 3 if D[0, ∞) is endowed with the uniform topology over bounded intervals or the J 1 topology.
In order to continue our analysis, we need to define appropriate filtrations. Take the following filtrations, for n ≥ 1, In the sequel, we work with the filtrations F n t = σ F n t,1 , F n t,2 , augmented by including all the null sets for t ≥ 0 and n ≥ 1. Now, notice that the system dynamics (2.1) and (2.3) can be rewritten in the following form The process Y n (t) counts all the customers that are lost when all the servers (chargers) are busy up to time t in the n th system. Defining the operator M r (·) := N r (·)−(r·), where "r" indicates the rate of the Poisson process N r (·), the system dynamics take the following form: In order to derive appropriate equations (in the pre-limit) for the diffusion scaled processes, subtract and add the terms n ν+µ ν , n ν+µ ν t in (6.14) and the terms n, nµt, nνt in (6.15), and then divide both by √ n. Recalling that L n (Z n (t)) = Z n (t) ∧ M n , λ n = n(ν + µ), M n = λ n ν+µ + β √ n and observing that M n − n = β √ n and λ n t − n(ν + µ)t = 0, we obtain the following equations for the diffusion scaled processesQ n (·) andẐ n (·), shows that the processesM r (·) are martingales.
Proof. Fix n ≥ 1. The result for the processM n λ n (·) follows immediately by applying [30, Lemma 3.1]. Now, note that by the system dynamics (2.1), (2.3), the fact that Q n (t) = Z n (t) + C n (t), and (2.2) we have thatM Since N ν,1 (·), N ν,2 (·) are independent Poisson processes, [30, Lemma 3.1] implies thatM n ν,i (·) are F nmartingales for i = 1, 2. Observe thatM n ν,2 (·) is adapted to the filtration F n as the latter contains all the information about the processes Q n (·) and Z n (·) for fixed n. This is enough to determine the process C n (·) at any t ≥ 0. Using (6.22), the assumptions E [Z n (0)] , E [Q n (0)] < ∞, the inequalities Z n (t) ≤ Z n (0) + N λ n (t), (6.23) and adapting the proof in [30,Lemma 3.4], we obtain that for fixed n ≥ 1, the following moments conditions are satisfied Also, again by [30,Lemma 3.4], we derive that the following moments related to the Poisson processes are finite.  is an F n -martingale. Last, note that by (6.22), asM n µ (·) is adapted to the filtration F n , we obtain that M n ν t 0 Q n (s)ds is also an F n -martingale.
In order to apply the martingale central limit theorem, we first need to show that the corresponding fluid scaled processes converge to a deterministic function (step 3), i.e., Proof. We prove the fluid limits using the martingale representations (6.16) and (6.17). If the sequences {Ẑ n (·)} and {Q n (·)} are stochastically bounded in D[0, ∞) then by [30,Lemma 5.9], we have that where the function η : [0, ∞) → [0, ∞) is defined by η(t) ≡ 0. Note that the last limits are equivalent to (6.24) and (6.25). The diffusion scaled processesẐ n (·) andQ n (·) have the martingale representation (6.16) and (6.17). In order to prove that they are stochastically bounded in D[0, ∞), it is enough to show that the corresponding martingales are stochastically bounded in D[0, ∞); see [30,Lemma 5.5]. But by [30,Lemma 5.8] the martingales are stochastically bounded if the sequence of their predictable quadratic variations (6.18)-(6.21) are stochastically bounded in R for each t ≥ 0. To prove that the sequences of the predictable quadratic variations of the martingales in expressions (6.16) and (6.17) are stochastically bounded in R, we use (6.23), [30, Lemma 6.2] and the fact that for any t ≥ 0, Q n (t) n ≤ K n n < ∞. The predictable quadratic variation for the arrival process (6.18) is obviously bounded as it is a deterministic function. For (6.18), we have that The result for (6.20) follows by applying (6.23) and [30,Lemma 6.2]. Last, applying the inequality Q n (t) n ≤ K n n in (6.21) we obtain Before we move to the final step of the proof of Theorem 4.5, we make a remark for the fluid limit of the number of fully charged EVs in the , which we need later. Remark 6.1. Note that the diffusion scaled process {Ĉ n (·)} is also stochastically bounded in D[0, ∞) and the fluid limit is given by the difference ν+µ ν − 1 = µ ν . Now, we are ready to put all the pieces together, leading to the last step of the proof of Theorem 4.5.
Recalling thatM n r (·) := N f r (·)−(r·) √ n and dividing the last equation by √ n, we have that Observing that the quantity t 0 K n −Z n f (s) n ds is stochastically bounded and allowing n → ∞ in the last equation, we derivê where W 1 (t) and W 2 (t) are (independent) standard Brownian motions. Last, by the properties of Brownian motion, we get where W (t) is a standard Brownian motion.

Proof of Theorem 4.7
The rest of this section gives a proof of Theorem 4.7. We first show a bound for the process E n (·).
Proposition 6.4. Let T > 0. We have that for any > 0 there exists n such that for any n ≥ n , where L is a positive constant.
Before we proceed to the proof of Proposition 6.4, we show a preliminary result. Choosing q n = L(nT ) 1/4 + L log(nT ) and applying Lemma 6.5, we have that for any > 0 there exists n such that P sup 0≤t≤nT E M (t) ≤ L(nT ) 1/4 + L log(nT ) > 1 − , for n > n . This concludes the proof as is arbitrary.
By Proposition 6.4 follows that P (G n ) → 1, an n → ∞. In the sequel, we assume that ω ∈ G n .
Let T > 0 and ω ∈ G n . We have that as n → ∞. The result now follows by adapting the proof of Proposition 4.6.

Proofs for Section 4.2.3
Proof of Proposition 4.8. First, we note that it is enough to show the result for M = 1. Then we replace µ by µM ; see [41,Remark 5].
Scale the time by n, and add and subtract the means of the Poisson processes in (2.1) and (2.3) to obtain Q n (nt) = Q n (0) + (N λ n (nt) − λ n nt) − N ν n n By the overloaded assumption, it follows that √ n t 0 1 {Z n (s)=0} ds → 0, as n → ∞; [41]. Now, taking the limit n → ∞, we derive

Concluding Remarks
This paper proposes to model an electric vehicle charging model by using layered queueing networks. We develop several bounds and approximations for the number of uncharged EVs in the system and the probability that an EV leaves the charging station with fully charged battery. In the numerical examples, it seems that a modification of the approximation for a full parking lot leads to a good approximation. Further, the fluid approximation seems to be good in most cases and we believe that the diffusion approximation in the Halfin-Whitt regime will improve the fluid approximation. Unfortunately, the exact (or even numerical) solution of (4.8) seems very hard. From an application standpoint, it is important to remove various model assumptions. If parking and charging times are given by the (possibly dependent) generally distributed random variables B and D, we can develop a measure-valued fluid model by extending [21]. In addition, we can include in the model time-varying arrival rates, multiple EV types, and multiple parking lots, thus extending results in [32]. Moreover, the distribution grid (low-voltage network) plays a crucial role and it should be included in the model; see [2] for a heuristic approach. For another application in EV-charging including the distribution grid, see [9] where simulation results are presented for a Markovian model.