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 a customer (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 [25]. 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 toward the current grid infrastructure. This is illustrated in [24]; 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 [12]. 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; for example, 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 [24].
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 within 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. Finally, we assume that all available power is shared at the same rate to all cars that need charging. The available power that can be 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 steadystate 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 [22]. 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. Finally, 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 martingale arguments [32].
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,48]. Both papers consider a gradient scheduler to control delays. Next, [47] presents a queueing model for battery swapping while [40] is an early paper on a queueing analysis of EV charging, focusing on designing safety control rules (in terms 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 [39] 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 [38]. In [31], 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 [37]. In [30], 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 networks have been successfully applied in analyzing interactive networks in communication 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 [36,45]) 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 manufacturing systems, for example, [17,18]. 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 communication networks, for example, Web-based multitiered system architectures. In such environments, different applications compete for access to shared infrastruc-ture resources, both at the software level and at the hardware level. For background, see [42,43].
The paper is organized as follows: In Sect. 2, we provide a detailed model description-in particular, we introduce our stochastic model and we define the system dynamics. Next, in Sect. 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 Sect. 5. Finally, all proofs are gathered in Sect. 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 > 0, let D[0, T ] 2 := D[0, T ]×D[0, T ] be the two-dimensional Skorokhod space, i.e., the space of two-dimensional realvalued functions on [0, T ] that are right continuous with left limits endowed with the J 1 topology; cf. [11]. Observe that as all candidate limit objects we consider are continuous, we only need to work with the uniform topology. It is well known that the space (D[0, T ] 2 , J 1 ) is a complete and separate metric space (i.e., a Polish metric space) [7]. 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, 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. Finally, 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 charging, 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 of 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 Fig. 1.
Finally, 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 transition rates in the interior and on the boundary are shown in Fig. 2.

Alternative model description in the case of infinitely many parking spaces
Here, we give an alternative description of our model in the case where 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 Fig. 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 an 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 Charged 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 in 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 Fig. 4.

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 [11,32] 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 λ In other words, and by the properties of the Poisson process, we have and hence (2.1) describes the population in a well-known Erlang loss queue [29]. 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.
Finally, the process which describes the number of fully charged EVs is given by Observe that in the case K = ∞, (2.3) is reduced to the Erlang-A system [21,49]. 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. Note that the fraction of EVs that get fully charged (in steady-state) equals P s = , since in equilibrium νE C K M (∞) is the departure rate of fully charged EVs and νE Q K M (∞) is the departure rate of all EVs. Further, P s gives the probability that a vehicle leaves the charging station with a fully charged battery. For the general model (i.e., for K < ∞ and M < ∞) given in Sect. 2, define the steady-state 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 (3.1) 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.

Proposition 3.1 Let C K M (∞) and Q K M (∞) be the number of fully charged EVs and the total number of EVs in steady state for any K and any M. We have that
.

(3.2)
Moreover, an additional lower bound is given by The proof of this proposition makes use of coupling arguments and stochastic ordering of random variables, and it is given in Sect. 6. We now briefly present the solution of the balance equations for the three special cases described above.

Enough power for everyone
Assume that K is finite and that there is enough power for all EVs to be charged at a maximum rate, i.e., M = K . In this case, the allocation function takes the form L(Z (t)) = Z (t), and the balance equations can be solved explicitly and are given below.
Moreover, the probability of an empty system is given by

Enough parking spaces for everyone
In the second case, we assume that there are infinitely many parking spaces; i.e., K = ∞ and M < ∞. In this case, all EVs can find a free position and the process Z (·) can be modeled as a Markov process itself where its transition rates are given in Fig. 4. We see in the next proposition that the process Z (·) behaves as a modified Erlang-A model with M servers [49]. The main difference here is that EVs can leave the system even if they are in service (i.e., are getting charged).

A full parking lot
Finally, 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. Under 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 steady-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
In Sect. 5, we validate these bounds in the three regimes: moderately, critically, and overloaded. 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 Sect. 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. Finally, the proofs of the limit theorems are based on martingale arguments and are given in Sects. 6.2-6.5. We give a rigorous proof for Theorem 4.5 in Sect. 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 [22]. The main differences here are the finitely many servers in the system and that the state space consists of two regions: To obtain a nontrivial 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 Z n (·) 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, we shall see that if the initial state of the fluid model solution is z(0) ∈ [0, K ], then z(t) ≤ K for any t ≥ 0. This last statement ensures that 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 Sect. 6.2.

Remark 4.1
We point out that we can extend the previous proposition to the case This requires a modification of Definition 4.1. In particular, we need to replace the quantity 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 viewed as an approximation of the expected number of uncharged EVs in the system for the original (stochastic) model. Observing that the quantity E min D, B max 1, z * M is the actual sojourn time of an uncharged 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 [22,Eq. (4.1)].

Remark 4.2
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 measurevalued 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 [34], 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 * .

Proposition 4.3 The stationary fluid-scaled sequence of random variables Z n (∞)
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 λ( 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 [33]. 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 a fully charged battery in the fluid model. This is given by P s = P (D > B max{1, z * /M}), where z * is the unique solution of (4.3). Under our assumptions, an explicit expression for this probability can be found. That is, We now focus on the fluid approximation for the number of uncharged EVs when the parking lot is full (see Sect. 3.3). Analogously 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 = K n and M n = Mn, 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 to deriving (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: Let S n (B) be the sojourn time (in steady state) of an uncharged EV in the n th system. By Little's law, we have that By the discussion after Theorem 2.4 in [23], because we observe the system in steady state at time 0, the number of uncharged EVs hardly changes for large enough n. By the Applying the last equation in (4.5) and dividing 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 this last equation leads to (4.4).
We shall see in the numerical examples in Sect. 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 (the vector process). To do this, we follow the strategy set up in [32] using the martingale representation.
First, we work on the Halfin-Whitt regime (see Sect. 4.2.1). Using the "square-root staffing rule" to scale the system parameters, we extend [32, Theorem 7.1] and we obtain a limit which is a reflected two-dimensional 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.
The second asymptotic regime we consider here is an overloaded regime (Sect. 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.
Finally, in Sect. 4.2.3, we focus on the case where 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 the 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: In addition,Ŷ (·) is the unique nondecreasing nonnegative process such that (4.6) holds and ∞ 0 1 {Q(t)<κ} dŶ (t) = 0. Adapting [32,Sect. 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 satisfies the functional central limit theorem.
The proof of Theorem 4.5 is given in Sect. 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 Eqs. (6.15), (6.16), and Proposition 6.2. 3. Next, we show in Proposition 6.3 that the corresponding fluid-scaled processes converge weakly to deterministic functions. 4. Finally, the proof of Theorem 4.5 is completed by applying the martingale central limit theorem in [20] and Proposition 6.1.
Next, we focus on characterizing the joint steady-state distribution of the limit given by (4.6). Our approach is to find a functional equation which describes the joint steady-state distribution, the so-called basic adjoint relation. The next step is to use the BAR in order to obtain a key relation for the moment-generating function of the vector process. The piecewise linear drift and the existence of the reflection in (4.6) make 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.6): whereŶ (·) is defined in Theorem 4.5. Multiplying by (−1), adding and subtracting the terms κ and νκ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 X(·) is a semimartingale. Let G = {x := (x 1 , x 2 ) ∈ R 2 : x 2 > 0}. 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 (Fig. 5). β G Before we continue the analysis of deriving the BAR, we note some properties for the processŶ (·), which is known as the 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 [8, Lemma 3.1] and [28].
In the sequel, we focus on deriving a functional equation which characterizes the steady-state distribution π(·, ·) of the process {X(t), t ≥ 0}, provided that it exists. To handle the boundary of the space G, we define a measure σ on (Ḡ, B(Ḡ)) given by where the boundary measure σ is defined in (4.8) and L is the second-order operator, i.e., The next step is to derive a key relation between the moment-generating functions of π and σ . Let us define the two-dimensional moment-generating function (MGF) of π , for θ := (θ 1 , θ 2 ) ∈ R 2 , and θ · x := θ 1 x 1 + θ 2 x 2 . In the same way, we define the one-dimensional MGF of σ , 1 σ (d x).
Equation (4.10) is rather complicated due to the piecewise linear term and the existence of the boundary measure. Although the analysis of (4.10) is beyond the scope of the current paper, we conjecture that the Wiener-Hopf method [13] and boundary value techniques [14] may be applied.
It turns out that (4.10) remains quite complicated even if we assume K = ∞, i.e., no boundary measure. Contrary to the one-dimensional case [9], 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 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 is shown in Fig. 6. It is well known thatQ(∞) follows a normal distribution with zero mean and variance ν+μ ν . For a discussion on this topic, see [15]. In the sequel, we move in different asymptotic regimes.

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 Sect. 3.3). We then show that the diffusion-scaled vector process for the original system collapses to the onedimensional 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 timescaling), 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 timescaling.
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 processes with rates ν and μ, respectively. For any t ≥ 0, we have that where Z f (0) ≤ K almost surely.
Next, we 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 (·).

Proposition 4.6 Define the scaled processẐ
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 by following [9,Eq. (3)]. Define the following truncated normal probability density functions: 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.14) By (4.13), it follows that Applying the last equation in (4.14) yields 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 in the original system in an overloaded regime is now given by and, by using (4.12), 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 the case where 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 Sect. 2.2.1 for an alternative model description when K = ∞. In this case, the limit is a two-dimensional OU process with reflection. Then, making an overloaded assumption (for the uncharged EVs), we derive a twodimensional 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 √ n(1 − ρ n ) = c, 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.

WZ (·) and WQ(·) are driftless, univariate Brownian motions such that E WZ (t)WQ(t) = t/2.
Observe that the limit process in the last proposition depends on the reflection at zero, which makes the calculation of the joint distribution hard. We next consider an overloaded regime for the number of uncharged EVs. In this regime, the processZ n (·) will not reach zero for large enough n [41]. To this end, let λ, μ, M be fixed with λ > μM and ν n = 1/n. Modifying slightly the scaled processes, i.e., we are able to show the following proposition. . The diffusion limit satisfies the following twodimensional 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.  (D[0, ∞), J 1 ) which is a complete and separate metric space [44, Corollary 3.1]. Then, Lemma 4.10 follows and holds true for any deviating sequence R n instead of K n −λn √ n . Finally, 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.

Remark 4.3 The sequence {Q n o (t), t ≥ 0} is stochastically bounded as it converges in distribution in
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 Σ = . Note that Σ is indeed a covariance 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 a fully charged battery (the success probability). In all the numerical examples, we solve the flow balance equation (3.1) using standard numerical methods and we let ν = 1. Finally, the relative error is calculated by the following formula: RE =

|E[Z (∞)]−E[Z ap (∞)]| E[Z (∞)]
100%, where E [Z (∞)] denotes the expected number of EVs in the original system by solving the two-dimensional Markov process and E Z ap (∞) denotes the expected number of uncharged EVs for the aforementioned approximations.
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 and for μ = 1. 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 Tables 2, 3, and 4, we present the relative error between the expected number of uncharged EVs for the original system and the modified fluid approximation given in Eq.   Table 2 Evaluation of the modified fluid approximation for μ = 1/2  Table 3 Evaluation of the modified fluid approximation for μ = 1 values of the charging rate: μ = 1/2, 1, 3/2. Not surprisingly, the relative error is much smaller in this case. For high values of λ and K the relative error is approximately 2-5% 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 the case of a full parking lot (see Sect. 3.3) and the diffusion approximation in an overloaded regime given by (4.16). To improve the approximations, we directly modify them by replacing the parameter K by the expected total number of EVs in the original system, i.e., λ (1 − B(λ/ν, K ))/ν. In particular, the modified E Z f (∞) is calculated based on the stationary distribution in Proposition 3.4 after replacing K by λ (1 − B(λ/ν, K ))/ν. Moreover, the modified diffusion approximation in an overloaded regime is given by where E Ẑ f (∞) is given in (4.17). Tables 5, 6, and 7 give the relative error for the modified E Z f (∞) and μ = 1/2, 1, 3/2, respectively. As we expect, it decreases as λ and K increase for all values of the charging rate μ. Furthermore, this approximation results in small relative errors in all regimes (< 3%). The (prelimit) approximation E Z f (∞) is better than the modified diffusion approximation in (4.16), as we see in Tables 8, 9, and 10.
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.3) (dashed line), and the modified diffusion approximation in (4.16) (dash-dot line). Note that the modified bound in (4.16) becomes  Table 7 Evaluation of the modified E Z f (∞) for μ = 3/2   Table 10 Evaluation of the modified diffusion approximation in an overloaded regime for μ = 3/2 where the modified E Z f (∞) is as we discussed earlier. In Figs. 7,8,9,10,11,12,13,14,15,16,17, and 18, the vertical axes give the probability that an EV leaves the parking lot with a fully charged battery (the 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 and μ = 1. In the moderate 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. Finally, note that the modified bound (3.3) does not give a lower bound for the original system. However, it seems to be the best approximation for all the cases, even in the moderate regime. Recall that the probability (in stationarity) that an EV leaves with a full battery is given by P s = .
As the parking times do not depend on M, we have that (as the total number of EVs remains the same for both systems). We now proceed to show the last inequality using coupling arguments. Fix a sample path ω ∈ Ω. Assume for simplicity that Z K M (0) = Z K K (0) ≤ M, and take identical arrival, charging, and parking times for both systems. Note that the invariant distributions do not depend on the initial conditions, and so what we assume on the initial conditions is without loss of generality. Further, define T := inf{t ≥ 0 : Z K M (t) = M}. Note that the arrival process is the same for both systems as the total number of EVs remains the same for both systems. That is, Z K M (t) = Z K K (t) for t ≤ T since the charging rate is the same for both systems. Moreover, using the fact that Z K M (t) ≤ K and M ≤ K , the following inequality holds for any t > T : . Hence, the lower bound in (3.2) follows. 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 are 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 rates μ 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 x > 0, P Q λ 2 (t) ≤ x ≤ P Q λ 1 (t) ≤ x and P Z λ 2 (t) ≤ x ≤ P Z λ 1 (t) ≤ x . The last relation is equivalent to Q λ 1 (t) ≤ st Q λ 2 (t) and Z λ 1 (t) ≤ st Z λ 2 (t), for t ≥ 0. In the sequel, we see that Z f (·) can arise as the limit of Z λ (·) as λ → ∞, assuming that Q λ (0) Taking λ → ∞ and using the continuous mapping theorem, we have that Z λ (·) and where the last equality follows by (4.11). Furthermore, Z λ (·) is nondecreasing. That is, Z λ (t) ≤ st Z f (t) for any t ≥ 0 and, by the existence of the stationary distributions, we obtain Z λ (∞) ≤ st Z f (∞). By the last inequality, it follows that and hence (3.3).

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 [29, Sect. 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 Markov process is unique, we shall show that p e (q, z) for z ≤ q satisfies the flow balance equation (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 Note that the balance equations can be simplified to Applying z = M − 1 in the first equation of (6.9), we obtain and recursively we have that Changing the variable z = M − i in the last equation yields Working analogously, by the second equation of (6.9) we derive which leads to Finally, π f (M) is determined by the normalization equation K z=0 π f (z) = 1.

Proof of Proposition 4.1
In this proof, we use martingale arguments. Define the 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 [32], 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 to finding the system full under the fluid scaling. By [26,Proposition 4.11], this is a weakly convergent sequence. Further, by [35,Propositions 6.17 and 6.21] and by the On the other hand, if λ/ν < K , we have that q(t) < K for t > 0. Now, by the proof of [26,Theorem 3.6] (see [26,Eqs. (3.34) and (3.40)]), it turns out that Further, observing thatM n r (·) are zero mean martingales with respect to the filtration F n t for any n ∈ N (cf. [32]) 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 if z(0) − z 2 ≥ 0. So, for any t ≥ 0, it follows that 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.10) 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. Supposing that there exists t * such that z(t * ) > M, we have that by assumption z(0) − z 1 > 0, so the last inequality leads to Now, by the assumption z(0) ≤ M, we obtain M−z 1 z(0)−z 1 ≥ 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 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 words, we show that (6.11) [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.11). That is, z 1 is the unique solution of (6.11). On the other hand, if z 1 > M [i.e., it is not solution of (6.11)], then z 2 > M. That is, z 2 is the unique solution of (6.11). Finally, if z 1 = M, then we have that z 2 = z 1 = M. In any case, (6.11) 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 words, 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 n is tight for any convergent subsequence, there exists a further subsequence, say Zn(∞) n , such that Zn(∞) n d →z * , asn → ∞. We now have that, for any t ≥ 0, and soz * in an invariant point. By the uniqueness of the invariant point, we derivē z * = z * . This concludes the proof.

Proof of Theorem 4.5
We start the analysis by establishing a continuity result, which can be proved by using results in [32].
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 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: (6.14) In order to derive appropriate equations (in the prelimit) for the diffusion-scaled processes, subtract and add the terms n ν+μ ν , n ν+μ ν t in (6.13) and the terms n, nμt, nνt in (6.14), 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 (·): and n (s)ds −Ŷ n (t), (6.16) whereM n r (·) := M r (·) √ n and the scaling for the processŶ n (·) is analogous. The following proposition shows that the processesM r (·) are martingales. Proposition 6.2 Under the assumptions E Z n (0) < ∞ and E Q n (0) < ∞, we have that the processesM n λ n (·),M n μ (·),M n ν,1 (·), andM n ν (·) are square-integrable martingales with respect to the filtration F n := {F n t , t ≥ 0}. Their associated predictable quadratic variations, denoted by · , are 20) and E M n r (t) < ∞, for t ≥ 0 and n ≥ 1.
Proof Fix n ≥ 1. The result for the processM n λ n (·) follows immediately by applying [32, 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 that Since N ν,1 (·) and N ν,2 (·) are independent Poisson processes, [32, Lemma 3.1] implies thatM n ν,i (·) are F n -martingales 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.21), the assumptions E Z n (0) , E Q n (0) < ∞, the inequalities Z n (t) ≤ Z n (0) + N λ n (t), Q n (t) ≤ Q n (0) + N λ n (t), (6.22) and adapting the proof in [32,Lemma 3.4], we obtain that, for fixed n ≥ 1, the following moment conditions are satisfied: Also, again by [32,Lemma 3.4], we derive that the following moments related to the Poisson processes are finite:  is an F n -martingale. Finally, note that by (6.21), asM n μ (·) is adapted to the filtration F n , we obtain thatM 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.,Z  Applying the martingale central limit theorem in [20], we have that where W λ (·), W μ (·), W ν,1 (·), and W ν (·) are (nonindependent) standard Brownian motions. It is essential to observe that, by (6.21) and Remark 6.1, we have that where now W λ (·), W μ (·), W ν,i (·), i = 1, 2, are independent standard Brownian motions. Furthermore, by the properties of Brownian motion, we obtain where WẐ (·) and WQ(·) are nonindependent standard Brownian motions. Further, we have that In addition, by [32,Theorem 7.3], we know that (Q n (·), Y n (·)) satisfies a onedimensional reflection mapping; see [11,Sect. 6.2] for background on the reflection mapping. That is, Y n (·) is the unique nondecreasing nonnegative process such that Q n (t) ≤ K n , (6.15) holds and which is equivalent to

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, it follows that P (G n ) → 1, as n → ∞. In the sequel, we assume that ω ∈ G n .
Proof of Theorem 4.7 Rewriting (4.15) for the n th system and using the assumption P (E n (0) = 0) = 1 yields Let T > 0 and ω ∈ G n . We have that as n → ∞. The result now follows by adapting the proof of Proposition 4.6. Now, taking the limit n → ∞ and using the reflection mapping [11], we derive dQ(t) = − (cμ +Q(t))dt + 2μdWQ(t), where 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 and we can write √ μW (t) where W i are independent standard Brownian motions.

Concluding remarks
This paper proposes modelling an electric vehicle charging system 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 a 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.9) 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 [22]. In addition, we can include in the model time-varying arrival rates, multiple EV types, and multiple parking lots, thus extending results in [34]. Moreover, the distribution grid (low-voltage network) plays a crucial role and it should be included in the model; see [4] for a heuristic approach. For another application in EV charging including the distribution grid, see [10], where simulation results are presented for a Markovian model.