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

TL;DR

This paper analyzes a battery swapping network for electric vehicles, proposing a capacity provisioning rule that ensures near-optimal resource utilization and quality of service through complete resource pooling.

Contribution

It introduces a capacity provisioning rule for battery swapping stations that achieves complete resource pooling and near-optimal performance in a Markovian queueing model.

Findings

Capacity provisioning leads to near-optimal resource utilization.

State-space collapse implies performance as if resources are pooled.

Guarantees good quality-of-service levels for EV users.

Abstract

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.

Equations

$$\text{Total \# batteries in system}=r+\sum _{j=1}^S{B}_i.$$

$$F_i\le B_i+G_i,i=1,\dots ,S.$$

$$\displaystyle\left\{\begin{array}[]{ll}B_{i}^{r}=p_{i}\left(\frac{\lambda r}{\mu}+\beta\sqrt{\frac{\lambda r}{\mu}}\right),&\beta\in\mathbb{R},\\ F_{i}^{r}=p_{i}\left(\frac{\lambda r}{\mu}+\gamma\sqrt{\frac{\lambda r}{\mu}}\right),&\gamma\leq\beta,\end{array}\right.$$

$$\frac{Q_i(t)}{p_i}<\frac{Q_j(t)}{p_j},$$

$$p_i=\sum _{m\in \mathcal{M}}\text{mathds}{1}_{\{i\in m\}}{s}_m/|m|.$$

$$\begin{array}{rl}{B}^r& =\frac{\lambda r}{\mu }+\beta \sqrt{\frac{\lambda r}{\mu }},\beta \in \mathbb{R},\\ F^r& =\frac{\lambda r}{\mu }+\gamma \sqrt{\frac{\lambda r}{\mu }},\gamma \le \beta ,\end{array}$$

$${\pi }_k^{(B^r,F^r,r)}=\mathbb{P}\left(Q(\infty )=k\right)$$

$$\displaystyle\pi_{k}^{(B,F,r)}=\left\{\begin{array}[]{ll}\frac{(\lambda r/\mu)^{k}}{k!}\pi_{0}^{(B,F,r)}&\textrm{if }0\leq k\leq\min\{B,F\},\\ \frac{\left(\lambda r/\mu\right)^{k}}{F!F^{k-F}}\pi_{0}^{(B,F,r)}&\textrm{if }F\leq k\leq B,\\ \frac{r^{B}r!}{(r+B-k)!}\frac{\left(\lambda/\mu\right)^{k}}{k!}\pi_{0}^{(B,F,r)}&\textrm{if }B\leq k\leq F,\\ \frac{r^{B}r!}{(r+B-k)!}\frac{\left(\lambda/\mu\right)^{k}}{F!F^{k-F}}\pi_{0}^{(B,F,r)}&\textrm{if }\max\{B,F\}\leq k\leq B+r,\end{array}\right.$$

$${\pi }_0^{(B,F,r)}={\left(\sum _{k=0}^F\frac{(\lambda r/\mu {)}^k}{k!}+\sum _{k=F+1}^{B-1}\frac{{\left(\lambda r/\mu \right)}^k}{F!F^{k-F}}+\sum _{k=B}^{B+r}\frac{r^{B}r!}{(r+B-k)!}\frac{{\left(\lambda /\mu \right)}^k}{F!F^{k-F}}\right)}^{-1}$$

$${\pi }_0^{(B,F,r)}={\left(\sum _{k=0}^B\frac{(\lambda r/\mu {)}^k}{k!}+\sum _{k=B+1}^{F-1}\frac{r^{B}r!}{(r+B-k)!}\frac{{\left(\lambda /\mu \right)}^k}{k!}+\sum _{k=F}^{B+r}\frac{r^{B}r!}{(r+B-k)!}\frac{{\left(\lambda /\mu \right)}^k}{F!F^{k-F}}\right)}^{-1}$$

$$\displaystyle\pi_{k}^{(B,r)}=\left\{\begin{array}[]{ll}\frac{(\lambda r/\mu)^{k}}{k!}\pi_{0}^{(B,r)}&\textrm{if }0\leq k\leq B\\ \frac{r!r^{B}}{(r+B)!}\binom{r+B}{k}\left(\frac{\lambda}{\mu}\right)^{k}\pi_{0}^{(B,r)}&\textrm{if }B\leq k\leq B+r,\end{array}\right.$$

$${\pi }_0^{(B,r)}={\left(\sum _{k=0}^{B-1}\frac{(\lambda r/\mu {)}^k}{k!}+\frac{r!r^B}{(r+B)!}\sum _{k=B}^{B+r}\left(\genfrac{}{}{0ex}{}{r+B}{k}\right){\left(\frac{\lambda }{\mu }\right)}^k\right)}^{-1}$$

$${\bar{Q}}^r(t)=\frac{Q^r(t)}{r},r\ge 1,t\ge 0.$$

$$\displaystyle\frac{d\bar{Q}(t)}{dt}=\left\{\begin{array}[]{ll}\lambda-\mu\bar{Q}(t)&\textrm{if }\bar{Q}(t)<\lambda/\mu,\\ \lambda^{2}/\mu-\lambda\bar{Q}(t)&\textrm{if }\bar{Q}(t)\geq\lambda/\mu\end{array}\right.$$

$$\underset{t\to \infty }{\lim }\bar{Q}(t)=\frac{\lambda }{\mu }.$$

$$Q^r(t)\approx r\bar{Q}(t),$$

$${\hat{Q}}^r(t)=\frac{Q^r(t)-\lambda r/\mu }{\sqrt{\lambda r/\mu }},r\ge 1,t\ge 0.$$

$$m(x)=-\lambda (x-\beta {)}^{+}-\mu \min \{x,\gamma \},$$

$$\displaystyle\hat{f}(x)=\left\{\begin{array}[]{ll}\alpha_{1}\frac{\phi(x)}{\Phi(\gamma)}&\textrm{if }x<\gamma,\\ \alpha_{2}\left(\gamma e^{-\gamma(x-\gamma)}\right)\left(1-e^{-\gamma(\beta-\gamma)}\right)^{-1}&\textrm{if }\gamma\leq x<\beta,\\ \alpha_{3}\sqrt{\frac{\lambda}{\mu}}\phi\left(\frac{x-(\beta-\frac{\mu}{\lambda}\gamma)}{\sqrt{\mu/\lambda}}\right)\Phi\left(-\sqrt{\frac{\mu}{\lambda}}\gamma\right)^{-1}&\textrm{if }x\geq\beta,\\ \end{array}\right.$$

$$r_1$$

$$r_2$$

$$r_3$$

$$\mathbb{P}(W>0)=\underset{r\to \infty }{\lim }\mathbb{P}\left(Q^r(\infty )\ge B^r\right)=\mathbb{P}\left(\hat{Q}(\infty )\ge \beta \right),$$

$$\mathbb{E}(Q_W^r)=\theta \mathbb{E}(W),$$

$$\theta =\lambda r-\lambda \mathbb{E}(Q_W^r),$$

$$\mathbb{E}(W)=\frac{\mathbb{E}(Q_W^r)}{\lambda (r-\mathbb{E}(Q_W^r))}.$$

$$\mathbb{E}(Q_W^r)$$

$$\sim \sqrt{\frac{\lambda r}{\mu }}{\int }_{\beta }^{\infty }(x-\beta )\hat{f}(x)dx$$

$${\rho }_{F^r}=1-O(1/\sqrt{r}),{\rho }_{B^r}=1-O(1/\sqrt{r}).$$

$$\mathbb{P}$$

$$={\left(1+\sqrt{\frac{\lambda }{\mu }}\frac{\varphi (\sqrt{\mu /\lambda }\gamma )}{\varphi (\gamma )}{e}^{\gamma (\beta -\gamma )}\frac{\Phi (\gamma )}{\Phi (-\sqrt{\mu /\lambda }\gamma )}+\sqrt{\frac{\lambda }{\mu }}\frac{\varphi (\sqrt{\mu /\lambda }\gamma )}{\gamma }\left(e^{\gamma (\beta -\gamma )}-1\right)\Phi {\left(-\sqrt{\frac{\mu }{\lambda }}\gamma \right)}^{-1}\right)}^{-1}.$$

$$\frac{\mathbb{E}(W)}{\sqrt{r}}\sim \frac{{\alpha }_3}{\sqrt{\lambda \mu }}\left(\sqrt{\frac{\mu }{\lambda }}\varphi \left(\frac{\mu }{\lambda }\gamma \right)\Phi {\left(-\sqrt{\frac{\mu }{\lambda }}\gamma \right)}^{-1}-\frac{\mu }{\lambda }\gamma \right),$$

Full text

1 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 [16]. 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 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 [12]. 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 [11, 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/\sqrt{r}$, and near-optimal resource utilization of $1-O(1/\sqrt{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 [5, 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 [16, 18, 17, 19, 15], 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. Load-balancing policies have attracted a lot of attention in recent years due to extremely relevant applications in large data centers, see [21] for an 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 stations 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 [6]. 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 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.

2 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\in\{1,\ldots,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 fully-charged 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 $\mu$, 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 $\lambda$. 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 define which station is relatively least loaded more precisely later in this section. Figure 1 illustrates the closed queueing model under this load-balancing arrival mechanism.

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

[TABLE]

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\in\{1,\ldots,S\}$. As a consequence, having more charging points creates no additional charging capacity, and can be bounded as

[TABLE]

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,\ldots,S$, with $G<\infty$ being a small fixed number.

The main quantity of interest in this paper is the number of batteries that are in need of charging, i.e. the aggregated number of batteries that are being charged at a charging point and the possible exchanged batteries that are waiting for an available charging point. We also refer this quantity as to as the queue length. Let $Q_{i}(t)$ denote the number of batteries in need of charging at station $i$ at time $t\geq 0$, and we write $Q(t)=(Q_{1}(t),\ldots,Q_{S}(t))$. 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.

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$. Under the policy where every arrival would choose randomly between the two stations in its direct vicinity, we observe that $p_{i}=\sum_{j=1}^{S}p_{ij}/2$ represents the effective arrival probability for every station $i=1,\ldots,S$. Therefore, for a system with $r$ cars, we set the capacity levels of the number of charging points and the number of spare batteries as

[TABLE]

for all $i=1,\ldots,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 $\gamma\leq\beta$ requirement in (2.4).

Since there are two types of resources at every station, i.e. charging points and spare batteries, one can consider two types of utilization levels. However, in view of (2.4), we see that the capacity levels of both resources are of the magnitude $p_{i}\lambda r/\mu+O(\sqrt{r})$, and hence the utilization levels of both resources are given by $Q(t)/(p_{i}\lambda r/\mu)(1+o(1))$. Using this observation, we define the relative occupancy level (load) of a station as $Q_{i}(t)/p_{i}$. We let our load-balancing policy prescribe that an EV in need of charging closest to station $i$ and $j$ at time $t\geq 0$ moves to station $i$ iff

[TABLE]

where ties are broken evenly. In our results, we show that this load-balancing policy ensures that the resource utilization levels at the different stations are approximately equal at all times. Consequently, this also ensures that the expected waiting times are approximately the same at every station at all times.

Remark 2.1**.**

Our modeling prescribes that every EV user can choose between two stations in its direct vicinity. We point out that this is done for simplicity, as it helps to describe our scaling and load balancing-policy in a clear and concise manner. We point out that our model and results extends naturally to the cases where some EV arrivals may always move to one station, and some EV arrivals choose from multiple stations. With respect to the modeling, this extension can be included as follows. Let $\mathcal{M}$ be the set of arrival types, where every $m\in\mathcal{M}$ is a set of stations that is in the direct vicinity of the EV user. Let $s_{m},m\in\mathcal{M}$ denote the probability that an EV arrival is of type $m$. Then, the effective arrival rate at any station $i\in\{1,...,S\}$ is given by

[TABLE]

In this extended setting, the scaling (2.4) for the number of resources and the load-balancing policy (2.5) remains the same.

3 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. An illustration of the closed queueing model is given in Figure 2. The square-root scaling rules reduces to

[TABLE]

where we suppress the subscript $1$ for the station number in this case.

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 different locations. 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 in 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.

At first glance, the single-station variant of our model may seem similar to the classic repair man model. This model and its QED-scaling implications are thoroughly treated in [10, 11], which mainly focus on the healthcare setting. We point out that there is a crucial difference: our single-station model includes spare batteries, causing none of $r$ cars to be waiting at the station as long as there are sufficient fully-charged spares available. If $B=0$, our model reduces to the repair man model with $r$ machines and $F$ repair men. Generally, however, the birth rates are different.

3.1 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 of whether the QED scaled provisioning rules (3.1) hold. More specifically, the queue length $\{Q(t),t\geq 0\}$ is a birth-death process with state space $Q(t)\in\{0,1,\ldots,B^{r}+r\}$ for all $t\geq 0$, with birth rate $\lambda(r-\left(Q(t)-B^{r}\right)^{+})$ and death rate $\mu\min\{Q(t),F^{r}\}$. Let

[TABLE]

denote the steady state distribution of the number of batteries in need of charging.

Lemma 3.1**.**

Suppose $S=1$, where the single swapping station has $F$ charging points and $B$ spare batteries, i.e. we disregard the scaling in(3.1). The steady state distribution is given by

[TABLE]

where

[TABLE]

if $F\leq B$, and

[TABLE]

*if $B\leq F$.

Remark 3.2**.**

In view of (2.1), we exclude the case that $F\geq B$ in our analysis further on in this paper. Yet, in an application where e.g. $G=\infty$ and hence $F\geq 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=\infty$ and $F=\infty$, Lemma 3.1 shows that

[TABLE]

where

[TABLE]

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.

3.2 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\geq 0$. We consider the fluid scaling

[TABLE]

The fluid-scaled process converges to a deterministic, continuous monotone process with a single fixed steady state value.

Proposition 3.3**.**

Suppose $S=1$ and scaling rules (3.1) hold. If $\bar{Q}^{r}(0)\rightarrow\bar{Q}(0)$ as $r\rightarrow\infty$ with $\bar{Q}(0)$ a finite constant, then $\bar{Q}^{r}\rightarrow\bar{Q}$ in distribution as $r\rightarrow\infty$, where $\bar{Q}$ satisfies the ODE

[TABLE]

and has the steady state value

[TABLE]

Proposition 3.3 implies that the number of batteries in need of charging can be approximated by

[TABLE]

where $\bar{Q}(t)=\lim_{r\rightarrow\infty}\bar{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})\approx\lambda r/\mu$ for some $t_{0}\geq 0$, then $Q^{r}(t)\approx\lambda r/\mu$ for all $t\geq 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 scaling

[TABLE]

This scaling provides more sensitive approximations, as it captures fluctuations of order $\sqrt{r}$. The diffusion-scaled 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). If $\hat{{Q}}^{r}(0)\rightarrow\hat{{Q}}(0)$ in distribution as $r\rightarrow\infty$, then $\hat{{Q}}^{r}\rightarrow\hat{{Q}}$ in distribution as $r\rightarrow\infty$. The process $\hat{{Q}}$ is a diffusion process with drift

[TABLE]

and constant infinitesimal variance $2\mu$. The steady state density of $\hat{Q}(\infty)=\lim_{t\rightarrow\infty}\hat{Q}(t)$ is given by

[TABLE]

where $\alpha_{i}=r_{i}/(r_{1}+r_{2}+r_{3})$, $i=1,2,3$ with

[TABLE]

Equation (3.18) in Theorem 3.4 is obtained by taking the limit of the scaled diffusion process (as $r\rightarrow\infty$), and finding its steady state distribution (as $t\rightarrow\infty$). 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\rightarrow\infty$) and next take the limit as $r\rightarrow\infty$. Fortunately, the following theorem shows that the order at which one takes the limit leads to the same result.

Theorem 3.5**.**

*If $S=1$ and (3.1) holds, the steady state distribution of the diffusion scaled process $\hat{Q}^{r}(\infty)$ converges in distribution to $\hat{Q}(\infty)$ as in Theorem 3.4.

The proof of Theorem 3.5 is given in Appendix A.2. As the order at which the limits are taken does not affect the result, we use the limiting process $\hat{Q}(\infty)$ to obtain approximations for the performance measures.

3.3 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\rightarrow\infty$, 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 fluctuations in arrival rate is of order $O(\sqrt{r})$, i.e. the arrival rate is $\lambda r-O(\sqrt{r})$ (with high probability). These small changes will therefore become negligible as $r\rightarrow\infty$. In other words, this argument implies that the PASTA property remains valid asymptotically. This notion can be formalized similarly as is done in [10]. Summarizing, if $W$ denotes the waiting time of an arriving EV user, then

[TABLE]

where $\hat{Q}(\infty)$ 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,

[TABLE]

where the throughput $\theta$ 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

[TABLE]

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

[TABLE]

In turn, the expected number of waiting cars can be derived directly using Theorem 3.4 and the observation $Q_{W}^{r}=(Q^{r}(\infty)-B^{r})^{+}$,

[TABLE]

as $r\rightarrow\infty$. We point out that $\mathbb{E}(Q_{W}^{r})$ is consequently of order $\Theta(\sqrt{r})$, and together with (3.19) this implies that $E(W)$ is of order $\Theta(1/\sqrt{r})$ and hence vanishes in the limit.

The resources will be fully utilized under (3.1) as $r\rightarrow\infty$. Theorem 3.4 implies that at most $O(\sqrt{r})$ charging points are not utilized, and the number of fully-charged batteries is also of order $O(\sqrt{r})$. Therefore, as $r\rightarrow\infty$,

[TABLE]

Theorem 3.6**.**

Suppose $S=1$, and the system is operating under (3.1). Then the following properties hold as $r\rightarrow\infty$. The waiting probability has a non-degenerate limit given by

[TABLE]

The expected waiting time behaves as

[TABLE]

with $\alpha_{i}$ are as in Theorem 3.4. Finally, the resource utilizations behave as

[TABLE]

The proof of Theorem 3.6 is given in Appendix A.3.

4 System behavior in case of multiple stations

When the number of stations $S\geq 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.5) 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 $\mathcal{S}$ of stations, there is at least one pair $(i,j)$ with $i\in\mathcal{S}$ and $j\not\in\mathcal{S}$ for which $p_{ij}>0$. This statement is equivalent to the following assumption.

Assumption 4.1**.**

*Let $G=(V,E)$ be a graph, where $V=\{1,\ldots,S\}$ and $E=\{(i,j):p_{ij}>0\}$. We assume that the graph $G$ is connected.

Remark 4.2**.**

For our results to follow through in the extended model as described in Remark 2.1, Assumption 4.1 needs to be updated as follows. Let $G=(V,E)$ be a graph, where $V=\{1,\ldots,S\}$ and $E=\{(i,j):i,j\in m,|m|\geq 2,m\in\mathcal{M}\}$. Then, we assume that the graph $G$ is connected. Note that if $m=2$ for every $m\in\mathcal{M}$, the setting as well as this assumption reduces to the original setting as described in this paper.

4.1 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 $\{\mathbb{X}^{r}(t),t\geq 0\}$ with

[TABLE]

where

• •

$A^{r}=\left(A_{ij}^{r};\{i,j\}\in E\right)$, where $A_{ij}^{r}(t)$ is the number of arrivals that are closest to stations $i$ and $j$ until time $t\geq 0$ in the $r$’th system;

• •

$A^{r}_{d}=\left(A_{ij,i}^{r};\{i,j\}\in E\right)$, where $A_{ij,i}^{r}(t)$ is the number of arrivals that are closest to stations $i$ and $j$ and are routed to station $i$ until time $t\geq 0$ in the $r$’th system;

• •

$Q^{r}=\left(Q_{j}^{r};1\leq j\leq S\right)$, where $Q_{j}^{r}(t)$ is the number of batteries in need of charging at time $t\geq 0$ in the $r$’th system;

• •

$Z^{r}=\left(Z_{j}^{r};1\leq j\leq S\right)$, where $Z_{j}^{r}(t)$ is the number of busy servers (charging points) at time $t\geq 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\geq 0$ in the $r$’th system;

• •

$T^{r}=\left(T_{j}^{r};1\leq j\leq S\right)$, where $T_{j}^{r}(t)$ is the aggregated time of all servers at station $j$ that were charging until time $t\geq 0$ in the $r$’th system;

• •

$D^{r}=\left(D_{j}^{r};1\leq j\leq S\right)$, where $D_{j}^{r}(t)$ is the number of service completions at station $j$ until time $t\geq 0$ in the $r$’th system;

• •

$L^{r}$, where $L^{r}(t)$ is the number of batteries that are positioned in an EV not waiting at a station in the $r$’th system at time $t\geq 0$.

Clearly, there are strong relations between the individual processes in $\mathbb{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 $\{\Lambda_{ij}(t),t\geq 0\}$ for all $\{i,j\}\in E$ be independent Poisson processes with rate $p_{ij}\lambda$ and $\{S_{j}(t),t\geq 1\}$ for all $1\leq j\leq S$ be independent Poisson processes with rate $\mu$. The system dynamics satisfy the following identities:

[TABLE]

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 [6]. 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 $\mathbb{D}^{d}[0,\infty]$ the $d$-dimensional Skorohod path space. For $x,y\in\mathbb{D}^{d}[0,\infty]$ and $T>0$, let

[TABLE]

where $|z|=\max_{i=1,\ldots,d}|z_{i}|$ for any $z=(z_{1},\ldots,z_{d})\in\mathbb{R}^{d}$. The space $\mathbb{D}^{d}[0,\infty]$ 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 $\{x_{n}\}\in\mathbb{D}^{d}[0,\infty]$ converges uniformly on compact sets (u.o.c) sets to $x\in\mathbb{D}^{d}[0,\infty]$ as $n\rightarrow\infty$ if for each $T\geq 0$,

[TABLE]

as $n\rightarrow\infty$. Moreover, we say that $t\geq 0$ is a regular point of a function $x$ if $x$ is differentiable at $t\geq 0$, and denote its derivative by $x^{\prime}(\cdot)$. We assume that the random variables in $\mathbb{X}^{r}$ live on the same probability space $(\Omega,\mathcal{F},\mathbb{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}(\cdot,\omega)$ for the sample path associated with $\omega\in\Omega$ for a stochastic process $X^{r}$.

4.2 Fluid limit

To capture the rough system dynamics, we consider the fluid-scaled process

[TABLE]

For each process $X^{r}$ in $\mathbb{X}^{r}$, we define similarly its fluid equivalent as $\bar{X}^{r}=X^{r}/r$ and its limiting process $\bar{X}=\lim_{r\rightarrow\infty}\bar{X}^{r}$. We adopt the definition of a fluid limit and its invariant state(s) from [9]. That is, we consider $\mathcal{A}\subset\Omega$ such that the FSLLN holds, i.e.

[TABLE]

u.o.c. as $r\rightarrow\infty$. Due to the FSLLN, we observe that one can choose $\mathcal{A}$ large enough such that $\mathbb{P}(\mathcal{A})=1$.

Definition 4.3**.**

*We call $\bar{\mathbb{X}}$ a fluid limit of $\{\mathbb{X}^{r}\}$ if there exists an $\omega\in\mathcal{A}$ and (sub)sequence $\{r_{l}\}$ with $r_{l}\rightarrow\infty$ as $l\rightarrow\infty$, such that $\bar{\mathbb{X}}^{r_{l}}(\cdot,\omega)$ converges u.o.c. to $\bar{\mathbb{X}}(\cdot,\omega)$. Moreover, let ${q}=(q_{1},\ldots,q_{S})$ be an invariant state of the fluid limits if for any fluid limit $\bar{\mathbb{X}}$, $\bar{Q}(0)=(\bar{Q}_{1}(0),\ldots,\bar{Q}_{S}(0))=(q_{1},\ldots,q_{S})=q$ implies that $\bar{Q}(t)=q$ for all $t\geq 0$. *

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 requiring $\bar{Q}^{r}(0)\rightarrow\bar{Q}(0)$ with $\bar{Q}(0)$ a finite constant, Definition 4.3 allows for $\bar{Q}(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 of $\bar{Q}(0)$. Moreover, there is a single unique invariant state given by $\lambda/\mu$. A similar result holds when $S\geq 2$.

Theorem 4.4**.**

Let $\{\mathbb{X}^{r}\}$ be a sequence of systems. Then the fluid limits exist, where each component is Lipschitz continuous. Each fluid limit $\bar{\mathbb{X}}$ satisfies the following equations for all $t\geq 0$:

[TABLE]

Also, for every $\{i,j\}\in E$, if $t$ is a regular point of $\bar{\mathbb{X}}$, then

[TABLE]

*Finally, there is a unique invariant state given by $q=(q_{1},\ldots,q_{S})$ with $q_{i}=p_{i}\lambda/\mu$ for $i=1,\ldots,S$.

The (uniqueness of the) invariant state result for the fluid limit is central for the existence of a properly-defined diffusion process as it states that if $\bar{Q}(0)=(p_{1}\lambda/\mu,\ldots,p_{S}\lambda/\mu)$, the fluid limits are time invariant. We present a proof of Theorem 4.4 in Appendix D.

4.3 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.4), 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 as

[TABLE]

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.

Theorem 4.5**.**

Suppose

[TABLE]

as $r\rightarrow\infty$, where $\hat{Q}(0)$ is a random vector. Then, for every $K^{r}=o(\sqrt{r})$ with $K^{r}\rightarrow\infty$ as $r\rightarrow\infty$, and for every $T>0$ and $\epsilon>0$,

[TABLE]

for every $i,j\in\{1,\ldots,S\}$ as $r\rightarrow\infty$. If, in addition for every $i,j\in\{1,\ldots,S\}$,

[TABLE]

then

[TABLE]

*for every $i,j\in\{1,\ldots,S\}$ as $r\rightarrow\infty$.

The proof of Theorem 4.5 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}_{\Sigma}(t),t\geq 0\}$ with $r\in\mathbb{N}$, where $Q^{r}_{\Sigma}(t)=\sum_{j=1}^{S}Q_{j}^{r}(t)$, and

[TABLE]

As the state space collapse implies that $\hat{Q}_{i}^{r}(t)\approx\hat{Q}_{j}^{r}(t)$ for all $i,j\in\{1,\ldots,S\}$ (for $t\geq K^{r}/\sqrt{r}$), we can approximate the queue length at an individual queue by

[TABLE]

for all $j=1,\ldots,S$. The limiting process of the total queue length can be derived using the SSC result.

Theorem 4.6**.**

Suppose $\hat{Q}^{r}(0)\rightarrow\hat{Q}(0)$ in distribution as $r\rightarrow\infty$, and

[TABLE]

for all $i,j\in\{1,\ldots,S\}$. Then, $\hat{Q}^{r}_{\Sigma}\rightarrow\hat{Q}_{\Sigma}$ in distribution as $r\rightarrow\infty$, where $\hat{Q}_{\Sigma}$ is a diffusion process with drift

[TABLE]

and constant infinitesimal variance $2\mu$. The steady state density $\hat{Q}_{\Sigma}(\infty)$ is given by

[TABLE]

where $\alpha_{i}=r_{i}/(r_{1}+r_{2}+r_{3})$, $i=1,2,3$ with

[TABLE]

Proof.

We observe that the steady state density is a direct consequence of the diffusion process [7]. What remains to show is that $\hat{Q}_{\Sigma}^{r}$ converges to the described diffusion process as $r\rightarrow\infty$. Equivalently, we need to show that

[TABLE]

where $\{W(t),t\geq 0\}$ is a standard Brownian motion. We note that due to the system identities,

[TABLE]

where

[TABLE]

and

[TABLE]

We observe that due to the FCLT and Theorem 4.4,

[TABLE]

where $\{\textrm{BM}_{A}(t),t\geq 0\}$ is Brownian motion with mean zero and variance $\mu$. Similarly, due to the FCLT and Theorem 4.4,

[TABLE]

where $\{\textrm{BM}_{D}(t),t\geq 0\}$ is an (independent) Brownian motion with mean zero and variance $\mu$. The sum of these two processes results is equal to (in distribution) a Brownian with mean zero and variance $2\mu$, which contributes to the $\sqrt{2\mu}\,dW(t)$ term in (4.25). Next, we observe that due to the system identities and the definition of the diffusion scaling,

[TABLE]

Since

[TABLE]

for all $t\in[0,T]$, and due to Theorem 4.5,

[TABLE]

where the sequence $\{\epsilon(r),r\in\mathbb{N}\}$ can be chosen such that $\epsilon(r)\rightarrow 0$ as $r\rightarrow\infty$. Then, as $r\rightarrow\infty$,

[TABLE]

This contributes to the first two terms in (4.25). 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.7**.**

Suppose the system is operating under (3.1). Then the following properties hold as $r\rightarrow\infty$ for all $i=1,\ldots,S$. The waiting probability has a non-degenerate limit given by

[TABLE]

The expected waiting time behaves as

[TABLE]

with $\alpha_{i}$ are as in Theorem 4.6. Finally, the resource utilizations behave as

[TABLE]

5 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. We first focus on a large-scale system to illustrate the implications of our results. Next, we also zoom in on a moderate-sized system that reflect a more realistic setting for an EV battery swapping infrastructure.

5.1 Large-scaled system

Throughout the experiments in this section, we consider a network with $5$ stations where the arrival probabilities are given by $p_{12}=p_{24}=p_{34}=p_{35}=0.1$, $p_{23}=0.4$ and $p_{45}=0.2$, see Figure 3. This results in an effective arrival probability $p=(p_{1},\ldots,p_{5})=(0.05,0.3,0.3,0.2,0.15)$ at the stations.

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 following assumptions. We assume that the battery swapping facility installed (relatively) fast charging points where recharging takes one hour on average ($\mu=1$), and that every EV user returns for recharging services after every $40$ hours on average ($\lambda=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. The effective loads at the stations in this case are

[TABLE]

and we note that due to the QED provisioning rule in (2.4), 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

[TABLE]

5.1.1 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.4. We validate this approximation for the above-described setting, with initial queue length $Q_{i}^{r}(0)=150$ for all stations $i=1,\ldots,5$. That is, only station 1 is initially overloaded, while all other station are underloaded. The equations in Theorem 4.4 together with the Lipschitz continuity describe a unique fluid limit with the given initial queue length. This yields the approximations

[TABLE]

In particular, in the case when the initial queue length is $Q_{i}^{r}(0)=150$ for all stations $i=1,\ldots,5$, this results in the approximations

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where $t_{1}\approx 0.1826$, $t_{2}\approx 0.3189$, $t_{3}=1.4$ and $t_{4}\approx 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}\approx 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.

5.1.2 Performance measures

Our results imply approximations for performance measures such as the waiting probability and waiting time, see Corollary 4.7. In particular, the state-space collapse result implies that the performance at all stations is approximately the same, and can be approximated by the closed-form expressions as given in Corollary 4.7.

In Figure 6, we plotted the waiting probabilities of all stations in the case of $2500000$ EV arrivals averaged over $20$ samples for the large-scaled system. We point out that the stair-type 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. This is also reflected in Figure 5, where a closer view suggests that the diffusion-scaled queue length at station 1 being smaller than the queue length at another station (often station 2 or station 3). 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.4) guarantee that a desired waiting probability is achieved.

5.2 Universality result for charging 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 7 and 8, we consider the system setting as described in Section 5.1. 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\geq 1$. Similarly to the setting with exponential charging times, the maximum difference between the diffusion-scaled queue length behaves quite erratically, see Figure 9. 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.

5.3 The role of system size

In the previous sections, we commented at a view points that 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 diffusion-scaled 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}\sqrt{\lambda r/\mu})\approx 0.5657$ in the described setting in Section 5.1. Theorem 4.5 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 10, which is monotonically decreasing in $r$. In addition, we observe that the average maximum distance, i.e. $1/T\int_{0}^{T}\max_{1\leq i<j\leq S}\{|Q_{i}(t)-Q_{j}(t)|\,\}dt$, is also monotonically decreasing in $r$, and is not excessively smaller than the maximum distance of the interval.

Summarizing, as the system size increases, the accuracy of the approximations improve. However, one can imagine that a real-life battery swapping infrastructure is not of the scale as discussed in the previous sections. Therefore, we consider how our results hold up in a more realistic for an EV battery swapping infrastructure.

5.4 Moderate-sized system

We consider the following setting. We have the same network structure as given in Figure 3, but with other arrival probabilities. More specifically, we assume $p_{12}=p_{23}=p_{24}=p_{25}=0.2$ and $p_{13}=p_{34}=0.1$, giving an effective arrival probability of $(p_{1},\ldots,p_{5})=(0.1,0.25,0.3,0.2,0.15)$. For the other parameters, we assume that recharging takes four hours on average ($\mu=0.25$) and every EV user returns for recharging services after every 50 hours on average ($\lambda=0.02$). We assume that there our infrastructure consists of a thousand electric vehicles ($r=1000$). The effective load at the stations are

[TABLE]

Note that $\lambda r/\mu=80$ and $\sqrt{\lambda r/\mu}=4\sqrt{5}\approx 8.94$. Relatively, their sizes are much closer to one another than when $r$ becomes larger, and this will also have its impact on the behavior.

Due to the smaller system size, we can expect that the fluctuations of around its fluid limit is relatively larger. Indeed, if one plots a sample path for this system with initial queue length $Q(0)=(0,0,0,0,0)$, we observe that the fluctuations to its corresponding fluid limit approximation are significant, see Figure 11 and 12. Moreover, in Figure 12, the queue lengths even seem to lie above the fluid limit results. We point out that this is a consequence of the high waiting probabilities for these parameter settings. Moreover, since the fluctuations are of a significant size (relatively), this leads to sample paths that appear quite far off (above) its fluid limit approximation at first glance.

To see whether a load-balancing effect still takes place for a moderate-sized system, we should consider the diffusion scaled queue lengths, see Figure 13. Numerous experiments, including the setting as in Figure 13 suggest that even for these settings, the effect of state-space collapse is very much visible. In other words, there is still a strong load-balancing effect present that leads to the occupation level at the different stations stay close to one another. In turn, this lead that performance, e.g. waiting probability and waiting times, are comparable at the different stations at all times.

Acknowledgment

This research is supported by the Netherlands Organisation for Scientific Research through the programmes Gravitation NETWORKS grant 024.002.003, MEERVOUD grant 632.003.002, and VICI grant 639.033.413. The work of Seva Shneer is supported by the Mathematical Center in Akademgorodok under agreement No. 075-15-2019-1675 with the Ministry of Science and Higher Education of the Russian Federation.

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

Proof of Proposition 3.3.

We observe that the $Q^{r}(\cdot)$ is a birth-death process with state space $\mathcal{Q}^{r}=\{0,1,\ldots,B^{r}+r\}$, arrival (birth) rate $\lambda_{r}(j)=\lambda\left(r-(j-B^{r})^{+}\right)$ and service (death) rate $\mu_{r}(j)=\mu\min\{j,F^{r}\}$ for all $j\in\mathcal{Q}^{r}$. The fluid-scaled process $\bar{Q}^{r}$ therefore has state space $\{0,1/r,\ldots,(B^{r}+r)/r\}$, and drift and infinitesimal variance functions

[TABLE]

and

[TABLE]

Taking the limit yields under scaling rule (3.1)

[TABLE]

and

[TABLE]

That is, we obtain a degenerate limiting diffusion, and hence the limiting initial point will yield a deterministic path for the limiting process $\bar{Q}$. The mean directly provides the ODE that the limiting process $Q$ satisfies. The limiting value (as $t\rightarrow\infty)$ of the limiting process is also correct due to the following reasoning. For all $t\geq 0$ for which $\bar{Q}(t)<\frac{\lambda}{\mu}$, we observe $d\bar{Q}(t)/dt>0$. Hence, if $\bar{Q}(0)<\lambda/\mu$, then the limiting process follows by the deterministic path

[TABLE]

Similarly, for all $t\geq 0$, if $\bar{Q}(t)>\lambda/\mu$, then $d\bar{Q}(t)/dt<0$. Hence if $\bar{Q}(0)>\lambda/\mu$ , then the limiting process follows the deterministic path

[TABLE]

Finally, if $\bar{Q}(0)=\lambda/\mu$, then $d\bar{Q}(t)/dt=0$ for all $t\geq 0$, which yields $\bar{Q}(t)=\lambda/\mu$ for all $t\geq 0$.

The second result involves the diffusion limit.

Proof of Theorem 3.4.

The infinitesimal mean of the centered scaled process is given by

[TABLE]

and the infinitesimal variance is given by

[TABLE]

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.18), where $\alpha_{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

Theorem 3.5 states that it is allowed to interchange the limits. To prove this result, we need to show that the steady state distribution (as $t\rightarrow\infty$) converges to the same distribution as $r\rightarrow\infty$ as the diffusion limiting result. First we derive the limiting steady state distribution conditioned to be in one of the intervals $[0,F^{r}]$, $[F^{r},B^{r}]$, and $[B^{r},B^{r}+r]$.

Lemma A.1**.**

Suppose $S=1$. Under scaling rules (3.1), the following limiting results hold as $r\rightarrow\infty$:

[TABLE]

Proof.

The first expression follows from the central limit theorem and the properties of the Poisson distribution,

[TABLE]

The second expression can be obtained using geometric series,

[TABLE]

The final expression again uses the central limit theorem and the properties of the Poisson distribution. That is,

[TABLE]

Since for every $z\in\mathbb{R}$,

[TABLE]

as $r\rightarrow\infty$, we obtain

[TABLE]

To obtain the limiting distribution of the steady state distribution as $r\rightarrow\infty$, 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.

Lemma A.2**.**

Under scaling rules (3.1), as $r\rightarrow\infty$,

[TABLE]

Proof.

This follows from the central limit and the properties of the Poisson distribution,

[TABLE]

Lemma A.3**.**

Under scaling rules (3.1) with $\gamma\neq 0$, as $r\rightarrow\infty$,

[TABLE]

If $\gamma=0$, then

[TABLE]

Proof.

If $\gamma\neq 0$, we observe that

[TABLE]

due to geometric series. Moreover,

[TABLE]

and due to geometric series,

[TABLE]

In addition,

[TABLE]

Finally, using Stirling’s approximation

[TABLE]

Combining the expressions yields

[TABLE]

If $\gamma=0$, using Stirling’s approximation, we directly obtain

[TABLE]

Lemma A.4**.**

Under scaling rules (3.1), as $r\rightarrow\infty$,

[TABLE]

Proof.

We observe that

[TABLE]

Moreover, applying Stirling’s approximation twice,

[TABLE]

Multiplying the two expression concludes the result.

Lemma A.5**.**

Under scaling rules (3.1), as $r\rightarrow\infty$:

[TABLE]

where

[TABLE]

Proof.

It follows from Lemmas A.2-A.4 that

[TABLE]

and applying Lemmas A.2-A.4 again yields

[TABLE]

Consequently, it follows that Theorem 3.5 holds.

Proof of Theorem 3.5.

This result follows directly from Lemmas A.1 and A.5.

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,

[TABLE]

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

[TABLE]

where $\alpha_{3}$ is as in Theorem 3.4. We remark that this expression is always positive: trivially for $\gamma\leq 0$, and also for $\gamma>0$ since for every $x>0$

[TABLE]

which can be shown by partial integration. Due to Little’s law, we note that

[TABLE]

which yields the result for the expected waiting time.

Finally, for the utilization levels, the scaling of $\hat{Q}(\infty)$ implies that

[TABLE]

Therefore, the utilization levels satisfy (3.20) and the result follows.

Appendix B Unlimited number of charging points for single station system

A related setting involves the case where there are an unlimited number of charging points ($F=\infty$) and also an unlimited number of swapping servers $G=\infty$. Comparing these results to the performance of a related system where $G<\infty$ and the charging points scaled as in (3.1) provides good insight into the significance of the effect of these resources. The number of spare batteries remains as before, i.e.

[TABLE]

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 $\{\bar{Q}^{r}(t),t\geq 0\}$, where $\bar{Q}^{r}(t)=Q^{r}(t)/r$ for all $r\geq 1$. The following theorem holds under the QED scaling.

Theorem B.1**.**

If $S=1$ and $\bar{Q}^{r}(0)\rightarrow\bar{Q}(0)$ in distribution with $\bar{Q}(0)$ a finite constant, then $\hat{Q}^{r}\rightarrow\bar{Q}$ in distribution as $r\rightarrow\infty$, where $\bar{Q}$ satisfies the ODE

[TABLE]

and has the unique steady state value

[TABLE]

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 $\mathcal{Q}^{r}=\{0,1,\ldots,B^{r}+r\}$, and arrival rates $\lambda_{r}(j)=\lambda\left(r-(j-B^{r})^{+}\right)$ and service rate $\mu_{r}(j)=j\mu$ for all $j\in\mathcal{Q}^{r}$. The fluid-scaled process $\bar{Q}^{r}$ therefore has state space $\{0,1/r,\ldots,(B^{r}+r)/r\}$, and drift and diffusion functions

[TABLE]

and

[TABLE]

Taking the limit yields

[TABLE]

and

[TABLE]

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 process $\bar{Q}$. The mean directly provides the ODE that the limiting process $\bar{Q}$ satisfies. What remains to show is the limiting value (as $t\rightarrow\infty)$ .

For all $t\geq 0$ for which $\bar{Q}(t)<\frac{\lambda}{\mu}$, we have

[TABLE]

and for all $t\geq 0$ for which $\bar{Q}(t)>\frac{\lambda}{\mu}$, we obtain

[TABLE]

Moreover, if $\bar{Q}(t)=\lambda/\mu$ for some $t\geq 0$, then $d\bar{Q}(t)/dt=0$. This implies that the initial limiting value $\bar{Q}(0)$ determines the unique deterministic path for the limiting process. That is, if $\bar{Q}(0)<\lambda/\mu$, then the limiting process satisfies

[TABLE]

a strictly increasing function that tends to $\lambda/\mu$ as $t\rightarrow\infty$. Similarly, if $\bar{Q}(0)>\lambda/\mu$, then the limiting process satisfies

[TABLE]

a strictly decreasing function that tends to $\lambda/\mu$ as $t\rightarrow\infty$. Finally, if $\bar{Q}(0)=\lambda/\mu$, then $d\bar{Q}(t)/dt=0$ and hence $\bar{Q}(t)=\lambda/\mu$ for all $t\geq 0$.

We now zoom in on the fluctuations around the fluid limit by taking the diffusion limit. We consider the centered scaled process

[TABLE]

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). If $\hat{Q}^{r}(0)\rightarrow\hat{Q}(0)$ in distribution as $r\rightarrow\infty$, then $\hat{Q}^{r}\rightarrow\hat{Q}$ in distribution as $r\rightarrow\infty$, where $\hat{Q}$ can be described as a two pieced-together Ornstein-Uhlenbeck (OU) process. More precisely, $\hat{Q}$ is a diffusion process with mean

[TABLE]

and constant infinitesimal variance $2\mu$. The steady state density $\hat{Q}(\infty)$ is given by

[TABLE]

where

[TABLE]

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

[TABLE]

which converges to

[TABLE]

as $r\rightarrow\infty$. The infinitesimal variance tends to

[TABLE]

as $r\rightarrow\infty$. 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.5), where $\alpha$ is the probability that the steady state process is above $\beta$. Equations (5) and (6) from [7] yield

[TABLE]

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\rightarrow\infty$ with respect to the limiting diffusion process. In other words, the steady state density function (B.5) is obtained by first taking the limit of the centered scaled process as $n\rightarrow\infty$, after which its steady state distribution is obtained by setting $t\rightarrow\infty$. The following result shows that interchanging the limits does not change the result.

Theorem B.3**.**

As $r\rightarrow\infty$,

[TABLE]

*in distribution, where $\hat{{Q}}(\infty)$ is the steady state distribution as in Theorem B.2.

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.12).

Lemma B.4**.**

Under scaling rule (B.1), as $r\rightarrow\infty$,

[TABLE]

Proof.

Note that

[TABLE]

The convergence of the first term follows directly from the CLT for a Poisson distributed random variable,

[TABLE]

We extract from the second term

[TABLE]

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

[TABLE]

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

[TABLE]

*as $r\rightarrow\infty$.

Proof.

Note that

[TABLE]

It follows from the proof of Lemma B.4 that, as $r\rightarrow\infty$,

[TABLE]

Rewriting the terms concludes the result.

Next, we prove Theorem B.3.

Proof of Theorem B.3.

By definition,

[TABLE]

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.11). Due to Proposition B.5, it suffices to show that for every $x<\beta$

[TABLE]

and for every $x\geq\beta$,

[TABLE]

Using the steady state distribution given in (3.11), we observe that for every $\beta$,

[TABLE]

as $r\rightarrow\infty$. For every $x\geq\beta$,

[TABLE]

Similarly to the proof of Proposition B.5, The CLT implies

[TABLE]

and

[TABLE]

as $r\rightarrow\infty$.

B.3 Performance measures

We consider the waiting probability, the waiting time and the utilization level of the spare batteries. We note that $\alpha$ in (B.5) corresponds to the probability that $\{\hat{Q}(\infty)\geq\beta\}$. 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\rightarrow\infty$,

[TABLE]

For the waiting time, we point out that in this case the identity

[TABLE]

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\rightarrow\infty$,

[TABLE]

Consequently, the waiting time of a car behaves as

[TABLE]

*as $r\rightarrow\infty$.

Proof.

We note that as $r\rightarrow\infty$,

[TABLE]

Elementary calculus yields

[TABLE]

Then as $r\rightarrow\infty$, we obtain

[TABLE]

where

[TABLE]

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 $\mathbb{E}(\textrm{\# Fully-Charged Batteries})=\Theta(\sqrt{r})$ and $\mathbb{E}(Q^{r})=\Theta(r)$.

Corollary B.8**.**

Under scaling rule (B.1),

[TABLE]

Appendix 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=\infty$ or equivalently $G\geq\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}\leq B^{r}$, this reduces to the main setting in this paper. Therefore, in this section, we consider only

[TABLE]

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 $\Theta(\sqrt{r})$ are too small to be of any impact on the fluid scale. There will however be a notable effect on the diffusion scale.

Theorem C.1**.**

Suppose the system operates under the QED scaling rule (3.1). If $\hat{Q}^{r}(0)\rightarrow\hat{Q}(0)$ in distribution as $r\rightarrow\infty$, then $\hat{Q}^{r}\rightarrow\hat{Q}$ in distribution as $r\rightarrow\infty$. The process $\{\hat{Q}(t),t\geq 0\}$ is a diffusion process with mean

[TABLE]

and constant infinitesimal variance $2\mu$. The steady state density $\hat{Q}(\infty)$ is given by

[TABLE]

where $\alpha_{i}=r_{i}/(r_{1}+r_{2}+r_{3})$, $i=1,2,3$ with

[TABLE]

Proof.

The infinitesimal mean of the centered scaled process is given by

[TABLE]

and the infinitesimal variance is given by

[TABLE]

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.6), where $\alpha_{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\rightarrow\infty$,

[TABLE]

Proof.

The first expression follows, similarly to the case of $F^{r}<B^{r}$, from the CLT and the properties of the Poisson distribution:

[TABLE]

The second expression follows as in the proof in Proposition 14, using the binomial distribution and the CLT:

[TABLE]

The final expression again uses the CLT and the properties of the Poisson distribution:

[TABLE]

Since for every $z\in\mathbb{R}$,

[TABLE]

as $r\rightarrow\infty$, the result follows immediately.

Lemma C.3**.**

Under scaling rules (3.1), as $r\rightarrow\infty$,

[TABLE]

where

[TABLE]

Proof.

Applying Lemmas C.4-C.6 below yields the result.

Lemma C.4**.**

Under scaling rules (3.1), as $r\rightarrow\infty$,

[TABLE]

Proof.

This follows from the CLT and the properties of the Poisson distribution.

Lemma C.5**.**

Under scaling rules (3.1), as $r\rightarrow\infty$,

[TABLE]

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\rightarrow\infty$,

[TABLE]

Proof.

We observe that

[TABLE]

Moreover, applying Stirling’s approximation twice,

[TABLE]

Multiplying the two expressions concludes the result.

C.3 Performance measures

Theorem C.7**.**

Under scaling rules (3.1),

[TABLE]

The expected waiting time behaves as

[TABLE]

with $\alpha_{i}$ as in Theorem C.1. Finally, the utilization levels behave as

[TABLE]

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,

[TABLE]

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

[TABLE]

Due to Little’s law, we have the identity

[TABLE]

which yield the results for the expected waiting time.

Finally, for the utilization levels, the scaling of $\hat{Q}(\infty)$ implies that

[TABLE]

Therefore, the utilization levels satisfy (3.20) and the result follows.

Appendix D Fluid limit proofs

The proof of Theorem 4.4 is similar to the proof of [9], but adapted appropriately to our system.

Proof of Theorem 4.4.

First, we show that the fluid limits exist, where all components are Lipschitz continuous. For this purpose, we show that for all $\omega\in\mathcal{A}$ the sequence $\{\bar{\mathbb{X}}(\cdot,\omega)\}$ has a convergent subsequence, where every component in the limiting process is Lipschitz continuous.

Fix $\omega\in\mathcal{A}$. We observe that for every $0\leq t_{1}<t_{2}$,

[TABLE]

Therefore, there exists a subsequence $\{r_{l}\}$ such that $\bar{T}_{j}^{r_{l}}(\cdot,\omega)$ converges u.o.c. to some $\bar{T}_{j}(\cdot,\omega)$ as $l\rightarrow\infty$ for every $j=1,\ldots,S$, which is Lipschitz continuous.

Using Lemma 11 in [1], Equation (4.4) and the fact that $\omega\in\mathcal{A}$, it follows that $\bar{D}^{r}_{j}(\cdot,\omega)$ also converges u.o.c. to $D_{j}(\cdot,\omega)$ for every $j=1,\ldots,S$. In fact, it follows that $D_{j}(\cdot,\omega)=\mu\bar{T}_{j}(\cdot,\omega)$, and is therefore also Lipschitz continuous.

Next, we consider the arrival processes. First, we observe that $L^{r}(t)\leq r$ for every $t\geq 0$. Therefore,

[TABLE]

for all $0\leq t_{1}<t_{2}$, and hence there exists a subsequence $\{r_{l}\}$ such that $\bar{Y}^{r_{l}}(\cdot,\omega)$ converges u.o.c. to some $\bar{Y}(\cdot,\omega)$ as $l\rightarrow\infty$ for every $j=1,\ldots,S$, which is again Lipschitz continuous.

Moreover, it follows from (4.2) for all $0\leq t_{1}<t_{2}$,

[TABLE]

As $\omega\in\mathcal{A}$ and the FSLLN applies, it follows from Theorem 12.3 in [3] that there is some subsequence $\{r_{l}\}$ such that $\bar{A}_{ij}^{r}(\cdot,\omega)$ converges u.o.c. as $l\rightarrow\infty$ to some process $\bar{A}_{ij}(\cdot,\omega)$. In particular, it holds for all $0\leq t_{1}<t_{2}$ that

[TABLE]

and $\bar{A}_{ij}$ is hence Lipschitz continuous. Similarly, we can show the same convergence result for the processes $\bar{A}_{ij,i}^{r}(\cdot,\omega)$ to $\bar{A}_{ij,j}(\cdot,\omega)$.

By (4.3), it follows also that $\{Q_{j}^{r_{l}}(\cdot,\omega)\}$ is precompact, which in turn implies that $\{Z_{j}^{r_{l}}(\cdot,\omega)\}$ is precompact due to (4.7). Moreover, $\{L^{r_{l}}(\cdot,\omega)\}$ is precompact by (4.8). In conclusion, the fluid limit exists with each component being 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. Suppose $\bar{\mathbb{X}}$ to be a fluid limit with corresponding $\omega\in\mathcal{A}$ and subsequence $\{r_{l}\}_{l\in\mathbb{N}}$. If for some $t>0$ we have that $\bar{Q}_{j}(t)/p_{j}>\bar{Q}_{i}(t)/p_{i}$, then it follows by the continuity of the fluid limit that there exists a $\delta>0$ such that

[TABLE]

for all $s\in[t-\delta,t+\delta]$. By definition of the fluid limit, it holds for large enough $r_{l}$ that

[TABLE]

for all $s\in[t-\delta,t+\delta]$. 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-\delta,t+\delta]$, and hence $\bar{A}^{\prime}_{ij,j}(t)=0$. Moreover, station $i$ receives all arrivals and by the FSLLN and (4.11),

[TABLE]

for all $t_{1}<t_{2}$ with $t_{1},t_{2}\in[t-\delta,t+\delta]$. It follows that $\bar{A}^{\prime}_{ij,i}(t,\omega)=p_{ij}\lambda\bar{L}(t,\omega)$ by (4.1).

Finally, we show that there is a unique invariant state given by $q=(p_{1}\lambda/\mu,\ldots,p_{S}\lambda/\mu)$. Introduce the function

[TABLE]

and write

[TABLE]

Trivially, $h(t)\geq 0$. Since $\bar{Q}(\cdot)$ is Lipschitz continuous, so is $h(\cdot)$, and hence it is differentiable almost everywhere. To show that $h(0)=0$ implies $h(t)=0$ for all $t\geq 0$, it therefore suffices to show that if $h(t)>0$ then $h^{\prime}(t)<0$ for every regular point $t$ of $\bar{\mathbb{X}}$. By (4.12), we observe that for every $j=1,\ldots,S$ and regular point $t\geq 0$,

[TABLE]

Due to (4.13)-(4.16), $\bar{D}^{\prime}_{j}(t)=\min\{\mu\bar{Q}_{j}(t),p_{j}\lambda\}$. In particular, we observe that $\bar{D}^{\prime}_{i}(t)/p_{i}=\bar{D}^{\prime}_{j}(t)/p_{j}$ for all $i,j\in\bar{S}_{\max}(t)$, as well as for all $i,j\in\bar{S}_{\min}(t)$. Due to Lemma 2.8.6 from [8], as $t$ is a regular point, it follows that

[TABLE]

for every $j,l\in\bar{S}_{\max}(t)$, as well as for every $j,l\in\bar{S}_{\min}(t)$. Due to (4.10), (4.11) and (4.18), we conclude that for every $j\in S_{\max}(t)$,

[TABLE]

On the other hand, for every $j\in S_{\min}(t)$,

[TABLE]

Observing that $\bar{D}^{\prime}_{j}(t)/p_{j}>\bar{D}^{\prime}_{i}(t)/p_{i}$ for every $j\in S_{\max}(t)$ and $i\in S_{\min}(t)$, we therefore conclude that if $h(t)>0$ with $t$ a regular point, then

[TABLE]

In other words, for every invariant state of the fluid limit, it must hold that $\bar{Q}_{i}/p_{i}(t)=\bar{Q}_{j}/p_{j}(t)$ for every $i\neq j$. We observe, in view of (4.12), that

[TABLE]

where for every $1\leq j\leq S$,

[TABLE]

and hence in view of (4.10) and (4.11),

[TABLE]

That is, $\bar{Q}^{\prime}_{j}(t)=0$ if and only if

[TABLE]

In view of (4.17), this occurs if and only if $\bar{Q}_{j}(t)=p_{j}\lambda/\mu$ for every $j=1,\ldots,S$.

Appendix E State space collapse proofs

To prove Theorem 4.5, we use a framework similar to that of [6], and [9]. The construction consists of several steps, which we lay out next.

Divide interval $[0,T]$ into $T\sqrt{r}$ intervals of length $1/\sqrt{r}$, indexed by $m$. In each interval, consider the hydrodynamically-scaled process of $\mathbb{X}$. For each of these interval, we

1. (a)

show the scaled process is ”almost” Lipschitz continuous; 2. (b)

show convergence to some hydrodynamic limiting process for a sufficiently large part of the state space; 3. (c)

derive the hydrodynamic limit equations. 2. 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. 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\sqrt{r}$ intervals of length $1/\sqrt{r}$, indexed by $m$. We write $p=(p_{1},\ldots,p_{S})$, and

[TABLE]

For the processes in $\mathbb{X}$, we introduce the following hydronimically-scaled variants. For $Q^{r}$, $Z^{r}$ and $L^{r}$, let

[TABLE]

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

[TABLE]

In other words, we track the increase of these processes during the interval $[m/\sqrt{r},m/\sqrt{r}+\sqrt{x_{r,m}}t/r]$. By definition of $x_{r,m}$, we note that

[TABLE]

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 $\sqrt{x_{r,m}}/r$ is very small for all $\omega\in\mathcal{A}$.

Lemma E.1**.**

Suppose $\hat{Q}^{r}(0)\rightarrow\hat{Q}(0)$ for some random vector $\hat{Q}(0)$, and let $M>0$ fixed. For every $\epsilon>0$ and $\omega\in\mathcal{A}$,

[TABLE]

*for $r$ large enough.

Proof.

Due to our fluid limit result in Theorem 4.4 and the definition of $x_{r,m}$ in (E.1), we observe that

[TABLE]

for $r$ large enough. Moreover, for $r$ large enough,

[TABLE]

We conclude that for every $m\leq\sqrt{r}T$,

[TABLE]

and

[TABLE]

For the hydrodynamically scaled process, the system identities translate to

[TABLE]

In order to show that $\mathbb{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 $\epsilon>0$, $M>0$ and $T>0$. For $r$ large enough, there exists a constant $N>0$ (only depending on $\lambda$, $\mu$, and $\{p_{ij};\{i,j\}\in E\}$) such that

[TABLE]

Moreover,

[TABLE]

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\}\in E$,

[TABLE]

From (E.14) and (E.17), we observe that $Y_{r,m}(t)\leq t$ and non-decreasing. Due to the properties of Poisson processes,

[TABLE]

Therefore,

[TABLE]

where the final to last inequality follows from Proposition 4.3 in [6]. Choosing $N=1/(\lambda\min_{\{i,j\}\in E}p_{ij})$ and applying the union bound twice, we obtain

[TABLE]

which yields (E.19).

The proof for (E.20) is completely analogous, but with minor adaptions as one uses $T_{j}^{r,m}(t)\leq p_{j}\lambda/\mu t$ instead of $Y^{r,m}(t)\leq 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\geq 1\}$ be independent exponentially distributed random variables with rate $p_{ij}\lambda$, representing the time that a car has before it needs recharging at either station $i$ or $j$. Let $\{v_{j}(l),l\geq 1\}$ be independent exponentially distributed random variables with rate $\mu$, representing the service requirement (recharging time) of a battery at station $j$. Define

[TABLE]

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

[TABLE]

Moreover, due to (4.2) and (4.4), we observe

[TABLE]

As in [9], we define for notational convention, for $b=(b_{1},b_{2})\in\mathbb{N}$,

[TABLE]

In view of (E.23) and (E.24), this yields the inequalities

[TABLE]

Using these processes, we first prove the following bounds:

[TABLE]

for $b=(1,0)$ and $b=(0,1)$. The proof is similar to that of (78) in [20]. We observe that the proof of (E.19) implies that in particular

[TABLE]

and hence also

[TABLE]

for $r$ large enough. Proposition 4.2 of [6] states

[TABLE]

Therefore, it follows that

[TABLE]

and

[TABLE]

Increasing $\epsilon$ appropriately, we obtain

[TABLE]

for both $b=(0,0)$ and $b=(1,0)$. Using the union bound yields

[TABLE]

for both $b=(0,0)$ and $b=(1,0)$. To conclude the proof for $b=(0,1)$ as well, we observe

[TABLE]

where the final inequality follows from Lemma 5.1 in [6] with

[TABLE]

The proof of (E.30) is analogous to (E.29), replacing the arrival processes by the service processes. Equations (E.21) and (E.22) are then a direct consequence of (E.27) and (E.28).

Using the previous result, we can show that $\mathbb{X}$ is almost Lipschitz continuous.

Proposition E.3**.**

Fix $\epsilon>0$, $M>0$ and $T>0$. For $r$ large enough,

[TABLE]

*where $N<\infty$ is constant (depending only on $\lambda$, $\mu$ and $\{p_{ij};\{i,j\}\in E\}$).

Proof.

This follows in a straightforward way from Lemma E.2 and the hydrondynamically-scaled system equations. That is, let $\mathcal{V}^{r}$ denote the intersection of the complements of the events given in equations (E.19)-(E.22), so $\mathbb{P}(\mathcal{V}^{r})\leq 1-N_{0}\epsilon$ 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 $\omega\in\mathcal{V}^{r}$, and for every $t_{1},t_{2}\in[0,T]$ and $m<\sqrt{r}T$,

[TABLE]

where $N_{1}$ and $N_{2}$ are only dependent on the system parameters (i.e. $\lambda$, $\mu$, $p$). Let $t_{1},t_{2}\in[0,T]$ with $t_{1}\leq t_{2}$. By definition of $\mathcal{V}^{r}$,

[TABLE]

and

[TABLE]

for $N$ as in Lemma E.2. Due to (E.10),

[TABLE]

In view of (E.12) and (E.10), we observe

[TABLE]

Due to (E.21),

[TABLE]

and similarly, due to (E.22),

[TABLE]

In view of (E.16),

[TABLE]

Finally, due to (E.17),

[TABLE]

We conclude that (E.31) is satisfied, as each process in $\mathbb{X}^{r,m}$ satisfies this property.

As is done in [6, 9], one can take $\epsilon$ appropriately small for every system. That is, for fixed $M>0$, $N>0$ and $T>0$, let

[TABLE]

where $\epsilon(r)\rightarrow 0$ as $r\rightarrow\infty$ is a sequence of positive real numbers. Moreover, in view of Lemma E.1, let for that same sequence $\{\epsilon(r)\}_{r\in\mathbb{R}}$,

[TABLE]

Let $\mathcal{K}^{r}$ denote the intersection of $\mathcal{K}_{0}^{r}$, $\mathcal{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 $\epsilon(r)$ if $\epsilon(r)\rightarrow 0$ sufficiently slow. We conclude that $\mathbb{P}(\mathcal{K}^{r})\rightarrow 1$ as $r\rightarrow\infty$.

Corollary E.4**.**

Fix $M>0$ and choose $N>0$ and $\epsilon(r)$ as above. Then,

[TABLE]

Following the framework of [6], we can use these results to state that the hydrodynamically scaled system convergences to a hydrodynamic limit. Fix $M>0$ and let $\tilde{E}$ be the set of right-continuous functions $x:[0,M]\rightarrow\mathbb{R}^{d}$ with left limits. Let

[TABLE]

Moreover, we set

[TABLE]

and

[TABLE]

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 $\mathcal{E}$ is a point $x\in\tilde{E}$ such that for all $\epsilon>0$ there exists a $r_{0}\in\mathbb{N}$ so that for every $r\geq r_{0}$ there is some $y\in E^{r}$ such that $\lvert x(\cdot)-y(\cdot)\rvert_{M}<\epsilon$. *

Since $|\mathbb{X}^{r,m}(0)|\leq 1$, the following result is a consequence of Proposition 4.1 in [6].

Corollary E.6**.**

Let $\tilde{E},E^{r},\mathcal{E}$ be as above. Fix $\epsilon>0$, $M>0$, $T\geq 0$, and choose $r$ large enough. Then, for $\omega\in\mathcal{K}^{r}$ and any $m<\sqrt{r}T$,

[TABLE]

*for some hydrodynamic limit $\mathbb{X}(\cdot)\in E^{\prime}$ of $\mathcal{E}$.

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 let $\tilde{\mathbb{X}}$ be a hydrodynamic limit of $\mathcal{E}$ over $[0,M]$. Then $\tilde{\mathbb{X}}$ satisfies the following equations:

[TABLE]

Remark E.8**.**

We cannot provide such general equations for $\tilde{Z}(\cdot)$ or $\tilde{L}(\cdot)$, since these limits depend on $x_{r,m}$. That is, the processes $\tilde{Z}^{r,m}(\cdot)$ and $\tilde{L}^{r,m}(\cdot)$ 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.

Proof of Proposition E.7.

Let $\tilde{\mathbb{X}}$ be a hydrodynamic limit of $\mathcal{E}$. For a given $\delta>0$, choose $(r,m)$ such that $\epsilon(r)\leq\delta$, and

[TABLE]

Due to (E.14), (E.15) and $\omega\in\mathcal{H}^{r}$, we derive

[TABLE]

From (E.21) and (E.22), we obtain

[TABLE]

and

[TABLE]

Equation (E.32) is a clear consequence of (E.10). Combining the above equations, we observe

[TABLE]

These bounds imply that any hydrodynamic limit satisfies equations (E.32)-(E.37). Finally, we still have to show (E.40). If for some $t\in[0,M]$,

[TABLE]

then by continuity of $\tilde{\mathbb{X}}$, there exists a $\eta>0$ such that this holds for all $s\in[t-\eta,t+\eta)$, and also

[TABLE]

Due to (E.18), this implies that $A_{ij,i}^{r,m}(s)$ is constant on $s\in[t-\eta,t+\eta])$. Therefore, its limit is also constant on $[t-\eta,t+\eta]$, and hence the derivative is zero. On the other hand, if

[TABLE]

then by continuity of $\tilde{\mathbb{X}}$, there exists a $\eta>0$ such that this holds for all $s\in[t-\eta,t+\eta])$, and

[TABLE]

Since $\tilde{A}^{\prime}_{ij,j}=0$, and due to (E.10) with limiting process (E.33),

[TABLE]

E.2 The SSC function

In this section, we introduce the state space collapse (SSC) function under which we show multiplicative state space collapse. The SSC function we use in our paper is $g:\mathbb{R}^{S}\rightarrow\mathbb{R}$, defined as

[TABLE]

where $q=(q_{1},\ldots,q_{S})$. We note that $g(\cdot)$ is a non-negative continuous function and satisfies

[TABLE]

for every $\alpha>0$.

Lemma E.9**.**

Suppose $g:\mathbb{R}^{S}\rightarrow\mathbb{R}$ is defined as in (E.41). Then,

[TABLE]

for every hydrodynamic model solution $\tilde{\mathbb{X}}$ satisfying $|\tilde{\mathbb{X}}(0)|\leq 1$, where

[TABLE]

*with $h>0$ some constant that depends only on $\lambda$ and $\{p_{ij},\{i,j\}\in E\}$. Moreover, if $g(\tilde{Q}(0))=0$ and $|\tilde{\mathbb{X}}(0)|\leq 1$, then $g(\tilde{Q}(t))=0$ for all $t\geq 0$.

Proof.

The proof relies heavily on the ideas used in the proof for the fluid limit. Let

[TABLE]

Since

[TABLE]

for any non-empty $\mathcal{I},\mathcal{J}\subset\{1,\ldots,S\}$ with $\mathcal{I}\cap\mathcal{J}=\emptyset$, we observe that $h<0$. For a hydrodynamic limiting process $\tilde{\mathbb{X}}$, let $H_{\tilde{\mathbb{X}}}(\cdot)$ be given by

[TABLE]

We note that this function is non-negative, and satisfies $H_{\tilde{\mathbb{X}}}(t)=0$ for all $t\geq g(\tilde{Q}(0))/h$.

To show that $g(\tilde{Q}(t))$ is bounded by this function, we note that it suffices to show that whenever $g(\tilde{Q}(t))>0$ with $t\geq 0$ being a regular point of $\tilde{\mathbb{X}}$,

[TABLE]

For this purpose, let

[TABLE]

and

[TABLE]

Due to Lemma 2.8.6 in [8], it holds for all $i,j\in S_{\max}(t)$ that $\tilde{Q}^{\prime}_{i}(t)/p_{i}=\tilde{Q}^{\prime}_{j}(t)/p_{j}$, and similarly, for all $i,j\in S_{\min}(t)$ it holds that $\tilde{Q}^{\prime}_{i}(t)/p_{i}=\tilde{Q}^{\prime}_{j}(t)/p_{j}$. Therefore, due to hydrodynamic limit equations (E.32)-(E.40) and the observation that there is at least one station $j\not\in S_{\max}(t)$ such that $\{i,j\}\in E$, it follows that

[TABLE]

for all $i\in S_{\max}(t)$. Similarly, for all $i\in S_{\min}$,

[TABLE]

We conclude that $g^{\prime}(\tilde{Q}(t))\leq-h$, and hence $g(\tilde{Q}(t))\leq H_{\tilde{\mathbb{X}}}(t)$ for all $t\geq 0$. In particular, if $g(\tilde{Q}(0))=0$ and $|\tilde{\mathbb{X}}(0)|\leq 1$, it follows from the definition of $H_{\tilde{\mathbb{X}}}(\cdot)$ that $g(\tilde{Q}(t))=0$ for all $t\geq 0$.

The first statement of the lemma follows since for every hydrodynamic model solution $\tilde{\mathbb{X}}$ satisfying $|\tilde{\mathbb{X}}(0)|\leq 1$, it holds that $g(\tilde{Q}(0))\leq 2/\min_{1\leq j\leq S}p_{j}$. Hence, $H_{\tilde{\mathbb{X}}}(\cdot)\leq H(\cdot)$ for every hydrodynamic model solution $\tilde{\mathbb{X}}$ satisfying $|\tilde{\mathbb{X}}(0)|\leq 1$.

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 $\epsilon>0$, $M>0$ and $T>0$. Then, for every $\omega\in\mathcal{K}^{r}$,

[TABLE]

for all $t\in[0,M]$ and $m<\sqrt{r}T$, where $H(\cdot)$ is as in Lemma E.9. Moreover, if $g(\hat{Q}(0))\rightarrow 0$ in probability as $r\rightarrow\infty$, then for all $\omega\in\mathcal{L}^{r}$ with

[TABLE]

it holds that

[TABLE]

and

[TABLE]

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.41). To do so, we first need to relate the hydrodynamic and diffusion scaling. That is, we observe that

[TABLE]

where

[TABLE]

with

[TABLE]

Corollary E.10 can be translated to the diffusion scaled process. Consider the SSC function $\hat{g}:\mathbb{R}^{S}\rightarrow\mathbb{R}$ defined as

[TABLE]

with $q=(q_{1},\ldots,q_{S})$.

Corollary E.11**.**

Fix $\epsilon>0$, $M>0$ and $T>0$. Then for $r$ large enough, and $\omega\in\mathcal{K}^{r}$,

[TABLE]

for all $t\in[0,T]$ with $m\in\mathbb{N}$ such that

[TABLE]

Also, for all $\omega\in\mathcal{L}^{r}$,

[TABLE]

Since $H(\cdot)$ is given as in (E.42), we observe that $H(t)=0$ for all $t\geq 2/(h\min_{1\leq j\leq S}p_{j})$. We would like to show that $(\sqrt{r}t-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.

Lemma E.12**.**

Suppose $M\geq 2(N+2)$ is fixed, and let

[TABLE]

Then for $r$ large enough,

[TABLE]

*for every $\omega\in\mathcal{K}^{r}$ and $t\in(My_{r,0}/\sqrt{r},T]$.

Proof.

For every $\omega\in\mathcal{K}^{r}$, by definition of the set,

[TABLE]

for $t_{1},t_{2}\in[0,M]$ and $m<\sqrt{r}T$. In particular, for $t_{2}=1/y_{r,m}$, $t_{1}=0$ and $\epsilon\leq 1$,

[TABLE]

Applying the reverse triangle inequality, we observe

[TABLE]

and similarly,

[TABLE]

Therefore, it always holds that

[TABLE]

For every $t\in(My_{r,0}/\sqrt{r},T]$, it follows by definition of $m_{r}(t)$ that

[TABLE]

In particular,

[TABLE]

where the last inequality follows since $M\geq 2(N+2)$.

Next, we show the main result of this section.

Theorem E.13**.**

Suppose $\hat{Q}^{r}(0)\rightarrow\hat{Q}(0)$ for some random vector $\hat{Q}(0)$. For every $T>0$, $\epsilon>0$ and $M<\infty$ with

[TABLE]

it holds that

[TABLE]

as $r\rightarrow\infty$. If, in addition, $\hat{g}(\hat{Q}^{r}(0))\rightarrow 0$ in probability as $r\rightarrow\infty$, then for every $T>0$,

[TABLE]

*as $r\rightarrow\infty$.

Proof.

Fix $\eta>0$ and note that by construction there exists a $r_{0}$ such that for all $r>r_{0}$,

[TABLE]

For every $\omega\in\mathcal{K}^{r}$, we have derived bounds that only require that $M$ is bounded. We note that for $M$ as in the statement of the theorem allows for every $t\in[0,T]$, $m/\sqrt{r}\leq t\leq(m+y_{r,m}M)/\sqrt{r}$ for some $m<\sqrt{r}M$. Moreover, it follows from Lemma E.12 and (E.42), that

[TABLE]

for all $t\in(My_{r,0}/\sqrt{r},T]$. In view of Corollary E.11, we obtain for every $\epsilon>0$,

[TABLE]

for all $t\in(My_{r,0}/\sqrt{r},T]$. Since for all $t\in[0,T]$,

[TABLE]

we obtain for every $\omega\in\mathcal{K}^{r}$,

[TABLE]

Note that for every $t\geq 0$,

[TABLE]

Moreover, since $\epsilon>0$ is arbitrary, we can conclude that (E.43) holds.

If $|\hat{Q}^{r}(0)|\overset{\mathbb{P}}{\rightarrow}0$, it follows from Corollary E.11, that for all $\omega\in\mathcal{L}^{r}$ and $t\in[0,My_{r,0}/\sqrt{r}]$

[TABLE]

Since $\epsilon>0$ is arbitrary, together with (E.43) and (E.45), we obtain (E.44).

Remark E.14**.**

Note that the bounds in Theorem E.13 are obtained for every fixed $T>0$. Yet, from the proof it is clear that one has the following slightly more general result. Suppose $\hat{Q}^{r}(0)\rightarrow\hat{Q}(0)$ for some random vector $\hat{Q}(0)$. For every $\epsilon>0$, $M<\infty$ as in Theorem E.13, and $t_{r}\in(My_{r,0}/\sqrt{r},\infty)$,

[TABLE]

as $r\rightarrow\infty$. If, in addition, $\hat{g}(\hat{Q}^{r}(0))\rightarrow 0$ in probability as $r\rightarrow\infty$, then for every $t_{r}\in(My_{r,0}/\sqrt{r},\infty)$,

[TABLE]

as $r\rightarrow\infty$. In other words, the interval over which the state space collapse is considered can also be chosen as a sequence of intervals indexed by $r$.

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.5. To do so, it suffices to show that the denominators in Theorem E.13 are bounded in a probabilistic sense. More specifically, $\lVert\hat{Q}^{r}(t)\rVert_{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.

Lemma E.15**.**

Suppose $\hat{Q}^{r}(0)\rightarrow\hat{Q}(0)$ for some random vector $\hat{Q}(0)$, and $M\in[0,\infty)$. Then,

[TABLE]

Proof.

Fix $\epsilon\in(0,1)$ small. First, note that

[TABLE]

By definition,

[TABLE]

Since $\hat{Q}^{r}(0)\rightarrow\hat{Q}(0)$ for some random vector $\hat{Q}(0)$,

[TABLE]

and due to the definition of $y_{r,0}$ and (E.45) for $t=0$, this implies

[TABLE]

To bound the first term in (E.48) as well, we observe that the queue length at some time is trivially bounded by

[TABLE]

We observe that $F_{j}^{r}\leq(1+\epsilon)\lambda r/\mu$ for $r$ large enough. Moreover, if $\{\Lambda(t),t\geq 0\}$ denotes a Poisson process with rate $\lambda$, then due to the properties of the Poisson process it holds that $\sum_{\{i,j\}\in E}\Lambda_{ij}(\cdot)\overset{d}{=}\Lambda(\cdot)$. In terms of the diffusion scaling, the above bound yields for all $t\geq 0$,

[TABLE]

Therefore, using this bound for $t=My_{r,0}/\sqrt{r}\leq\epsilon MK/\sqrt{r}$ and noting that Poisson processes are (non-decreasing) counting processes,

[TABLE]

Due to the LLN, we observe that

[TABLE]

for $\epsilon>0$ small enough (e.g. for $\epsilon<1-M/(M+\sqrt{\lambda/\mu}\min_{1\leq j\leq S}p_{j})$). Similarly, due to the LLN,

[TABLE]

for $\epsilon>0$ small enough (e.g. for $\epsilon<\min_{1\leq j\leq S}p_{j}/(M\sqrt{\lambda/\mu}+\min_{1\leq j\leq S}p_{j})$). We conclude that the first term in (E.48) converges to zero, i.e.

[TABLE]

and hence the result follows.

Next, we show that the process $\hat{Q}^{r}(\cdot)$ satisfies the compact containment property.

Proposition E.16**.**

Suppose $|\hat{Q}^{r}(0)|\rightarrow\hat{Q}(0)$ for some random vector $\hat{Q}(0)$. Then, for every $T>0$ and $\epsilon>0$,

[TABLE]

Proof.

Fix $\epsilon\in(0,1/3)$ small, and let $K>\max\{2|\beta|+2,2|\gamma|+2\}$. Introduce the sequence of stopping times

[TABLE]

and similarly,

[TABLE]

Clearly,

[TABLE]

In order to improve the readability of the proof, we first present a proof in the case when $\hat{g}(\hat{Q}^{r}(0))\rightarrow 0$ in probability as $r\rightarrow\infty$. We then comment on the changes needed to adapt the proof to the case when this condition does not necessarily hold.

Case I: $\hat{g}(\hat{Q}^{r}(0))\rightarrow 0$ in probability as $r\rightarrow\infty$.

Since $|\hat{Q}^{r}(0)|\rightarrow\hat{Q}(0)$ for some random vector $\hat{Q}(0)$, it holds that

[TABLE]

and hence we can assume that both $\hat{\tau}_{K}^{r}>\hat{T}_{K}^{r}>0$ and $\breve{\tau}_{K}^{r}>\breve{T}_{K}^{r}>0$ (for $K$ large enough). Moreover, for every $t\leq\min\{\breve{\tau}_{K}^{r},\hat{\tau}_{K}^{r}\}$ ,

[TABLE]

Next, we provide bounds for the two ways of crossing the boundary $K$ separately. First, we consider the first term in (E.50). We observe

[TABLE]

Due to Theorem E.13 and (E.47),

[TABLE]

To bound the first term, define the process $\{\hat{Q}^{r}_{\Sigma}(t),t\geq 0\}$ with

[TABLE]

We observe that

[TABLE]

Due to the system identities, we observe that for every $t\in[\hat{T}_{K}^{r},\hat{\tau}_{K}^{r}]$,

[TABLE]

We note that due to the properties of the Poisson process,

[TABLE]

with $\{\Lambda(t),t\geq 0\}$ an (independent) Poisson process with rate $\lambda$. Moreover, since for all $t\in[\hat{T}_{K}^{r},\hat{\tau}_{K}^{r}]$ it holds that $\hat{Q}^{r}_{j}(t)\geq\gamma$ for every $i\in\{1,\ldots,S\}$,

[TABLE]

Using the FCLT, we observe that

[TABLE]

as $r\rightarrow\infty$, where $\{\textrm{BM}(t),t\geq 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.51) and (E.52), for all $t\in[\hat{T}_{K}^{r},\hat{\tau}_{K}^{r}]$,

[TABLE]

We conclude that

[TABLE]

which converges to zero as $K\rightarrow\infty$ since $\epsilon\in(0,1/3)$.

The analysis of the second term in (E.50) uses similar arguments as the first term. We observe

[TABLE]

Again, due to Theorem E.13 and (E.47),

[TABLE]

Due to the system identities, we observe that for every $t\in[\breve{T}_{K}^{r},\breve{\tau}_{K}^{r}]$,

[TABLE]

Due to the definitions of the stopping times, we observe that for all $t\in[\breve{T}_{K}^{r},\breve{\tau}_{K}^{r}]$, it holds that $\hat{Q}_{j}^{r}(t)\leq\beta$ for every $j\in\{1,\ldots,S\}$, and hence $L^{r}(t)=r$. Therefore, due to the properties of the Poisson process,

[TABLE]

with $\{\Lambda(t),t\geq 0\}$ a Poisson process with rate $\lambda$. Moreover,

[TABLE]

Using the FCLT, we observe again that

[TABLE]

as $r\rightarrow\infty$, where we recall that $\{\textrm{BM}(t),t\geq 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.51) and (E.52), it holds for all $t\in[\breve{T}_{K}^{r},\breve{\tau}_{K}^{r}]$,

[TABLE]

We conclude that

[TABLE]

which also converges to zero as $K\rightarrow\infty$ since $\epsilon\in(0,1/3)$. Since this holds for both of the two summed probabilities in (E.50), we conclude that (E.49) holds.

Case II: general case, i.e. when we do not assume that $\hat{g}(\hat{Q}^{r}(0))\rightarrow 0$ in probability as $r\rightarrow\infty$.

Let $M\in[1,\infty)$ be fixed and satisfy the property as in E.13. Since $|\hat{Q}^{r}(0)|\rightarrow\hat{Q}(0)$ for some random vector $\hat{Q}(0)$ and due to Lemma E.15, it holds that

[TABLE]

Therefore, we can assume that both $\hat{\tau}_{K}^{r}>\hat{T}_{K}^{r}>My_{r,0}/\sqrt{r}$ and $\breve{\tau}_{K}^{r}>\breve{T}_{K}^{r}>My_{r,0}/\sqrt{r}$ (for $K$ large enough). The proof in this general case is then completely analogous to that in the previous one: $\lVert\hat{g}(\hat{Q}^{r}(t))\rVert_{\breve{\tau}_{K}^{r}}$ needs to be replaced with $\sup_{t\in(My_{r,0}/\sqrt{r},\breve{\tau}_{K}^{r}]}\lvert\hat{g}(\hat{Q}^{r}(t))\rvert$, and $\lVert\hat{g}(\hat{Q}^{r}(t))\rVert_{\hat{\tau}_{K}^{r}}$ – with $\sup_{t\in(My_{r,0}/\sqrt{r},\hat{\tau}_{K}^{r}]}\lvert\hat{g}(\hat{Q}^{r}(t))\rvert$.

Next, we prove our main result stated as in Theorem 4.5.

Proof of Theorem 4.5.

Equation (4.20) is a consequence of Theorem E.13 and Proposition E.16. To prove (4.19), note that for every $\epsilon>0$ and any sequence $\{K^{r},r\in\mathbb{N}\}$,

[TABLE]

Theorem E.13 and Proposition E.16 imply that for every $\epsilon>0$,

[TABLE]

Moreover, for any sequence $\{K^{r},r\in\mathbb{N}\}$ for which $K^{r}=o(\sqrt{r})$ with $K^{r}\rightarrow\infty$ as $r\rightarrow\infty$, it holds that

[TABLE]

by definition of $y_{r,0}$, (E.45) and since $\hat{Q}^{r}(0)\rightarrow\hat{Q}(0)$ for some random vector $\hat{Q}(0)$. We conclude that (4.19) holds as well.