Complete resource pooling of a load balancing policy for a network of battery swapping stations

To reduce carbon emission in the transportation sector, there is currently a steady move taking place to an electrified transportation system. This brings about various issues for which a promising solution involves the construction and operation of a battery swapping infrastructure rather than in-vehicle charging of batteries. In this paper, we study a closed Markovian queueing network that allows for spare batteries under a dynamic arrival policy. We propose a provisioning rule for the capacity levels and show that these lead to near-optimal resource utilization, while guaranteeing good quality-of-service levels for Electric Vehicle (EV) users. Key in the derivations is to prove a state-space collapse result, which in turn implies that performance levels are as good as if there would have been a single station with an aggregated number of resources, thus achieving complete resource pooling.


Introduction
A key challenge in the deployment and take up of electric vehicles by society is the provision of a scalable charging infrastructure. A viable solution is the development of a battery swapping network. Currently, there has been work done on the operation and control of a single battery swapping stations (for example [19]), but there is a clear gap within the literature when extending this to the operation a wider network of stations. In this paper, we introduce a novel stochastic network model describing a network of battery swapping stations which clearly addresses this need and provides a foundation for future studies. In addition, we carry out a detailed analysis of this model and obtained a number of novel insights into the operation of a battery swapping network.
A steady energy transition is taking place due to the de-carbonization of the economy, leading to many intrinsic challenges and research opportunities, of which an overview is given in [14] and [2]. There are numerous challenging problems caused by developments on the demand side. Examples include control problems in local, smart distribution grids, as well as managing increasing demand irregularities caused by e.g. electric vehicles. Modeling the behavior of individual agents and their interaction naturally leads to stochastic models.
Despite the apparent need of alternative energy sources in the transportation sector, the adoption of electrified vehicles has been slow initially due to various practical challenges, such as high purchase costs of an EV, battery life problems and long battery charging times [15]. A possible solution to address these issues is the construction and operation of a battery swapping infrastructure. The upfront costs of purchase of an EV can be significantly reduced when battery swapping station operators own and lease batteries to customers, the batteries can be charged more appropriately to prolong batteries' lifetime [19], and EV users can experience a fast exchange of batteries in contrast to long charging times. Beyond the consumer benefits, the centralized charging paradigm of battery swapping allows the deferment of huge network reinforcement works required to support charging at home by connecting the chargers to the medium voltage network. Furthermore, the aggregation of a large number of batteries at charging stations can provide a comprehensive range of flexibility services to transmission and distribution network service operators.
In this paper, we introduce a model for EVs utilizing battery swapping technology within the context of a fixed region. Within the region there are a number of charging/swapping stations and vehicles, in general, do not leave the region leading to the conservation of batteries. This leads us to introduce a class of closed Markovian queueing network model, which we use in a novel way to model the evolution of the battery population within a city.
With the advancement of smartphones and online technologies, a range of service providers will utilize these advancements to provide occupancy level information to customers to improve delay performance. In a battery swapping system, such information can motivate EV users to move to visit the most appealing location in the direct vicinity. In this paper, we integrate a load-balancing policy to incorporate this. An intrinsic problem is to establish suitable capacity levels that account for the inherent tradeoff between EV users' quality-of-service and operational costs. To the best of our knowledge, this is the first work that considers this question for a battery swapping system in a network framework under a dynamic arrival policy.
Adequately balancing service performance and resource utilization is very much in the spirit of the Quality-and-Efficiency-Driven (QED) regime known from asymptotic many-server queueing theory [10]. Typically, this gives rise to a square-root slack provisioning policy for the capacity levels and has been successfully implemented in many applications such as call centers [4,13,25], healthcare systems [12,24,22] and more. This policy leads to favorable performance for large systems: as the number of customers r grows large, the waiting probability tends to a value strictly between zero and one, the waiting time vanishes with a rate 1/ √ r, and near-optimal resource utilization of 1 − O(1/ √ r) is achieved. To inherit such properties for the battery swapping framework, we adopt a similar capacity level design policy for both the number of charging servers and the number of spare batteries relative to the expected offered load under the load-balancing arrival strategy.
To add to the agreeable properties of delay performance in the QED regime, the arrival strategy ensures that the relative charging loads at the different stations do not grow apart too much since arriving EV users always move to the least-loaded station. This phenomenon has been observed in a number of settings and is referred to as state space collapse, see [6,23] for an overview and [9] for work most closely related to this paper. In fact, when capacity levels are chosen appropriately, this effect is so strong that complete resource pooling takes place: the system behaves as if there is only a single station with an aggregated number of resources. It ensures that it is unlikely that EV users are waiting for a battery at one station, while another is readily available at any other station, even among those stations that are far from his direct vicinity.
The first main contribution of this paper is the introduction of a stochastic model for battery charging in a network setting. In recent years, there has been a growing amount of research on both the planning/design as well as the operation/scheduling in battery swapping systems, see [19] for an overview. Most papers employ robust optimization techniques to find optimal solutions for certain objectives, while little of the works focus on the quality-of-service for EV users. The exception are a collection of papers written by a set of authors [15,18,16,19,17], that use asymptotic analysis and Markov Decision Process techniques to propose suitable solutions. Whereas the focus in those papers is on issues arising in a single station, we propose a network setting to account for queue length correlations between stations.
Our second main contribution involves the novelty of our load-balancing arrival mechanism. Loadbalancing policies have attracted a lot of attention in recent years due to extremely relevant applications in large data centers, see [21] for an excellent literature overview. Typically, these systems comprise many single-server stations where a central dispatcher decides where to allocate incoming tasks. In contrast, our framework involves a network of (a fixed number of) multi-server station for which we introduce a unique feature: an arriving EV user restricts itself to move only to one of the stations in his direct vicinity. By appropriately setting the capacity levels according to the QED provisioning rule, we show that this constraint turns redundant in the sense that the resource pooling effect can still be achieved.
In this paper, we also make several theoretical contributions. Direct analysis of the steady-state distribution of the queue-length process is intractable under the load-balancing strategy in case of multiple stations. Instead, we resort to a fluid and diffusion limit approach. We derive the existence of the fluid limit and point out its unique invariant state. Using a diffusion-scaled queue length process, we zoom in on the fluctuations around the invariant state. We prove a state space collapse (SSC) result by showing that in the limit (as the number of EVs grows larger) the diffusion-scaled queue lengths tend to become arbitrarily close almost instantaneously and stay that way for any fixed interval. This property can be exploited to derive the limiting queue length behavior at every station, and show that it implies the complete resource pooling effect. The derivations of our results rely heavily on the framework developed by Dai and Tezcan [9], that in turn can be seen as an extension of [5]. We adapt their framework to incorporate a closed network setting under the novel load-balancing policy.
The introduction of the novel framework within this paper acts as a foundation for a substantial research programme in the modeling of battery swapping networks. This will provide practitioners with a better understanding of how such networks should be designed and operated from both the perspective of quality of service requirements but also from an economic viewpoint. This can be carried out by enriching the model, here we highlight a few possible directions we consider important and challenging future steps. Each of these will provide a detailed insight a specific aspect of such systems. Firstly, the inclusion of multiple customer types to model a range of car brands within the network using different battery systems. Secondly, there is a delay between the moment an EV user consults queue length information and the actual arrival due to transportation time. As is perceived in health care settings and bike-sharing systems, this can have a considerable effect on the queue length behavior. A third enhancement would be to incorporate a time-inhomogeneous demand rate to better simulate the expected diurnal variation. This will lead to a varying amount of slackness in the capacity within the QED regime. Finally, there is substantial underlying variability in the fluctuating energy prices, which sharply rise whenever the energy grid is more strained and vice versa. A battery swapping infrastructure will be sensitive to these prices changes and can provide an indispensable asset for in supporting a stable grid in the future. Since it can relieve strain during peak moments by deferring the moment of charging or even deplete batteries providing energy to the grid. It is beyond the scope of this paper to provide efficient and adequate provisioning rules in these challenging settings, yet offers intriguing avenues to pursue in future research. The main insight provided in the present study is the effectiveness of simple load balancing policies, and while the model is parsimonious, this insight is useful in, at least, the planning stage of a swapping network.
The remainder of this paper is organized as follows. In Section 2 we describe the battery swapping network and its corresponding load-balancing arrival mechanism. In Section 3 we present the fluid and diffusion results in the special case of a single station, and generalize these results for the multiple stations setting in Section 4. Our results imply approximations for certain performance measures, which we validate through several simulation experiments described in Section 5.

Model Description
In this paper, we consider a queueing network with S battery swapping stations and r EVs. Each EV has one battery (collection) providing the energy for the car to drive. Every station i ∈ {1, ..., S} has three types of assets: F i charging points, B i spare batteries and G i swapping servers. Whenever there is an EV arrival at a station, a swapping server takes out the almost depleted battery and exchanges it for a fullycharged one if available. The swapping time is relatively very short (with respect to charging times), and therefore we assume it to occur instantaneously. Batteries in need of charging are being recharged whenever a charging point is available, and we assume every recharge to take an exponential amount of time with rate µ, independent of everything else. Whenever a battery is fully charged, it is placed in an EV immediately if one is waiting, and otherwise stocked for a future EV arrival. After receiving a fully-charged battery, the EV requires recharging after an exponential amount of time with rate λ. With probability p ij stations i and j are in the EV user's direct vicinity. We assume that EV users consult some online device, and are motivated to move to the station that is relatively least-loaded (ties are broken evenly).
We point out that batteries are always exchanged, and therefore the number of batteries physically present at a station can never be below this station's number of spare batteries. In fact, this observation implies that the queueing model is closed, where the total number of batteries is given by Another observation concerns the role of the swapping servers. Whenever a battery is taken out of the EV, it cannot move from the swapping server until an exchange of batteries has taken place. Therefore, no more than B i + G i batteries can be charged simultaneously at a station i ∈ {1, ..., S}. As a consequence, having more charging points creates no additional charging capacity, and can be bounded as In addition, we assume that the number of such expensive swapping technologies is small at every station, i.e. G i < G for all i = 1..., S, with G < ∞ being small fixed number.
The main quantity of interest in this paper is the number of batteries that are in need of charging, which we will also refer to as the queue length. Let Q i (t) denote the number of batteries in need of charging at station i at time t ≥ 0, and denote Q(t) = (Q 1 (t), ..., Q S (t)). The queue length also defines the occupancy level (load) of a station. More specifically, let p i = S j=1 p ij /2, i = 1, ..., S be the effective arrival probability at station i. An EV in need of charging closest to station i and j at time t ≥ 0 moves to station i iff where ties are broken evenly. Figure 1 illustrates the closed queueing model under this load-balancing arrival mechanism.
To achieve favorable performance levels, we propose an associated QED-scaled capacity level for the resources at the stations. More specifically, we consider a sequence of systems indexed by the number of cars r, where we write a superscript r for processes and quantities to stress the dependency on r. For a system with r cars, we set the capacity levels of the number of charging points and the number of spare batteries as for all i = 1, ..., S. We remark that the bound for the number of charging points originates from (2.1). Since the number of swapping servers is fixed and small and the number of cars r grows large, this condition reduces to the γ ≤ β requirement in (2.3).
Besides the queue length process, we focus on three performance measures in this paper: the waiting probability of an arbitrary EV, its expected waiting time and the resource utilization levels of the stations. As the role of swapping servers is non-existent in this framework, we consider the resources of the swapping stations to be the charging points and the spare batteries. We define the utilization level of the charging points to be the fraction of charging points that are busy with charging, and the utilization level of the spare batteries to be the fraction of batteries that are not fully-charged with respect to the total number of batteries at the station. In steady state, the latter corresponds to the fraction of time at a station that a battery is expected to wait for an arriving EV. Figure 1: Illustration of closed queueing network with multiple stations.

System behavior in case of a single swapping station
When there is only a single battery swapping station, all EVs simply move to this station with probability one. The system reduces to a closed network where batteries are in two possible locations: either positioned in a car in no need of charging, or at the station. The square-root scaling rules reduces to where we suppress the subscript 1 for the station number in this case. An illustration of the closed queueing model is given in Figure 2 The notable advantage of a single station is that all resources are assembled at one entity, and inherently no resources are unavailable for being at a different location. There is also a considerable upside in terms of the analysis: since there is no routing policy anymore, the queue length process becomes a simple birth-death process for which the steady state distribution is easily derived. Yet, the steady-state distribution provides little qualitative insight on the queue length behavior, and in particular, the behavior of the process when it has not reached steady state yet. Therefore, we resort to fluid and diffusion limits which in practice serve as good approximations for moderate to large-scale systems. This allows us to provide approximations for the performance measures of our interest, e.g. the waiting probability and the expected waiting time.

Steady state distribution
As the queue length process is a birth-death process, it is straightforward to derive the steady-state distribution of the queue-length process by standard theory for Markov chains, irrespective whether the QED scaled provisioning rules (3.1) hold. More specifically, the queue length denote the steady state distribution of the number of batteries in need of charging. Lemma 3.1 Suppose S = 1, where the single swapping stations has F r charging points and B r spare batteries. In order to simplify formulas below, we are going to B for B r and F for F r . The steady state distribution is given by

Remark 3.2
In view of (2.1), we exclude the case that F ≥ B in our analysis further on in this paper. Yet, in an application where e.g. G = ∞ and hence F ≥ B possibly holds, we point out that the distribution can be derived similarly. That is, all EVs that arrive at the station find an available swapping server, and the swapping servers do not pose any restriction on the number of batteries that can be charged simultaneously.
Only the number of charging points bounds the charging rate. One can also consider the QED provisioning rule in this case, which we treat in Appendix C. Moreover, if both G = ∞ and F = ∞, Lemma 3.1 shows that Also in this particular case one can pose a QED provisioning rule for the number of spare batteries alone, and derive the asymptotic properties. We treat this case in Appendix B

Limiting queue length behavior
Due to the curse of dimensionality, it is very challenging to gain a qualitative insight in the (transient) behavior of processes in large-scale systems. Therefore, we resort to fluid and diffusion limits to provide good approximations for the behavior in the actual system when r is large. Recall that Q r (t) corresponds to the queue length process (the number of batteries in need of charging) under the scaling rules (3.1) with r cars at time t ≥ 0. We consider the fluid scalinḡ The fluid-scaled process converges to a deterministic, continuous monotone process with a single fixed steady state value.
and has the steady state value lim t→∞Q (t) = λ µ .
Proposition 3.3 implies that the number of batteries in need of charging can be approximated by whereQ(t) = lim r→∞Q r (t) is a solution of an ODE. It describes the approximate (possible) transient behavior before reaching steady state. The proof of Proposition 3.3 is given in Appendix A.1.
We point out that whenever the queue length is near its steady state value, it remains close to its steady state value from that time onward. That is, if Q r (t 0 ) ≈ λr/µ for some t 0 ≥ 0, then Q r (t) ≈ λr/µ for all t ≥ t 0 . From that point on, the fluid limit becomes a rather rough estimate for the number of batteries in need of charging that allows for further investigation on the fluctuations around this value. Therefore, we turn our focus to the diffusion scalinĝ This scaling provides more sensitive approximations, as it captures fluctuations of order √ r. The diffusionscaled process will tend to a piecewise linear Ornstein-Uhlenbeck processes, with a steady state distribution that can be expressed analytically. The proof can be found in Appendix A.1.
Theorem 3.4 Suppose S = 1 and the system operates under (3.1). IfQ r (0) →Q(0) in distribution as r → ∞, then Q r →Q in distribution as r → ∞. The processQ is a diffusion process with drift and constant infinitesimal variance 2µ. The steady state density ofQ(∞) = lim t→∞Q (t) is given bŷ

9)
where α i = r i /(r 1 + r 2 + r 3 ), i = 1, 2, 3 with r 1 = 1, Equation (3.9) in Theorem 3.4 is obtained by taking the limit of the scaled diffusion process (as r → ∞), and finding its steady state distribution (as t → ∞). However, in order to obtain a good approximation of the steady state distribution with a fixed number of cars r, it is arguably more reasonable to consider the steady state distribution of the scaled diffusion process (as t → ∞) and next take the limit as r → ∞.
Fortunately, the following theorem shows that the order at which one takes the limit leads to the same result. The proof of Theorem 3.5 is given in Appendix A.2. As the order at which the limits are taken does not impact the result, we use the limiting processQ(∞) to obtain approximations for the performance measures.

Performance measures
Typical performance measures for the QoS level for the EV users include the waiting probability and the expected waiting time. We view the efficiency-level for the station by the resources utilization. Typically, the QED regime in many-server systems causes the waiting probability to tend to a non-degenerate limit as r → ∞, the waiting time to vanish, while the resource utilization tends to one. These features also appear in our system under the proposed QED scaling.
Due to the PASTA (Poisson Arrivals See Time Averages) property, in open queueing systems where the arrival process is a time-homogeneous Poisson process, the steady state value of any quantity is the same as at arrival instances. In particular, the waiting probability equals the steady state probability that the number of fully-charged batteries is zero, or equivalently, the number of batteries in need of charging is at least B. Unfortunately, the arrival process in our closed setting is state-dependent. Yet, Theorem 3.4 shows that the changes in arrival rate is only of order O( √ r), i.e. the arrival rate is λr − O( √ r) (with high probability).
These small changes will therefore become negligible as r → ∞. In other words, this argument implies that the PASTA property remains valid asymptotically. This notion can be formalized similarly as is done in [11]. Summarizing, if W denotes the waiting time of an arriving EV user, then whereQ(∞) is as in Theorem 3.4.
The key concept to derive the expected waiting time is Little's law, stating that the long-term average number of waiting cars, denoted by Q r W , equals the long-term throughput multiplied by the average waiting time. In other words, E(Q r W ) = θE(W ), where the throughput θ can be viewed as the long-term average rate at which EVs arrive, and hence also leave the battery swapping station. We can express the throughput as θ = λr − λE(Q r W ), since the long-term average number of batteries not in need of charging is in fact the expected number of cars not waiting at the station in this closed system. Therefore, it follows that . (3.10) In turn, the expected number of waiting cars can be derived directly using Theorem 3.4 and the observation as r → ∞. We point out that E(Q r W ) is consequently of order Θ( √ r), and together with (3.10) this implies that E(W ) is of order Θ(1/ √ r) and hence vanishes in the limit.
The resources will be fully utilized under the proposed QED policy as r → ∞. Theorem 3.4 implies that that at most O( √ r) charging points are not utilized, and the number of fully-charged batteries is also of order O( √ r). Therefore, as r → ∞, Theorem 3.6 Suppose S = 1, and the system is operating under (3.1). Then the following properties hold as r → ∞.
The waiting probability has a non-degenerate limit given by The expected waiting time behaves as with α i are as in Theorem 3.4. Finally, the resource utilizations behave as The proof of Theorem 3.6 is given in Appendix A.3.

System behavior in case of multiple stations
When the number of stations S ≥ 2, the analysis of system behavior needs to account for the underlying routing mechanism of arriving EVs. Whenever an EV is in need of recharging, stations i and j are in its direct vicinity with probability p ij , and it chooses to move the station i if (2.2) holds. For a resource pooling effect to occur, we require that there is a sufficient number of pairs (i, j) for which p ij > 0. For example, if the network consists of four stations with p 12 = p 34 = 1/2, there are no arrivals that can choose between one station in the set {1, 2} and another in the set {3, 4}. Therefore, possible discrepancies in queue lengths are not levelled by the arrival mechanism between these two sets. Therefore, we assume that for every non-empty set S of stations, there is at least one pair (i, j) with i ∈ S and j ∈ S for which p ij > 0. This statement is equivalent with the following assumption.
We assume that the graph G is connected.

System dynamics
There are many processes that are of interest in this system, and in particular, the queue length process at each station. In our analysis, we consider {X r (t), t ≥ 0} with is the number of arrivals that are closest to stations i and j until time t ≥ 0 in the r'th system; is the number of arrivals that are closest to stations i and j and are routed to station i until time t ≥ 0 in the r'th system ; is the number of batteries in need of charging at time t ≥ 0 in the r'th system; is the number of busy servers (charging points) at time t ≥ 0 in the r'th system; • Y r , where Y r (t) is the aggregated time of all cars that are not waiting at some station until time t ≥ 0 in the r'th system; is the aggregated time of all servers at station j that were charging until time t ≥ 0 in the r'th system; is the number of service completions at station j until time t ≥ 0 in the r'th system; is the number of batteries that are positioned in an EV not waiting at a station in the r'th system at time t ≥ 0.
Clearly, there are strong relations between the individual processes in X r . For example, there is a routing policy that dictates where a car in need of a full battery drives to in order to swap its battery. This notion is captured by the arrival processes A r (the classification of the different arrival types) and A r d (the routing decision). To generate the arrival and service completion processes, we introduce a set of independent Poisson processes. Let {Λ ij (t), t ≥ 0} for all {i, j} ∈ E be independent Poisson processes with rate p ij λ and {S j (t), t ≥ 1} for all 1 ≤ j ≤ S be independent Poisson processes with rate µ. The system dynamics satisfy the following identities.
We refer to these equations as the system identities, and they prove to be central for deriving our results. The derivations use the framework set out in [9], which in turn is based on [5]. We adopt much of the notation and definitions in this paper, and before stating our main results, we repeat them for the purpose of self-containment. For each positive integer d, we denote by is endowed with the J 1 topology, and the weak convergence in this space is considered with respect to this topology. We say a sequence of functions as n → ∞. Moreover, we say that t ≥ 0 is a regular point of a function x if x is differentiable at t ≥ 0, and denote its derivative by x (·). We assume that the random variables in X r live on the same probability space (Ω, F, P). Often, we consider sample paths of stochastic processes, and whenever we want to make the dependence on the sample path explicit, we write X r (·, ω) for the sample path associated with ω ∈ Ω for a stochastic process X r .

Fluid limit
To capture the rough system dynamics, we consider the fluid-scaled process For each process X r in X r , we define similarly its fluid equivalent asX r = X r /r and its limiting process X = lim r→∞X r . We adopt the definition of a fluid limit and its invariant state(s) from [9]. That is, we consider A ⊂ Ω such that the FSLLN holds, i.e.
Due to the FSLLN, we observe that one can choose A large enough such that P(A) = 1.

Definition 4.2
We callX a fluid limit of {X r } if there exists an ω ∈ A and (sub)sequence {r l } with r l → ∞ as l → ∞, such thatX r l (·, ω) converges u.o.c. toX(·, ω). Moreover, let q = (q 1 , ..., q S ) be an invariant state of the fluid limits if for any fluid limitX, In Proposition 3.3, we focus on the fluid-scaled queue length process only for S = 1, and the sequence r l = l. Instead of requiringQ r (0) →Q(0) withQ(0) a finite constant, Definition 4.2 allows forQ(0) to be random. Proposition 3.3 implies that in case that S = 1, the fluid limits exist and are deterministic, (Lipschitz) continuous paths that depend only on the realization ofQ(0). Moreover, there is a single unique invariant state given by λ/µ. A similar result holds when S ≥ 2.
The (uniqueness of the) invariant state result for the fluid limit is central for the existence of a properlydefined diffusion process as it states that ifQ(0) = (p 1 λ/µ, ..., p S λ/µ), the fluid limits are time invariant. We present a proof of Theorem 4.3 in Appendix D.

Diffusion limit
Due to the policy governing which station a car drives to in order to replace a battery, one observes the so-called load-balancing effect. By setting the number of resources as in (2.3), this load-balancing effect is so strong that in fact complete resource pooling occurs. In other words, the system behaves as if there is a single large swapping station where the number of resources equals the aggregated total of the individual stations. This appealing consequence ensures that there are no idle resources at one station, while at another there are possible long waiting lines of cars that are waiting for a battery exchange.
The key concept to derive this effect is to show a state space collapse (SSC) result. That is, we consider the diffusion-scaled queue length process defined aŝ In our model, the SSC result states that (almost instantaneously) the diffusion-scaled queue length processes are arbitrarily close at all stations, and stay close during any fixed interval.
The proof of Theorem 4.4 is given in Appendix E. This result reveals that instead of considering the individual queue length processes, it suffices to track the total queue length process instead. More specifically, define the sequence of random processes {Q r As the state space collapse implies thatQ r i (t) ≈Q r j (t) for all i, j ∈ {1, ..., S} (for t ≥ K r / √ r), we can approximate the queue length at an individual queue by λr µ for all j = 1, ..., S. The limiting process of the total queue length can be derived using the SSC result.
and constant infinitesimal variance 2µ. The steady state densityQ Σ (∞) is given bŷ where Proof. We observe that the steady state density is a direct consequence of the diffusion process [7]. What remains to show is thatQ r Σ converges to the described diffusion process as r → ∞. Equivalently, we need to show that where {W (t), t ≥ 0} is a standard Brownian motion. We note that due to the system identities, We observe that due to the FCLT and Theorem 4.3, where {BM A (t), t ≥ 0} is Brownian motion with mean zero and variance µ. Similarly, due to the FCLT and Theorem 4.3, where {BM D (t), t ≥ 0} is an (independent) Brownian motion with mean zero and variance µ. The sum of these two processes results is equal to (in distribution) a Brownian with mean zero and variance 2µ, which contributes to the √ 2µ dW (t) term in (4.22). Next, we observe that due to the system identities and the definition of the diffusion scaling, for all t ∈ [0, T ], and due to Theorem 4.4, where the sequence { (r), r ∈ N} can be chosen such that (r) → 0 as r → ∞. Then, as r → ∞, This contributes to the first two terms in (4.22). Applying the continuous mapping theorem concludes the result.
Another consequence of the state space collapse result is that the waiting probabilities and expected waiting times are equal at all stations, as well as the resource utilization levels. In fact, it exhibits the same behavior as if there would be a single station due to the complete resource pooling effect. Corollary 4.6 Suppose the system is operating under (3.1). Then the following properties hold as r → ∞ for all i = 1, ..., S. The waiting probability has a non-degenerate limit given by The expected waiting time behaves as with α i are as in Theorem 4.5. Finally, the resource utilizations behave as

Simulation experiments
The results presented are given in an asymptotic regime where the charging times are exponentially distributed. In this section, we conduct simulation experiments to evaluate the quality of our approximations and the robustness of the state-space collapse result against more general charging times. Throughout the experiments, we consider a network with 5 stations where the arrival probabilities are given by  As a battery swapping infrastructure currently does not exist yet in real-life, there is no (significant) data that can be exploited to obtain useful parameter choices. Instead, we discuss an adequate provisioning strategy under the assumption that recharging takes one hour on average (µ = 1), and that every EV user returns for recharging services after every 40 hours on average (λ = 0.025). In addition, we stress that our results are based on an asymptotic regime, and therefore require the system to be sufficiently large for the approximation to become meaningful. We allow for (at least) r = 50000 EV users in this infrastructure in order for the number of resources to be sufficiently large for the asymptotic results to become meaningful. The effective loads at the stations in this case are p j λr µ and we note that due to the QED provisioning policy in (2.3), the numbers of charging points and spare batteries are close to these values. Obviously, the number of resources are integer values, and in our simulation experiments we choose

State space collapse for exponential charging times
A first-order approximation for the queue length process is implied by the fluid result in Theorem 4.3. We validate this approximation for the above-described setting, with initial queue length Q r i (0) = 150 for all stations i = 1, ..., 5. That is, only station 1 is initially overloaded, while all other station are underloaded. The equations in Theorem 4.3 together with the Lipschitz continuity describe a unique fluid limit with the given initial queue length. This yields the approximations In particular, in the case when the initial queue length is Q r i (0) = 150 for all stations i = 1, ..., 5, this results in approximations where t 1 ≈ 0.1826, t 2 ≈ 0.3189, t 3 = 1.4 and t 4 ≈ 1.4758. The times t i , i = 1, 2, 4 correspond to the times where two stations (approximately) have the same relative queue lengths, and t 3 is the moment where the number of EVs in need of recharging is (approximately) equal to the number of stations/spare batteries. A sample path comparison with its fluid approximation (dotted lines) is graphically illustrated in Figure 4. We observe that the fluid limit approximations capture the typical values of the actual queue length process quite accurately. We observe apparent fluctuations around its approximation, and we note that these become relatively small as r grows large.
To observe the state-space collapse, we plot the same sample path by its diffusion scaling, see Figure 5. Indeed, around t 4 ≈ 1.4758 the diffusion-scaled queue lengths appear to become close and remain nearly equal to one another after this time. In addition, as time moves, on the diffusion-scaled queue lengths fluctuate around zero.
The differences between the diffusion-scaled queue lengths are small and fluctuate erratically among each other when one would zoom in on this domain. Obviously, the differences between the diffusionscaled queue lengths are not arbitrarily small since r is finite. Even if the diffusion-scaled queue lengths at all stations are the same, and an arriving EV moves to station 1, this causes a discrepancy of 1/(p 1 λr/µ) ≈ 0.5657 in the described setting. Theorem 4.4 implies that the distance between the queue lengths become smaller as the number of EV users r grows large. To illustrate this notion, we consider the maximum difference between the queue lengths over a finite interval T = 1 in Figure 6, which is monotonically decreasing in r. In addition, we observe that the average maximum distance, i.e. 1/T is also monotonically decreasing in r, and is not excessively smaller than the maximum distance of the interval.
It turns out that the maximum distance is usually caused by the diffusion-scaled queue length at station 1 being smaller than the queue length at another station (often station 2 or station 3). This discrepancy is also reflected in the waiting time probabilities of the stations. In Figure 7, we plotted the waiting probabilities of all stations in the case of 2500000 EV arrivals averaged over 20 samples. We point out that the stairtype effect appearing in the waiting probabilities is due to the ceiling of the number of resources at the stations. Moreover, as r is finite and we use the ceiling function, the waiting times are not all exactly equal, which is most apparent for station 1. As r grows large, the waiting probabilities do grow closer and move near to their asymptotic expressions. Still, the waiting probabilities are typically below their asymptotic expressions. This implies that the provisioning rules (2.3) guarantee that a desired waiting probability is achieved.

Non-Markovian harging time distribution
In order to be able to rigorously prove the state-space and consequential results, we assumed exponential charging times in our framework. Yet, extensive simulation experiments suggest that these results hold for any charging time distribution with finite mean and variance. In Figure 8 and 9, it appears that similar   behavior occurs on fluid scale in case of deterministic and uniformly distributed charging times as for the exponential case.
When the queue lengths are initially zero, the system behaves close to its invariant state for t ≥ 1. Similarly to the setting with exponential charging times, the maximum difference between the diffusion-scaled queue length behaves quite erratically, see Figure 10. Still, the differences are very small, and grow smaller as r grows larger suggesting that state-space collapse also holds in this setting. That is, the system behaves similarly to the situation when there would be a single station with an aggregated number of charging points and spare batteries and a charging time distribution as at the individual stations. Consequently, performance measures such as waiting probability and expected waiting time are approximately equal to their equivalents in a single-station system.

A Fluid and diffusion limit proofs for single station system
In this appendix, we provide the proofs for the results in case that S = 1.

A.1 Fluid and diffusion limits
First, we derive the fluid limit.
The second result involves the diffusion limit.
Using Equations (28) and (33) in [7], this implies that the limiting process is a piecewise-linear diffusion process with steady state density given by (3.9), where α i , i = 1, 2, 3 is the probability that the steady state process is in that interval. Equations (5) and (6) from [7] yield the steady state probabilities.

A.2 Steady state limits
Proof. The first expression follows from the central limit theorem and the properties of the Poisson distribution, The second expression can be obtained using geometric series, The final expression again uses the central limit theorem and the properties of the Poisson distribution. That is, Since for every z ∈ R, To obtain the limiting distribution of the steady state distribution as r → ∞, it remains to derive the asymptotic behavior of the probabilities that the queue length process is in the intervals [0, F r ], [F r , B r ], and [B r , B r + r]. To do so, we use the following three lemmas.  Proof. This follows from the central limit and the properties of the Poisson distribution, λr µ k k! e −λr/µ = P Pois(λr/µ) < λr/µ + γ λr/µ → Φ(γ).

Proof. It follows from Lemmas A.2-A.4 that
and applying Lemmas A.2-A.4 again yields , .
Consequently, it follows that Theorem 3.5 holds.

A.3 Performance measures
Using the results of Theorem 3.4, one can derive the asymptotic behavior of the performance measures, as described in Section 3.3.
Proof of Theorem 3.6. Since the PASTA property holds asymptotically, we find that the probability an arriving car has to wait for a fully-charged battery is asymptotically equivalent to the probability that the system is in a state where more than B r batteries are charging at the same time, which gives the first asymptotic relation. Then, where r i , i = 1, 2, 3 are as in Theorem 3.4, from which the result follows directly.
For the expected waiting time, we consider the expected number of waiting cars at the charging station. We observe that where α 3 is as in Theorem 3.4. We remark that this expression is always positive: trivially for γ ≤ 0, and also for γ > 0 since for every x > 0 which can be shown by partial integration. Due to Little's law, we note that which yields the result for the expected waiting time.
Finally, for the utilization levels, the scaling ofQ(∞) implies that Therefore, the utilization levels satisfy (3.11) and the result follows.

B.1 Fluid and diffusion limits
Although our main focus is performance measures in steady state, it is instructive to see the transient behavior for this large-scale system. To that end, we explore the fluid and diffusion limits.
Consider the fluid-scaled process {Q r (t), t ≥ 0}, whereQ r (t) = Q r (t)/r for all r ≥ 1. The following theorem holds under the QED scaling.
Theorem B.1 If S = 1 andQ r (0) →Q(0) in distribution withQ(0) a finite constant, thenQ r →Q in distribution as r → ∞, whereQ satisfies the ODE and has the unique steady state valueQ Proof. We follow the framework of [7] for birth-death processes. We observe that the Q r is a birth-death process with state space Q r = {0, 1, ..., B r + r}, and arrival rates λ r (j) = λ (r − (j − B r ) + ) and service rate µ r (j) = jµ for all j ∈ Q r . The fluid-scaled processQ r therefore has state space {0, 1/r, ..., (B r + r)/r}, and drift and diffusion functions Taking the limit yields That is, we obtain a degenerate limiting diffusion (since the limiting infinitesimal variance is zero), and hence the limiting initial point will yield a deterministic path for the limiting processQ. The mean directly provides the ODE that the limiting processQ satisfies. What remains to show is the limiting value (as t → ∞) .
We now zoom in on the fluctuations around the fluid limit by taking the diffusion limit. We consider the centered scaled processQ Under the QED scaling rule (B.1), we find the following diffusion limiting process.
Theorem B.2 Suppose the system operates under the QED scaling rule (B.1). IfQ r (0) →Q(0) in distribution as r → ∞, thenQ r →Q in distribution as r → ∞, whereQ can be described as a two pieced-together Ornstein-Uhlenbeck (OU) process. More precisely,Q is a diffusion process with mean and constant infinitesimal variance 2µ. The steady state densityQ(∞) is given bŷ Proof. The proof is similar to the techniques used in the proof of the fluid limits. The infinitesimal mean of the centered scaled process is given by λ r λr/µ + x λr/µ + µ r λr/µ + x λr/µ λr/µ → 2λr λr/µ = 2µ, as r → ∞. Using Equation (28) from [7], this implies that the limiting process is a piecewise-linear diffusion process with steady state density given as in (B.3), where α is the probability that the steady state process is above β. Equations (5) and (6) from [7] yield

B.2 Steady state limits
We point out that the fluid and diffusion limits are useful for studying the system in transience. The steady state limiting distribution is then obtained by setting t → ∞ with respect to the limiting diffusion process. In other words, the steady state density function (B.3) is obtained by first taking the limit of the centered scaled process as n → ∞, after which its steady state distribution is obtained by setting t → ∞. The following result shows that interchanging the limits does not change the result. Before moving to the proof of this result, we need to understand the limiting behavior of the steady state probability that more than B r batteries are in need of recharging. To do so, we first determine the asymptotic behavior of the normalization constant in (3.6).

Lemma B.4 Under scaling rule
The convergence of the first term follows directly from the CLT for a Poisson distributed random variable, We extract from the second term Bin(B r + r, λ/(λ + µ)) − (B r + r) λ λ+µ (B r + r) λµ where we used the CLT for a sum of independent Bernoulli distributed random variables. Finally, what is left of the second term is given by where we used Stirling's approximation twice and a Taylor expansion for the logarithm term. Combining the three expressions yields the results.
Using Lemma B.4, we can derive the asymptotic behavior of having more than B r batteries in need of recharging in steady state.

Proposition B.5
Under scaling rule (B.1), the probability that an arriving car has to wait for a battery behaves as Proof. Note that It follows from the proof of Lemma B.4 that, as r → ∞, Rewriting the terms concludes the result.

Proof of Theorem B.3. By definition,Q
and we want to show that changing the order of the limits yields the same distribution. The steady state distribution for a system with r cars is given by (3.5). Due to Proposition B.5, it suffices to show that for every x < β and for every x ≥ β, Using the steady state distribution given in (3.5), we observe that for every β, as r → ∞. For every x ≥ β, Similarly to the proof of Proposition B.5, The CLT implies and

B.3 Performance measures
We consider the waiting probability, the waiting time and the utilization level of the spare batteries. We note that α in (B.3) corresponds to the probability that {Q(∞) ≥ β}. Although the PASTA property does not hold in this closed system, it does hold asymptotically, providing an appropriate approximation for the waiting probability.

Corollary B.6
Under scaling rule (B.1), as r → ∞, For the waiting time, we point out that in this case the identity still holds due to Little's law. Using this relation, we can obtain the expected waiting time.

Proposition B.7
Under scaling rule (B.1), as r → ∞, Consequently, the waiting time of a car behaves as Proof. We note that as r → ∞,

Elementary calculus yields
Then as r → ∞, we obtain .
Simply rewriting the terms yields the result.
Finally, the utilization level follows directly from the scaling of the diffusion scaled process. More specifically, the scaling implies that E(# Fully-Charged Batteries) = Θ( √ r) and E(Q r ) = Θ(r).

C Unlimited number of swapping servers for single station system
In this appendix, we consider a setting with a single station where the swapping servers pose no condition, i.e. G = ∞ or equivalently G ≥ max{1, F r − B r }. That is, every battery of an EV that is at the station can be taken to a charging point (if available). We point out that for all F r ≤ B r , this reduces to the main setting in this paper. Therefore, in this section, we consider only

C.1 Diffusion limits
The steady state distribution is given in Lemma 3.1. Moreover, the fluid limit result of Proposition 3.3 and its proof remain valid in this case, since the fluctuations of order Θ( √ r) are too small to be of any impact on the fluid scale. There will however be a notable effect on the diffusion scale.
Using Equation (28) and (33) from [7], this implies that the limiting process is a piecewise-linear diffusion process with steady state density given as in (C.3), where α i , i = 1, 2, 3 is the probability that the steady state process is in that interval. These steady state probabilities follow directly from Equations (5) and (6) in [7].

C.2 Steady state limits
The main goal of this section is to show that the limiting steady state distribution is the same process as the diffusion limiting process in steady state. Again, we observe that the form of the formula is different in the intervals [0, F ], [B r , F ], and [F, B r + r]. First we derive the limiting steady state distribution conditioned to be in one of these intervals.

Lemma C.2
Under scaling rules (3.1), as r → ∞, x < β, Proof. The first expression follows, similarly to the case of F r < B r , from the CLT and the properties of the Poisson distribution: The second expression follows as in the proof in Proposition 14, using the binomial distribution and the CLT: The final expression again uses the CLT and the properties of the Poisson distribution: Since for every z ∈ R, as r → ∞, the result follows immediately.

Lemma C.4
Under scaling rules (3.1), as r → ∞, Proof. This follows from the CLT and the properties of the Poisson distribution.
Lemma C.5 Under scaling rules (3.1), as r → ∞, Proof. The proof is similar to that of Lemma B.4, using the properties of the binomial distribution and the CLT, and Stirling's approximation.
Lemma C.6 Under scaling rules (3.1), as r → ∞, Proof. We observe that Moreover, applying Stirling's approximation twice, Multiplying the two expressions concludes the result.

C.3 Performance measures
Theorem C.7 Under scaling rules (3.1), The expected waiting time behaves as with α i as in Theorem C.1. Finally, the utilization levels behave as

Proof.
Since the PASTA property holds asymptotically, we find that the probability an arriving car has to wait for a fully-charged battery is asymptotically equal to the probability that the system is in a state where more than B r batteries are charging at the same time, which gives the first asymptotic relation. Thus, where r i , i = 1, 2, 3 are as in Theorem C.1, from which the result follows directly.
For the expected waiting time, note that Due to Little's law, we have the identity which yield the results for the expected waiting time.
Therefore, the utilization levels satisfy (3.11) and the result follows.

D Fluid limit proofs
The proof of Theorem 4.3 is similar to the proof of [9], but adapted appropriately to our system.
Proof. First, we show that the fluid limits exist, where all components are Lipschitz continuous. For this purpose, we show that for all ω ∈ A the sequence {X(·, ω)} has a convergent subsequence, where every component in the limiting process is Lipschitz continuous.
Next, we consider the arrival processes. First, we observe that L r (t) ≤ r for every t ≥ 0. Therefore, for all 0 ≤ t 1 < t 2 , and hence there exists a subsequence {r l } such thatȲ r l (·, ω) converges u.o.c. to somē Y (·, ω) as l → ∞ for every j = 1, ..., S, which is again Lipschitz continuous.
Fluid equations (4.10)-(4.17) follow from the FSLLN results and applying Lemma 11 of [1]. Equation (4.18) requires additional arguments. SupposeX to be a fluid limit with corresponding ω ∈ A and subsequence {r l } l∈N . If for some t > 0 we have thatQ j (t)/p j >Q i (t)/p i , then it follows by the continuity of the fluid limit that there exists a δ > 0 such that . By definition of the fluid limit, it holds for large enough r l that In this case, the routing policy states that all arrivals of type {i, j} move to station station i. Therefore, A r l ij,j remains constant on [t − δ, t + δ], and henceĀ ij,j (t) = 0. Moreover, station i receives all arrivals and by the FSLLN and (4.11), for all t 1 < t 2 with t 1 , t 2 ∈ [t − δ, t + δ]. It follows thatĀ ij,i (t, ω) = p ij λL(t, ω) by (4.1).
Finally, we show that there is a unique invariant state given by q = (p 1 λ/µ, ..., p S λ/µ). Introduce the function and writeS Trivially, h(t) ≥ 0. SinceQ(·) is Lipschitz continuous, so is h(·), and hence it is differentiable almost everywhere. To show that h(0) = 0 implies h(t) = 0 for all t ≥ 0, it therefore suffices to show that if h(t) > 0 then h (t) < 0 for every regular point t ofX. By (4.12), we observe that for every j = 1, ..., S and regular point Due to (4.13)-(4.16),D j (t) = min{µQ j (t), p j λ}. In particular, we observe thatD i (t)/p i =D j (t)/p j for all i, j ∈S max (t), as well as for all i, j ∈S min (t). Due to Lemma 2.8.6 from [8], as t is a regular point, it follows that i:{i,j}∈EĀ ij,j (t) p j = k:{k,l}∈EĀ kl,l (t) p l for every j, l ∈S max (t), as well as for every j, l ∈S min (t). Due to (4.10), (4.11) and (4.18), we conclude that for every j ∈ S max (t), On the other hand, for every j ∈ S min (t), Observing thatD j (t)/p j >D i (t)/p i for every j ∈ S max (t) and i ∈ S min (t), we therefore conclude that if h(t) > 0 with t a regular point, then In other words, for every invariant state of the fluid limit, it must hold thatQ i /p i (t) =Q j /p j (t) for every i = j. We observe, in view of (4.12), that where for every 1 ≤ j ≤ S,D j (t) p j = µ min{Q j (t)/p j , λ/µ}, and hence in view of (4.10) and (4.11), That is,Q j (t) = 0 if and only if µ min{Q j (t)/p j , λ/µ} = p j λL(t).

E State space collapse proofs
To prove Theorem 4.4, we use a framework similar to that of [5], and [9]. The construction consists of several steps, which we lay out next.
1. Divide interval [0, T ] into T √ r intervals of length 1/ √ r, indexed by m. In each interval, consider the hydrodynamically-scaled process of X. For each of these interval, we (a) show the scaled process is "almost" Lipschitz continuous; (b) show convergence to some hydrodynamic limiting process for a sufficiently large part of the state space; (c) derive the hydrodynamic limit equations.
2. Relate the hydrodynamic scaling to the diffusion scaling, using a SSC function to deal with complications regarding the range of the time variable. Transferring the results appropriately, we show multiplicative SSC with respect to the SSC function.
3. Using a compact containment condition, we show that this implies strong SSC.

E.1 Hydrodynamic scaling and its limiting process
In order to introduce the hydrodynamic scaling, we use a diffusion scaling for the values of the process but we slow the process down in time in order to analyze what occurs initially (what would happen instantaneously on diffusive scale). That is, we divide the interval [0, T ] in T √ r intervals of length 1/ √ r, indexed by m. We write p = (p 1 , ..., p S ), and For the processes in X, we introduce the following hydronimically-scaled variants. For Q r , Z r and L r , let the deviations of these processes with respect to their fluid limits. For the processes A r , A d , Y r , T r and D r , we introduce In other words, we track the increase of these processes during the interval [m/ √ r, m/ √ r + √ x r,m t/r]. By definition of x r,m , we note that |X r,m (0)| ≤ 1, which will be a required compactness property when we prove convergence to a hydrodynamic limit. Moreover, due to our fluid limit results, we can show that √ x r,m /r is very small for all ω ∈ A. Proof. Due to our fluid limit result in Theorem 4.3 and the definition of x r,m in (E.1), we observe that for r large enough. Moreover, for r large enough, We conclude that for every m ≤ √ rT , For the hydrodynamically scaled process, the system identities translate to A r,m ij,i (t) can only increase when In order to show that X r,m is almost (with the exception of certain events) Lipschitz continuous, we would like to exclude these certain events, i.e. show that such events are unlikely to occur. Lemma E.2 Fix > 0, M > 0 and T > 0. For r large enough, there exists a constant N > 0 (only depending on λ, µ, and {p ij ; {i, j} ∈ E}) such that Proof. First, we note that due to the memoryless property which both the arrival and service completion processes satisfy, the choice of m is irrelevant and thus can be made arbitrarily. To prove (E.19), we start by showing that for every {i, j} ∈ E, From (E.14) and (E.17), we observe that Y r,m (t) ≤ t and non-decreasing. Due to the properties of Poisson processes, Therefore, where the final to last inequality follows from Proposition 4.3 in [5]. Choosing N = 1/(λ min {i,j}∈E p ij ) and applying the union bound twice, we obtain which yields (E.19).
The proof for (E.20) is completely analogous, but with minor adaptions as one uses T r,m j (t) ≤ p j λ/µt instead of Y r,m (t) ≤ t. We conclude that (E. 19) and (E.20) show that the hydrodynamically scaled arrival process and service completion process are almost Lipschitz continuous.
In order to prove (E.21) and (E.22), we introduce the following processes. Let {u ij (l), l ≥ 1} be independent exponentially distributed random variables with rate p ij λ, representing the time that a car has before it needs recharging at either station i or j. Let {v j (l), l ≥ 1} be independent exponentially distributed random variables with rate µ, representing the service requirement (recharging time) of a battery at station j. Define the aggregated interarrival time of n cars that will choose between stations i and j, and the total service requirement of n batteries at station j, respectively. We observe the identities Λ ij (t) = max{n : U ij (n) ≤ t}, S j (t) = max{n : V j (n) ≤ t}, t ≥ 0.
Using the previous result, we can show that X is almost Lipschitz continuous.
Proof. This follows in a straightforward way from Lemma E.2 and the hydrondynamically-scaled system equations. That is, let V r denote the intersection of the complements of the events given in equations (E.19)-(E.22), so P(V r ) ≤ 1 − N 0 with N 0 the number of equations in Lemma E.2. We note that in order to prove the proposition, it suffices to show that for every ω ∈ V r , and for every t 1 , t 2 ∈ [0, T ] and m < √ rT , where N 1 and N 2 are only dependent on the system parameters (i.e. λ, µ, p). Let t 1 , t 2 ∈ [0, T ] with t 1 ≤ t 2 . By definition of V r , for N as in Lemma E.2. Due to (E.10), In view of (E.12) and (E.10), we observe Due to (E.21), and similarly, due to (E. 22), In view of (E. 16), Finally, due to (E.17), We conclude that (E.31) is satisfied, as each process in X r,m satisfies this property.
As is done in [5,9], one can take appropriately small for every system. That is, for fixed M > 0, N > 0 and T > 0, let where (r) → 0 as r → ∞ is a sequence of positive real numbers. Moreover, in view of Lemma E.1, let for that same sequence { (r)} r∈R , Let K r denote the intersection of K r 0 , H r , and the complements of the events in Lemma E.2. We note that Lemma E.1, Lemma E.2 and Proposition E.3 continue to hold for the sequence (r) if (r) → 0 sufficiently slow. We conclude that P(K r ) → 1 as r → ∞. Following the framework of [5], we can use these results to state that the hydrodynamically scaled system convergences to a hydrodynamic limit. Fix M > 0 and letẼ be the set of right-continuous functions Moreover, we set We remark that these definitions are not related to E, the set of all possible pairs of stations where cars can move to.
Definition E.5 A hydrodynamic limit of E is a point x ∈Ẽ such that for all > 0 there exists a r 0 ∈ N so that for every r ≥ r 0 there is some y ∈ E r such that |x(·) − y(·)| M < .
Finally, to conclude this section, we derive the equations that are satisfied by any hydrodynamic limit.
Proposition E.7 Let M > 0 be fixed, and letX be a hydrodynamic limit of E over [0, M ]. ThenX satisfies the following equations:Ã Remark E. 8 We cannot provide such general equations forZ(·) orL(·), since these limits depend on x r,m . That is, the processesZ r,m (·) andL r,m (·) converge to a limit, but the limiting process may differ for different m. In the proof, we specify the limiting equations of these processes as well.
To show that g(Q(t)) is bounded by this function, we note that it suffices to show that whenever g(Q(t)) > 0 with t ≥ 0 being a regular point ofX, For this purpose, let Due to Lemma 2.8.6 in [8], it holds for all i, j ∈ S max (t) thatQ i (t)/p i =Q j (t)/p j , and similarly, for all i, j ∈ S min (t) it holds thatQ i (t)/p i =Q j (t)/p j . Therefore, due to hydrodynamic limit equations (E.32)-(E.38) and the observation that there is at least one station j ∈ S max (t) such that {i, j} ∈ E, it follows thatQ Similarly, for all i ∈ S min , We conclude that g (Q(t)) ≤ −h, and hence g(Q(t)) ≤ HX(t) for all t ≥ 0. In particular, if g(Q(0)) = 0 and |X(0)| ≤ 1, it follows from the definition of HX(·) that g(Q(t)) = 0 for all t ≥ 0.
This result implies that the hydrodynamically scaled queue length almost satisfies this property as well. The next result is an immediate consequence of Corollary E.6. Corollary E.10 Fix > 0, M > 0 and T > 0. Then, for every ω ∈ K r , g(Q r,m (t)) ≤ H(t) + for all t ∈ [0, M ] and m < √ rT , where H(·) is as in Lemma E.9. Moreover, if g(Q(0)) → 0 in probability as r → ∞, then for all ω ∈ L r with L r = K r ∩ g(Q r,0 (0)) ≤ it holds that g(Q r,0 (t)) M ≤ , and lim r→∞ P(L r ) = 1.

E.3 Multiplicative state space collapse
The goal of this section is to show multiplicative state space collapse for the SSC function defined in (E.39).
To do so, we first need to relate the hydrodynamic and diffusion scaling. That is, we observe that where y r,m = x r,m r = max p λ µQ Corollary E.10 can be translated to the diffusion scaled process. Consider the SSC functionĝ : R S → R defined asĝ with q = (q 1 , ..., q S ).
Since H(·) is given as in (E.40), we observe that H(t) = 0 for all t ≥ 2/(h min 1≤j≤S p j ). We would like to show that ( √ rt − m)/y r,m can be chosen large enough to obtain a very small upper bound, and use that property to show multiplicative state space collapse.
For every t ∈ (M y r,0 / √ r, T ], it follows by definition of m r (t) that √ rt ≥ m r (t) − 1 + y r,mr(t)−1 M.
In particular, where the last inequality follows since M ≥ 2(N + 2).
Next, we show the main result of this section.
Proof. Fix η > 0 and note that by construction there exists a r 0 such that for all r > r 0 , Note that for every t ≥ 0, Moreover, since > 0 is arbitrary, we can conclude that (E.41) holds.

E.4 Strong state space collapse
Although Theorem E.13 shows multiplicative state space collapse, our interest lies in the strong state space collapse as is stated in Theorem 4.4. To do so, it suffices to show that the denominators in Theorem E.13 are bounded in a probabilistic sense. More specifically, Q r (t) T should satisfy the compact containment property. Before doing so, we prove a result that shows that even if the diffusion-scaled queue lengths are initially not necessarily close to one another, these queue lengths do not explode for a sufficiently short period of time. Proof. Fix ∈ (0, 1) small. First, note that P Q r (t) M yr,0/ √ r > K ≤ P Q r (t) M yr,0/ √ r > K; max |Q r (0)|, y r,0 ≤ K + P max |Q r (0)|, y r,0 > K . To bound the first term in (E.46) as well, we observe that the queue length at some time is trivially bounded by We observe that F r j ≤ (1 + )λr/µ for r large enough. Moreover, if {Λ(t), t ≥ 0} denotes a Poisson process with rate λ, then due to the properties of the Poisson process it holds that {i,j}∈E Λ ij (·) d = Λ(·). In terms of the diffusion scaling, the above bound yields for all t ≥ 0, for > 0 small enough (e.g. for < min 1≤j≤S p j /(M λ/µ + min 1≤j≤S p j )). We conclude that the first term in (E.46) converges to zero, i.e.
Next, we show that the processQ r (·) satisfies the compact containment property. Clearly, P Q r (t) T > K ≤ P (τ r K ≤τ r K ≤ T ) + P (τ r K ≤τ r K ≤ T ) .

(E.48)
In order to improve the readability of the proof, we first present a proof in the case whenĝ(Q r (0)) → 0 in probability as r → ∞. We then comment on the changes needed to adapt the proof to the case when this condition does not necessarily hold.
(E.50) Due to the system identities, we observe that for every t ∈ [T r K ,τ r K ], We note that due to the properties of the Poisson process, with {Λ(t), t ≥ 0} an (independent) Poisson process with rate λ. Moreover, since for all t ∈ [T r K ,τ r K ] it holds thatQ r j (t) ≥ γ for every i ∈ {1, ..., S}, Using the FCLT, we observe that as r → ∞, where {BM(t), t ≥ 0} is a Brownian motion with zero mean and finite variance (independent of K). Finally, by the definitions of the stopping times, and in view of (E.49) and (E.50), for all t ∈ [T r K ,τ r K ], We conclude that lim r→∞ P τ r K ≤τ r K ≤ T ; ĝ(Q r (t)) τ r K max{ Qr (t) τ r K , 1} ≤ ≤ P sup 0≤s≤t≤T BM(t) − BM(s) − γµ(t − s) ≥ 1 − 3 2 K , which converges to zero as K → ∞ since ∈ (0, 1/3).
Again, due to Theorem E.13 and (E.45), lim r→∞ P g(Q r (t)) τ r K max{ Qr (t) τ r K , 1} Due to the system identities, we observe that for every t ∈ [T r K ,τ r K ], Q r Σ (t) =Q r Σ (T r K ) + {i,j}∈E A r ij (t) − A r ij (T r K ) λr/µ − S j=1 D r j (t) − D r j (T r K ) λr/µ .
Due to the definitions of the stopping times, we observe that for all t ∈ [T r K ,τ r K ], it holds thatQ r j (t) ≤ β for every j ∈ {1, ..., S}, and hence L r (t) = r. Therefore, due to the properties of the Poisson process, with {Λ(t), t ≥ 0} a Poisson process with rate λ. Moreover, Using the FCLT, we observe again that as r → ∞, where we recall that {BM(t), t ≥ 0} is a Brownian motion with zero mean and finite variance (independent of K). Finally, by definition of the stopping times, and in view of (E.49) and (E.50), it holds for all t ∈ [T r K ,τ We conclude that which also converges to zero as K → ∞ since ∈ (0, 1/3). Since this holds for both of the two summed probabilities in (E.48), we conclude that (E.47) holds.
Case II: general case, i.e. when we do not assume thatĝ(Q r (0)) → 0 in probability as r → ∞.
Therefore, we can assume that bothτ r K >T r K > M y r,0 / √ r andτ r K >T r K > M y r,0 / √ r (for K large enough).
The proof in this general case is then completely analogous to that in the previous one: ĝ(Q r (t)) τ r K needs to be replaced with sup t∈(M yr,0/ √ r,τ r K ] |ĝ(Q r (t))|, and ĝ(Q r (t)) τ r K -with sup t∈(M yr,0/ √ r,τ r K ] |ĝ(Q r (t))|. Next, we prove our main result stated as in Theorem 4.4.