Strong Embeddings for Transitory Queueing Models

TL;DR

This paper develops strong embedding theorems for nonstationary, non-Markovian transitory queueing models, enabling accurate diffusion approximations of their performance metrics despite inherent complexity.

Contribution

It introduces strong embedding results in the Komlos-Major-Tusnady framework for complex queueing models, facilitating precise sample path approximations.

Findings

Established strong embedding theorems for transitory queueing models

Provided error bounds for diffusion process approximations

Enhanced understanding of performance metrics in nonstationary queues

Abstract

In this paper we establish strong embedding theorems, in the sense of the Komlos-Major-Tusnady framework, for the performance metrics of a general class of transitory queueing models of nonstationary queueing systems. The nonstationary and non-Markovian nature of these models makes the computation of performance metrics hard. The strong embeddings yield error bounds on sample path approximations by diffusion processes, in the form of functional strong approximation theorems.

Equations

$$G_n(t):=\frac{1}{n}\sum _{i=1}^n{1}_{\{T_i\le t\}}.$$

$$A_n(t)=\sum _{i=1}^{n{G}_n(t)}{\zeta }_i.$$

$$r_n(G):=\underset{t\in [0,\infty )}{\sup }\left|G^{(n)}(t)-G(t)\right|=O\left(\frac{1}{\sqrt{n}}\right),$$

$$W_n(t)=\sum _{i=1}^{A_n(t)}{V}_i,$$

$$M_n(t):=\sup \left\{0\le m\le n:\sum _{i=1}^m{V}_i\le c_{n}t\right\}.$$

$$H_n(t)=\{\begin{array}{ll}npG(t)+\sqrt{n}\left(p{B}_{G(t)}^{\text{br},n}+\sqrt{p(1-p)}{\hat{B}}_{G(t)}\right),& \text{under Assum. refassum:arriv-dropouts},\\ np(G(t)+r_n(G))+\sqrt{n}\left(p{B}_{G(t)}^{\text{br},n}+\sqrt{p(1-p)}{\hat{B}}_{G(t)}\right),& \text{under Assum. refassum:arriv-dropouts-2,}\end{array}$$

$$\underset{t\in [0,\infty )}{\sup }\left|A_n(t)-H_n(t)\right|\stackrel{\text{a.s.}}{=}O\left(n^{1/4}\sqrt{\log n}\right).$$

$$\frac{1}{\sqrt{np}}\left(H_n(t)-npt\right)$$

$$=\sqrt{np}\left(G^{(n)}(t)-t\right)+\sqrt{p}O(1)+\frac{1}{\sqrt{p}}\left(p{B}_{G(t)}^{\text{br},n}+\sqrt{p(1-p)}{\hat{B}}_{G(t)}\right)$$

$$\stackrel{D}{=}\sqrt{np}\left(G^{(n)}(t)-t\right)+\sqrt{p}O(1)+(\sqrt{p}+\sqrt{1-p}){\int }_0^t\sqrt{G^{\prime }(s)}d{\bar{B}}_s-\sqrt{p}G(t){\bar{B}}_1,$$

$$=\sqrt{p}\hat{G}(t)+\sqrt{p}O(1)+\left(\sqrt{p}+\sqrt{1-p}\right){\bar{B}}_t-\sqrt{p}tZ.$$

$$R_n(t)=\sqrt{n}\sigma B_{pG(t)}+\mu H_n(t)$$

$$\underset{t\in [0,\infty )}{\sup }\left|W_n(t)-R_n(t)\right|\stackrel{\text{a.s.}}{=}O\left(n^{1/4}\sqrt{\log n}\right).$$

$$\underset{t\in [0,\infty )}{\sup }\left|\varphi (W_n-c_n\cdot \text{id})(t)-\varphi (R_n-c_n\cdot \text{id})(t)\right|\stackrel{\text{a.s}}{=}O\left(n^{1/4}\sqrt{\log n}\right),$$

$$Q_n(t)=A_n(t)-M_n(D_n(t)),$$

$$I_n(t):=t-D_n(t).$$

$$X_n(t)=H_n(t)-\frac{c_{n}t}{\mu }+\sqrt{n}\frac{\sigma }{\mu }{B}_{E_n(t)},$$

$$E_n(t)=\{\begin{array}{ll}\frac{c_{n}t}{n\mu }+{\inf }_{s\le t}\left(pG(s)-\frac{c_n}{n}\frac{s}{\mu }\right),& \text{under Assum. refassum:arriv-dropouts},\\ \frac{c_{n}t}{n\mu }+p{r}_n(G)+{\inf }_{s\le t}\left(pG(s)-\frac{c_n}{n}\frac{s}{\mu }\right),& \text{under Assum. refassum:arriv-dropouts-2}.\end{array}$$

$$\underset{t\in [0,\infty )}{\sup }\left|Q_n(t)-\varphi (X_n)(t)\right|\stackrel{\text{a.s.}}{=}O\left(n^{1/4}\sqrt{\log n}\right),$$

$$\underset{t\in [0,\infty )}{\sup }\left|Q_n(t)-\varphi (X_n)(t)\right|\stackrel{\text{a.s.}}{=}O\left(n^{1/4}\sqrt{\log c_n}\right).$$

$$\mathbf{P}\left(\frac{S_n}{\sqrt{n}}\le y\right)⟶\Phi (y)\text{as}n\to \infty ,$$

$$S_n(t)=\frac{1}{\sqrt{n}}\left(S_{[nt]}+X_{[nt]+1}+\left(nt-[nt]\right)\right).$$

$$h(S_n)\stackrel{d}{⟶}h(B),$$

$$\frac{S_k-B_k}{\sqrt{n\log \log n}}\stackrel{\text{a.s.}}{⟶}0,\text{as}k\to \infty .$$

$$\underset{0\le t\le 1}{\sup }\frac{S_n(t)-\frac{1}{\sqrt{n}}{B}_{nt}}{\sqrt{\log \log n}}\stackrel{\text{a.s.}}{⟶}0.$$

$$F_n(t)=\frac{1}{n}\sum _{i=1}^n{1}_{\{U_i\le t\}},t\in [0,1].$$

$${\alpha }_n(t)=\sqrt{n}\left(F_n(t)-t\right),$$

$$B_t^{\text{br}}=B_t-t{B}_1,t\in [0,1],$$

$$\underset{0\le t\le 1}{\sup }\left|{\alpha }_n(t)-B_t^{\text{br}}\right|\stackrel{a.s.}{=}O\left({\left(\frac{\log n}{n}\right)}^{1/4}{\left(\log n\log \log n\right)}^{1/4}\right).$$

$$\underset{1\le k\le n}{\sup }\left|S_k-B_k\right|=O\left(\log n\right).$$

$$\mathbf{P}\left(\underset{1\le k\le n}{\sup }\left|\sum _{i=1}^n{X}_i-B_n\right|>C\log n+x\right)<K{e}^{-\lambda x},$$

$$\mathbf{P}\left(\underset{s\in [0,1]}{\sup }\sqrt{n}\left|{\alpha }_n(s)-B_n^{br}(s)\right|>C\log n+x\right)<K{e}^{-\lambda x},$$

Full text

Authors supported by the National Science Foundation through grant CMMI/1636069.

1. Introduction

In this paper we establish strong embedding theorems, in the sense of the Komlós-Major-Tusnády (KMT) framework, for the performance metrics of a general class of transitory queueing models [HoJaWa2016a, BeHoLe2019]. Transitory queueing models assume a large, but finite population of customers arrive at the system over some time horizon. Examples of such systems include hospital surgery departments and clinics, subscription-based services such as video and game streaming, app-based ride-sharing/transportation and delivery services. In each of these cases, the pool of potential customers is known to the service provider a priori – due to appointments that are handed out to patients ahead of time in the healthcare examples, and to subscriptions/sign-ons in the case of streaming and app-based services. However, the finite pool of customers implies that transitory models are nonstationary, both in the sense that they are purely transient in nature and also because the model parameters can vary temporally. This makes the computation of the performance metrics rather difficult. Consequently, we seek to approximate the performance metric stochastic processes by simpler ones that capture their most vital temporal features. The strong embedding theorems in this paper yield probabilisitic error bounds between the discrete-event performance metric processes and simpler diffusion process approximations in terms of the population size $n$. Our results, therefore, will provide practitioners and engineers with a sense of how large the population size must be, and how the service effort must be concomitantly scaled, so that diffusion approximations can be confidently used for performance analysis and optimization.

Strong approximations were first used for studying time homogeneous queueing models in [Ro1980]; see the survey paper [Gl1998] for a comprehensive introduction to the use of strong approximations to $G/G/1$ queueing models in heavy-traffic. In general, the transient analysis of queueing models is rather complicated, and therefore a number of approximations have been developed in appropriate scaling regimes, typically by certain types of reflected diffusion processes [ChYa2001]. As queueing models can be expressed (approximately) as functionals of random walks, strong approximations are particularly useful in this application context since the driving random walks can be directly replaced by approximating Brownian motion processes. Strong approximation analysis yields rates of convergence and, consequently, rigorous justification of the heavy-traffic approximation on a sample path basis. Our results provide similar insights for a class of nonstationary queueing models under a population acceleration scaling framework.

We assume that the offered load to the queueing system is ‘time-of-day’ dependent and displays long-range correlations. The modeling and analysis of transitory queues is, in general, quite complicated and we operate under the simplifying assumption that the time-of-day and correlative effects are present solely in the traffic characteristics and that the service requirements of the arriving customers are independent and identically distributed. We propose a two-variable traffic model labeled RS$(G,p)$ wherein we impute the $i$th arriving customer (out of $n$) with the random variable tuple $(T_{i},\zeta_{i})$, where $T_{i}$ takes values in $[0,\infty)$ and $\zeta_{i}$ is binary. $T_{i}$ models the (potential) arrival epoch of customer $i$ and $\sum_{i=1}^{n}\zeta_{i}$ is the number of customers who actually enter the queue; here ‘RS’ stands for ‘randomly scattered.’ We assume that the tuples are independent and identically distributed over the population, and that $T_{1}$ follows a distribution $G$ and $\mathbb{E}\zeta_{1}=p$. We also assume that the service requirements are generally distributed with finite moment generating function in the neighborhood of zero and independent of the tuple. Consequently we label this the  $\textnormal{RS}(G,p)/G/1$ queue. We make the following contributions in this paper:

1. (1)

We prove functional strong approximation theorems (FSATs) for the workload and queue length performance metric processes of the  $\textnormal{RS}(G,p)/G/1$ queue in Theorems 2.10 and 2.11 (respectively). These FSATs yield sample path error bounds between the performance metrics and nonstationary reflected Brownian Bridge processes. The nonstationary Brownian Bridge processes reflect the fact the offered load is time-of-day dependent and has long range correlations. 2. (2)

The proofs of the FSATs are consequences of non-asymptotic functional strong embedding theorems (FSETs) proved for the RS$(G,p)$ traffic process in Proposition 5.11, the workload process in Proposition 6.3 and the queue length process in Proposition 7.6 that yield exponential probability bounds as a function of the population size. 3. (3)

Our proof of the non-asymptotic probabilistic bounds require Dvoretzky-Kiefer-Wolfowitz (DKW) [DKW] style inequalities for Brownian motion randomly time-changed by a stochastic pure jump process, proved in Proposition 4.12. This generalized DKW inequality may be of independent interest and useful in proving bounds for other types of models.

Commentary on main results

Our analysis leans on strong approximations for empirical processes and random walks [CsRe2014], but also requires substantial innovation. The FSAT in Theorem 2.10 is a consequence of Proposition 6.3 where we prove a strong embedding result for the workload process of a  $\textnormal{RS}(G,p)/G/1$ queue, under the assumption that the service times possess finite moment generating functions in a neighborhood around zero. We show that, with high probability, for a given fixed population size $n$ the sample paths of the workload process can be approximated by those of a reflected Brownian bridge process with time dependent drift and diffusion coefficients. Indeed, we show that the error rate is $O(n^{1/4}\sqrt{\log n})$. Next, the FSAT to the queue length process of the  $\textnormal{RS}(G,p)/G/1$ queue in Theorem 2.11 follows from Proposition 7.6. Paralleling the result in Proposition 6.3, we show that the approximating process is a reflected Brownian bridge process with time inhomogeneous drift and diffusion coefficients. However, the drift and diffusion coefficients are scaled versions of those observed in Proposition 6.3. We note that the analysis of the queue length strong embedding theorem is significantly more involved. The proofs of these results requires a careful construction of a DKW-style inequality for a time-changed Brownian motion process, which we did not find in the literature (see Proposition 4.12). Again, we show that the error rate for the queue length process is $O(n^{1/4}\sqrt{\log n})$.

Relation with prior transitory analyses.

Note that the RS($G,p$) model affords flexibility for modeling service systems where the pool of potential customers is known a priori. This typically includes systems where customers ‘subscribe’ to the service ahead of time – for example, clinics and surgical departments in hospitals where patients are given appointment times, or video and game streaming services with subscribing customers, or ridesharing and food delivery services where the pool of customers are those who have downloaded the smartphone app. In each of these cases, the service provider has knowledge of who the potential customers are, but not all customers will use the service on a given ‘day’. The randomized arrivals in the RS($G,p$) model accounts for this effect, which is ignored in the  $\Delta_{(i)}/G/1$ model where $\sum_{i=1}^{n}\zeta_{i}=n$ (rendering this variable redundant). Note that the RS($G,p$) model can be extended to a periodic traffic setting, as done in [GlHo2017], and the performance metric approximations can still be used in that setting.

The bibliography on the  $\Delta_{(i)}/G/1$ model now includes pointwise limit results [Ne2013, Lo1994], functional strong laws and central limits [HoJaWa2016a, HoJaWa2016b, GlHo2017, BeHoLe2019] and large deviations principles [Ho2017, GlHo2017]. In the population acceleration scaling limit, the results in [HoJaWa2016a, HoJaWa2016b] show that the limiting diffusion for the workload and queue length processes are regulated through a directional derivative reflection map [MaRa2010]. This limit can be recovered by the FSATs in Theorems 2.10 and  2.11, though the result in [HoJaWa2016a, HoJaWa2016b] holds under the weaker condition that the service requirements have two finite moments. However, extracting performance measures (such as moments of the workload/queue length) from the directional derivative reflected process is incredibly hard. Indeed, in [BeHoLe2019], a different critical scaling is used to show that the queue length converges to a reflected Brownian motion with parabolic drift when the arrival epoch distribution $G$ is exponential. On the other hand, in [GlHo2017] a special “critical” load condition is used to prove that the workload process is approximated by a reflected Brownian motion process. These limit processes can be recovered automatically from the FSATs proved in this paper, albeit at the cost of stronger conditions on the service requirements. We note, however, that with effort it is possible to extend the FSATs to cases where only $m>2$ moments are available.

Relation to FSATs for nonstationary models.

There is a large and growing literature on nonstationary queueing models covering the whole range of problems that confront the modeling of nonstationary service systems. A crucial difference between this large body of work and the growing literature on transitory models is that the former implicitly assumes an infinite population of customers, while transitory models are exclusively finite population. We cannot possibly do justice to the large body of work on nonstationary models; see [Wh2018] for a recent review. Instead, we focus on strong approximation results that are most closely related. To the best of our knowledge, strong approximations have been proved almost exclusively for Markovian nonstationary models; note that the literature on strong approximations for stationary queueing networks is far more extensive. The most influential papers in this genre are [MaMa1995, MaMaRe1998], where the important uniform acceleration scaling regime was introduced. In the former, strong approximations for Markov processes were leveraged to prove an FSAT (and consequently functional strong laws and central limit theorems) for an isolated time-varying Markovian single-server queue. This analysis was significantly generalized in the latter paper to include multi-server queueing networks with abandonment. In [MaPa1995], strong approximations were leveraged to prove functional limits for state-dependent, nonstationary Markovian queues. More recently [KoPe2018] consider nonstationary Markovian arrival processes (MAPs) as models of the traffic, and develop a bespoke Poisson representation of the MAP process. They then exploit the strong approximations in [MaMaRe1998] to prove functional strong laws and central limit theorems. All of these results are premised on the availability of strong approximation results for Markov processes (see, for instance, [EtKu2009, Chapter 7]). However, the performance metric processes for the  $\textnormal{RS}(G,p)/G/1$ queue are not Markov (though, of course, one could do state-space enlargement) and we therefore choose to develop the strong approximation results ‘from scratch.’ What is also somewhat remarkable is the fact that we are able to leverage strong approximation results proved for stationary random walks and empirical processes to study nonstationary queueing models without making explicit Markovian assumptions. We believe the methods highlighted in this paper can be used for analyzing other nonstationary stochastic models (such as nonstationary many-server queues, networks of nonstationary queues and even nonstationary multi-class queues).

The rest of the paper is organized as follows. We start with preliminaries and main results in Section 2. In Section 3 we provide a brief primer on the strong approximation methodology, particularly the coupling arguments that underly the KMT construction. We do so to make the paper self-contained and since the KMT construction is recondite and not widely understood. Next, we present the DKW-style inequality for controlling the error between the Brownian motion and a counterpart process stochastically time-changed by a jump process in Section 4. Section 5 presents strong embeddings for the RS$(G,p)$ traffic process. The strong embeddings for the workload and queue length processes are proved in Section 6 and Section 7 (respectively). We end with commentary and conclusions in Section 8.

2. Preliminaries and Main Results

2.1. A mechanistic model of queueing

Consider a single server, infinite buffer queue that is non-preemptive, non-idling, and starts empty. Service follows a first-come-first-served (FCFS) schedule. Let $n$ be the nominal number of customers applying for service. Customers independently sample an arrival epoch $T_{i}$, $i=1,\ldots,n$, from a common distribution function. In addition, all customers independently sample identical Bernoulli random variables $\zeta_{i}$, $i=1,\ldots,n$. Customer $i$ chooses to turn up at time $T_{i}$ only if $\zeta_{i}=1$; we call this the “dropout” variable. The arrival process is the cumulative number of customers that have arrived by time $t$. Let Bern($p$) represent the Bernoulli probability distribution with parameter $p$.

Assumption 2.1**.**

For every $n\geq 1$, let $T_{1},\ldots,T_{n}$ be iid samples from a general distribution with distribution function G. Denote $G_{n}$ to be the empirical distribution function given by:

[TABLE]

Let $\zeta_{1},\ldots,\zeta_{n}$ be iid samples from Bern($p$). Then the arrival process $A_{n}$ is given by:

[TABLE]

Remark 2.2*.*

We call $A_{n}$ in (2) as the RS$(G,p)$ traffic model. The $\Delta_{(i)}/G/1$ model introduced in [HoJaWa2016a] is a special case of Assumption 2.1, corresponding to $p=1$.

Remark 2.3*.*

Note that it is possible to consider other ways of modeling a random number of arrivals. However, the dropout model considered here is a mechanistic way of describing the traffic. Note that the model assumes each user will sample a potential time to arrive and a binary indicator that the customer will actually enter the queue at that time. See Section 8 for further discussion.

Remark 2.4*.*

Observe that while the nominal number of arrivals is $n$, the actual number of arrivals realized is random. This traffic model provides a mechanistic description of nonstationary arrivals: since the distribution $G$ is non-uniform (in general), the expected number of arrivals per-unit time $\mathbf{E}[A_{n}(t)]/t$ can be seen to equal $npG(t)/t$, by an application of Wald’s identity. This can be seen as a surrogate of an arrival rate that is clearly time-varying; note that we have not assumed that the distribution is differentiable and consequently defining the rate as the derivative of $\mathbf{E}[A_{n}(t)]$ is inappropriate. A crucial point to note is that this time dependency arises from microscopic behavior as opposed to a posited time dependency in the rate function. This stands in contrast with the vast majority of nonstationary models proposed in the literature where the model description starts with posited time-varying rate functions.

Sometimes it is useful to consider arrivals from a general distribution which in turn approaches the limiting distribution $G$ as $n\to\infty$.

Assumption 2.5**.**

For every $n\geq 1$, let $T_{1},\ldots,T_{n}$ be iid samples from a general distribution with distribution function $G^{(n)}$ which satisfies the following condition:

[TABLE]

for some strictly increasing and Lipschitz continuous distribution function $G$. In addition, assume that each $G^{(n)}$ is Lipschitz continuous and the Lipschitz coefficient increases at most polynomially in $n$. The arrival process $A_{n}$ is now defined similar to (2).

Remark 2.6*.*

For simplicity and ease of presentation we will assume that arrivals are supported on $[0,\infty)$, that is, $G(0)=G^{(n)}(0)=0$.

Next, let $\{V_{i},i\geq 1\}$ be a sequence of independent and identically distributed non-negative random variables. $V_{i}$ represents the service requirement in time units of the $i^{\textnormal{th}}$ potential customer who turns up into the system. We also assume that the sequence is independent of the arrival times $T_{i}$, $i=1,\ldots,n$ and the corresponding indicators of turning up $\zeta_{i}$, $i=1,\ldots,n$.

Assumption 2.7**.**

For every $n\geq 1$, let $V_{1},\ldots,V_{n}$ be iid samples from a distribution which admits existence of a moment generating function in a neighborhood of zero. Let $\mu$ and $\sigma$ denote the mean and standard deviation respectively of this distribution. Let

[TABLE]

denote the cumulative offered load to the system till time $t$.

We assume that the server efficiency is $c_{n}$, i.e. it completes $c_{n}$ jobs in unit time. Let $M_{n}(t)$ be the “truncated” renewal process counting the number of jobs that the server can complete by time $t$ if working continuously with efficiency $c_{n}$ (notice that only $n$ jobs arrive to the system):

[TABLE]

2.2. Functional strong approximations

In this section, we list the main results proven in the sequel. Strong approximation results are usually stated in terms of versions of the random variables we wish to approximate. In our case, we require versions of the random arrival times $T_{i}$, the indicators of turning up $\zeta_{i}$ and the service times $V_{i}$. In order to avoid repetition we do not mention this crucial requirement in the following theorem statements. However the same version suffices for each theorem below. Let us also note that is often customary in the literature to assume that the underlying probability space is rich enough to support the random variables as well as the approximating stochastic processes. Our first result provides a strong embedding for the arrival process. Its proof follows from the forthcoming Proposition 5.11.

Theorem 2.8**.**

There exists a Brownian motion $\hat{B}$, a Brownian bridge $B^{\textnormal{br},n}$ such that if $H_{n}$ be defined as:

[TABLE]

then

[TABLE]

Remark 2.9*.*

It is useful to contrast Theorem 2.8 with the setting in [Wh2016]. In the latter, traffic is modeled through a sequence of time-changed stochastic counting processes $\{A^{n}(t):=(N\circ\Lambda^{n})(t)\}$, where $N$ is a stationary stochastic counting process that satisfies an FCLT and $\Lambda^{n}$ is a posited cumulative arrival rate function that is assumed to be such that $\hat{\Lambda}_{n}(t):=n^{-1/2}(\Lambda^{n}(nt)-nt)$ satisfies $\hat{\Lambda}_{n}(t)\to\hat{\Lambda}(t)$ uniformly on compact sets of $[0,\infty)$ as $n\to\infty$, for some deterministic limit function $\hat{\Lambda}$. [Wh2016, Theorem 3.1] shows that the scaled traffic process $\hat{A}_{n}(t):=n^{-1/2}(A^{n}(nt)-nt)$ converges to a limit $B+\hat{\Lambda}$, where $B$ is a Brownian motion. A vital advantage of such a traffic model is that the stochasticity and the non-stationarities/time-dependencies are completely separated from each other in the limit.

On the other hand we do not see such a clean separation in $H_{n}$ immediately. However, suppose that Assumption 2.5 holds with $G(t)=t$ on $[0,1]$, and $n^{1/2}(G^{(n)}(t)-t)\to\hat{G}(t)$ uniformly on compact sets of $[0,\infty)$ as $n\to\infty$. Then, using the fact that $B^{\textnormal{br},n}_{t}\stackrel{{\scriptstyle D}}{{=}}\tilde{B}^{n}_{t}-t\tilde{B}_{1}^{n}$ for a standard Brownian motion process $\tilde{B}^{n}$ and  [Ok2013, Theorem 8.5.2] it follows that

[TABLE]

where $\bar{B}$ is a standard Brownian motion process. Theorem 2.8 immediately shows that, for the arrival process $A_{n}(t)$ at fixed $t\in[0,1]$, $(np)^{-1/2}(A_{n}(t)-npt)\Rightarrow\sqrt{p}\hat{G}(t)+(\sqrt{p}+\sqrt{(1-p)})\bar{B}_{t}+\bar{Z}$ as $n\to\infty$, where $\bar{Z}$ is a Gaussian random variable with mean $O(\sqrt{p})$ and standard deviation $\sqrt{p}t$. This is reminiscent of the limit in [Wh2016, Theorem 3.1], and shows that our framework can recover a separation of the macroscopic time-dependencies and the mesoscopic stochasticity. The setting in [Wh2016] is important as it forms the basis for a whole series of works around nonstationary queueing models (see the survey [Wh2018]). We also note that a more rigorous weak limit analysis for a specific choice of $G^{(n)}$ is presented in [GlHo2017].

Our next major result proves strong embeddings for the workload process. In particular, for the cumulative load to the system we have the following result, which follows from the forthcoming Propositions 6.1 and 6.3.

Theorem 2.10**.**

Along with the Brownian motion $\hat{B}$ and Brownian bridge $B^{\textnormal{br},n}$ as considered in Theorem 2.8, there exists a Brownian motion $B$ such that if $R_{n}$ be defined as:

[TABLE]

then

[TABLE]

Let $\phi$ be the reflection map functional given by $\phi(f)(t):=f(t)-\inf_{u\leq t}f(u).$ Then the total remaining workload at time $t$ can be expressed as $\phi(W_{n}-c_{n}\cdot\textnormal{id})(t)$ and this satisfies:

[TABLE]

where $\textnormal{id}:x\mapsto x$ is the identity map.

Finally Theorems 2.8 and 2.10 are used to prove a strong embedding for the queue length process, $Q_{n}$ that includes both any customer in service and all waiting customers. Recall that the queue length $Q_{n}(t)$ at time $t$ is the difference between the number of arrivals and the number of job completions before time $t$. Denoting by $D_{n}(t)$ the amount of time the queue stays busy till time $t$, the queue length can be expressed as:

[TABLE]

Finally, the idle time process of the server is given by

[TABLE]

The following theorem is a consequence of Proposition 7.6.

Theorem 2.11**.**

Let $B$, $\hat{B}$ be the Brownian motions $B^{\textnormal{br},n}$ the Brownian bridge processes as considered in Theorems  2.8 and 2.10. Let

[TABLE]

where

[TABLE]

Then the queue length $Q_{n}(t)$ satisfies:

[TABLE]

if $c_{n}=O(n^{p})$ for some $p>0$ and $\liminf_{n}c_{n}>0$. Else we have

[TABLE]

Remark 2.12*.*

Observe that the queue length spends more time near zero as the server efficiency becomes super polynomial in $n$, resulting in a greater approximation error.

Remark 2.13*.*

Theorems 2.10 and 2.11 show that the scaled workload process $Z_{n}:=\phi(W_{n}-c_{n}\cdot\textnormal{id})/n$ and the scaled queue length process $Q_{n}/n$ are both closely approximated by nonstationary reflected Brownian motion (RBM) processes on a sample path basis. These theorems also imply the results in [HoJaWa2016a, HoJaWa2016b, BeHoLe2019] where functional strong laws and central limit theorems were proved for the scaled processes when $p=1$. More importantly, the proofs of Propositions 6.1, 6.3 and 7.6 show that the spatial scale of the processes

3. Strong Embeddings: A Primer

Let $X_{1},X_{2},\ldots$ be iid random variables from a distribution with mean [math] and variance $1$. Let $S_{n}=\sum_{i=1}^{n}X_{i}$ denote the $n^{\textnormal{th}}$ partial sum. Then the classical central limit theorem states that:

[TABLE]

where $\Phi$ is the central normal cdf. Equation (8) states that the distribution of $\frac{S_{n}}{\sqrt{n}}$ approaches that of a standard normal as $n\to\infty$. A stochastic process analog of (8) was proved in [donsker]. Let the stochastic process $\{S_{n}(t);t\in[0,1]\}$ be constructed as follows for each $n\in\mathbb{N}$:

[TABLE]

Then $\{S_{n}(t),t\in[0,1]\}$ converges in distribution to $\{B(t),t\in[0,1]\}$ as $n\to\infty$, where $B$ is a standard Brownian motion. More precisely,

[TABLE]

for every continuous functional $h:C(0,1)\to\mathbb{R}$. Heuristically, equations (9) and (10) imply that for $n$ large enough $S_{[nt]}+X_{[nt+1]}(nt-[nt])$ is close in distribution to $\sqrt{n}B_{t}$. Utilizing the scaling property of Brownian motion and observing that $X_{[nt+1]}$ is negligible compared to $S_{[nt]}$ (for large $n$), we can concur that $S_{k}$ is approximately close to $B_{k}$ for all $k\in\{1,\ldots,n\}$. A bound on the difference of the two was provided in [strassen], who showed the existence of a probability space containing versions of all associated random variables and processes such that:

[TABLE]

Equation (11) can be restated in the following form:

[TABLE]

A close associate of the partial sums $S_{n}$ are the empirical distribution functions corresponding to a sample of iid random variables. Consider for simplicity a random sample $U_{1},U_{2},\ldots$ of iid $U[0,1]$ random variables. The empirical cdf is then given by:

[TABLE]

Observe that the random quantities $\mathbf{1}_{\{U_{i}\leq t\}}$ are iid with mean $t$ and variance $t(1-t)$. After proper scaling, and considering our previous discussion we expect the empirical process $\alpha_{n}$ given by:

[TABLE]

to be close to a normal random variable with variance $t(1-t)$. We also expect a convergence result akin to (10) in the process level. Recall that the standard Brownian bridge $B^{\textnormal{br}}$ is a stochastic process that may be defined as:

[TABLE]

for a Brownian motion $B$. Since $B^{\textnormal{br}}$ is a Gaussian process and $\textnormal{Var}(B_{t}^{\textnormal{br}})=t(1-t)$, $B^{\textnormal{br}}$ is a possible candidate for the stochastic process approximating the empirical process. Indeed this was proved to be true in a result analogous to (12) in [brillinger], who showed the existence of a probability space containing versions of all associated random variables and processes such that:

[TABLE]

This result immediately implies the analogue to (10), i.e., $\{\alpha_{n}(t),t\in[0,1]\}$ converges in distribution to $\{B_{t}^{\textnormal{br}},t\in[0,1]\}$.

Equations (12) and (13) are insightful and provide a rate of convergence of the partial sums and the empirical processes. However, these are not the best rates of convergence one can achieve. It was shown by  Komlós-Major-Tusnády in [KMT], that when $X_{i}$ is allowed to have a finite moment generating function in a neighborhood of [math]:

[TABLE]

A similar rate is enjoyed by the empirical processes of uniforms. These two results are stated below in Theorems 3.1 and 3.2, along with the novel construction (also known as the Hungarian method) of $X_{i}$’s and $U_{i}$’s from the Brownian motion and Brownian bridge respectively. A new and different approach in proving such embedding results has been provided in [chatterjee] for the simple symmetric random walk. We will use the terminology strong embedding for coupling an arbitrary random variable $W$ with a Gaussian random variable $Z$ so that $W-Z$ has exponentially decaying tails at the appropriate scale. Theorems 3.1 and 3.2 thus provide strong embeddings to the partial sums $S_{n}$ and the empirical processes $\alpha_{n}$. As alluded to in the Introduction we will apply these results to obtain strong embeddings for the performance metrics of a $\textnormal{RS}(G,p)/G/1$ queue.

3.1. Strong embedding of the random walk

We present the KMT theorem for the strong embedding of the random walk. Proof ideas and construction can be found in the Appendix.

Theorem 3.1**.**

Let $F$ be a distribution function with mean [math] and variance $1$. In addition, suppose the moment generating function corresponding to $F$, $R(t)=\mathbb{E}(e^{tX})$, $X\sim F$, exists in a neighborhood of [math]. Then, given a Brownian motion $B$, and using it, one can construct a sequence of random variables $X_{1},X_{2},\ldots$ which are independent and identically distributed to $F$. Furthermore, the partial sums of $X_{i}$’s are strongly coupled to the Brownian motion $B$ in the following sense. For every $n\in\mathbb{N}$ and $x>0$:

[TABLE]

where $C$, $K$ and $\lambda$ are positive constants depending only on $F$.

3.2. Strong embedding of the empirical process

We present the strong embedding result for the empirical process. Proof ideas and construction can be found in the Appendix.

Theorem 3.2**.**

There exists a probability space with independent $U[0,1]$ random variables $U_{1},U_{2},\ldots$ and a sequence of Brownian bridges $B_{1}^{\textnormal{br}},B_{2}^{\textnormal{br}},\ldots$ such that for all $n\geq 1$ and $x\in\mathbb{R}$:

[TABLE]

for some constants $C$, $K$ and $\lambda$. Here the empirical process $\alpha_{n}$ is given by

[TABLE]

and

[TABLE]

Remark 3.3*.*

The constants $C$, $K$ and $\lambda$ in Theorem 3.2 can be chosen as $C=100$, $K=10$ and $\lambda=1/50$. See [CsRe2014, Theorem 4.4.1] for more details.

Remark 3.4*.*

The KMT construction relies on the generation of a sample of $n$ uniforms $U_{1},\ldots,U_{n}$ from a Brownian bridge $B^{\textnormal{br},n}$. It can be seen from the construction that having obtained $\{U_{1},\ldots,U_{n}\}$, one is unable to obtain another $U_{n+1}$ such that the new set $\{U_{1},\ldots,U_{n+1}\}$ satisfies (16) with the same Brownian bridge. Instead it would be necessary to redo the construction. This necessitates the need for a different Brownian bridge $B^{\textnormal{br},n}$ for every $n$.

4. Control of Time-changed Brownian motion

Our analyses in subsequent sections provide strong embedding results for several queue length characteristics to corresponding diffusion approximations. In order to achieve those results we need a strong control over the difference between a Brownian motion evaluated at several $n$-level stochastic quantities and their corresponding fluid limits as $n$ goes to infinity (for example, the empirical distribution of arrival epochs against the true arrival distribution). In this section we present general results on bounding the difference between Brownian motion evaluated at some stochastic jump process and its fluid limit. The forthcoming Proposition 4.12 is rather general and might be of independent interest. We start by stating an assumption on the fluid limit.

Assumption 4.1**.**

For each $n\geq 1$, let $\xi_{n}:[0,\infty)\mapsto\mathbb{R}$ be a bounded Lipschitz continuous function; that is, there exists $c_{\xi_{n}}>0$ such that

[TABLE]

for all $s,t\in[0,\infty)$.

We also impose regularity conditions on the stochastic jump process along which our Brownian motion will be evaluated. These are collected in the following assumption.

Assumption 4.2**.**

Let ${\{T_{n}\}}_{n\geq 1}$ be a sequence of non-decreasing positive numbers. Let $\Xi_{n}:=\{\Xi_{n}(s);s\in[0,L_{n}]\}$ be a stochastic pure jump process defined on a common probability space $(\Omega,\mathcal{F},\mathbf{P})$ for every $n\geq 1$, such that almost surely the number of its jumps in $[0,L_{n}]$ is bounded above by $kn^{m}$ for some fixed positive constants $k$ and $m$. In addition, assume that:

[TABLE]

Denote $D=\limsup_{n}\sup_{s}|\Xi_{n}(s)|$.

In order to obtain a non-asymptotic probabilistic bound on the difference $|B_{\Xi_{n}(s)}-B_{\xi_{n}(s)}|$, where $\xi_{n}$ and $\Xi_{n}$ are introduced in Assumptions 4.1 and 4.2 respectively, we impose further conditions on the distribution of $|\Xi_{n}-\xi_{n}|$. In particular, we require a Dvoretzky-Kiefer-Wolfowitz (DKW) style inequality [DKW] for the tail distribution of $|\Xi_{n}-\xi_{n}|$.

Assumption 4.3**.**

For every $n\geq 1$, let $\xi_{n}$ and $\Xi_{n}$ be as considered in Assumptions 4.1 and 4.2. Let there be constants $k_{0}$, $k_{1}$, $k_{2}$, $k_{3}$ and $0<\gamma<4$ such that the following inequality holds for every $\varepsilon>0$:

[TABLE]

In addition, denote $\alpha_{n}$ by:

[TABLE]

In the forthcoming Proposition 4.6 and Proposition 4.7 we will show that the Assumptions 4.1, 4.2 and 4.3 are satisfied for the arrival process given in (2) and the truncated renewal process given in (5). In order to prove these two lemma’s, we first recall a few facts on sub-exponential random variables.

Lemma 4.4**.**

Let $X_{1},\ldots,X_{n}$ be iid copies of a random variable with mean $\mu$ such that there exist parameters $(\nu,m)$ satisfying:

[TABLE]

Then the following holds true:

1. (i)

[TABLE] 2. (ii)

[TABLE]

Proof.

(i) The result follows from the usual considerations for sub-exponential random variables (see e.g. [Wainwright, Section 2.1.3]). The main ingredient is a Chernoff-type approach to obtain:

[TABLE]

Optimization of the right hand side followed by a premultiplication by $2$ to obtain the two-sided tail bound yields the desired result (19).

(ii) Observe that $M_{k}=\sum_{i=1}^{k}X_{i}-k\mu$ is a martingale. In addition, $x\mapsto e^{\lambda x}$ is a convex function. Consequently $e^{\lambda M_{k}}$ is submartingale. Thus applying Doob’s martingale inequality we obtain:

[TABLE]

thus reducing our considerations to (21). The same arguments carry forward and we obtain (20). ∎

The assumption (18) in Lemma 4.4 holds for every random variable $X$ with a finite moment generating function in a neighborhood of [math]. This is a consequence of the following lemma.

Lemma 4.5**.**

Let $X$ be a random variable with mean $\mu$, whose moment generating function exists in a neighborhood of [math]. Then we have:

[TABLE]

where $\nu=\sqrt{2\textnormal{Var}(X)}$ and $m$ is given by the condition

[TABLE]

Proof.

Observe that the moment generating function of $X$ satisfies:

[TABLE]

Noticing $\mathbf{E}X=\mu$ and by Cauchy-Schwarz inequality we have:

[TABLE]

Hence for all $\lambda$ satisfying $\mathbf{E}[e^{2\lambda|X-\mu|}]<4$ we have:

[TABLE]

However, $1+\lambda^{2}\textnormal{Var}(X)\leq e^{\lambda^{2}\textnormal{Var}(X)}$ and consequently:

[TABLE]

This yields (22). ∎

Proposition 4.6**.**

Let Assumption 2.1 hold with $G$ being the uniform distribution function on $[0,1]$. Let $\hat{A}_{n}$ denote the corresponding arrival process in (2). Then Assumptions 4.1, 4.2 and 4.3 hold with

[TABLE]

and $L_{n}=1$ for all $n$.

Proof.

Note that $\xi_{n}$ is bounded and Lipschitz with $c_{\xi_{n}}=p$. Next, notice that $\hat{A}_{n}$ has at most $n$ jumps in $[0,1]$. In addition:

[TABLE]

It remains to prove a DKW type inequality for the difference $|\frac{\hat{A}_{n}(t)}{n}-pt|$. To that effect observe that

[TABLE]

where $X_{i}=\frac{\chi_{i}-p}{\sqrt{p(1-p)}}$ are iid random variables with mean [math] and variance $1$. Consequently, we have

[TABLE]

From the standard DKW inequality for empirical distributions [DKW], the second term has the standard exponentially decreasing bound given by:

[TABLE]

For the first term, observe that the $X_{i}$’s have a finite moment generating function $E[e^{\lambda X_{i}}]$ for all $\lambda$. Hence appealing to Lemma 4.4 and Lemma 4.5 we obtain:

[TABLE]

Consequently there exist constants $k_{1}$ and $k_{2}$ such that

[TABLE]

∎

Proposition 4.7**.**

Let Assumption 2.7 hold and $M_{n}$ be given by (5). Then for any sequence of non-decreasing positive reals $L_{n}$, Assumptions 4.1, 4.2 and 4.3 hold with

[TABLE]

Proof.

Note that $\xi_{n}$ is bounded and Lipschitz with $c_{\xi_{n}}=\frac{c_{n}}{n\mu}$. Next notice that $M_{n}$ has at most $n$ jumps in $[0,L_{n}]$. In addition

[TABLE]

It remains to prove a DKW-style inequality for the difference $|\frac{M_{n}(t)}{n}-(\frac{c_{n}}{n}\frac{t}{\mu})\wedge 1|$. Observe that for any $L_{n}$ positive:

[TABLE]

where $N(t)=\inf\{m\geq 0:\sum_{i=1}^{m}V_{i}>t\}$. By a change of variable the first term on the right hand side has a simpler representation upon which we have:

[TABLE]

where $\tilde{N}_{s}=\inf\{m\geq 0:\sum_{i=1}^{m}\frac{V_{i}}{\mu}>t\}$. From [Ho84, Lemma] and observing that $\tilde{N}(\frac{S_{n}}{\mu}+)=n$, we notice that

[TABLE]

where $\tilde{U}(s)=\sum_{i=1}^{[s]}\frac{V_{i}}{\mu}$. In addition

[TABLE]

where $\tilde{S}_{k}=\sum_{i=1}^{k}\frac{V_{i}}{\mu}$. Combining equations (26), (27) and (28) we obtain for all $\varepsilon>0$:

[TABLE]

Observe that $V_{i}$’s have a finite moment generating function in a neighborhood of [math]. Hence appealing to Lemma 4.4 and Lemma 4.5 we have for all $\varepsilon>0$:

[TABLE]

for some constants $k_{2}$ and $k_{3}$. Similarly, Lemma 4.4 and Lemma 4.5 also imply for all $\varepsilon>0$:

[TABLE]

Consequently we have constants $k_{1}$, $k_{2}$ and $k_{3}$ such that for all $\varepsilon>0$:

[TABLE]

∎

In the following lemma we obtain an upper bound on the expected value of $\alpha_{n}$ as denoted in (17). This result would be utilized in the forthcoming Lemma 4.11.

Lemma 4.8**.**

Let Assumptions 4.1, 4.2 and 4.3 hold. Then there exists a constant $C^{\prime}$ such that:

[TABLE]

Proof.

Since $\xi_{n}$ is bounded, without loss of generality let us assume

[TABLE]

As a consequence notice:

[TABLE]

Thus we have:

[TABLE]

where the last step is obtained by a change of variable. Now breaking the integral into parts and using Assumption 4.3 we obtain from above:

[TABLE]

for some constant $C^{\prime}$, where the penultimate step is obtained by a change of variable ($n^{\gamma}s^{2}\mapsto z$) and then use of the fact that the integral $\int_{0}^{\infty}\exp(-4k_{2}z\wedge k_{3}\sqrt{z}n^{\gamma/2})z^{-3/4}dz$ is bounded yields our desired result (30). ∎

Assumption 4.9**.**

Let $\xi_{n}$ and $\Xi_{n}$ be as considered in Assumptions 4.1 and 4.2. Let a standard Brownian motion $B$ be defined on the same probability space $({\Omega},{\mathcal{F}},\mathbf{P})$. Then let $\hat{f_{s}}$ be the Gaussian process defined on $[0,L_{n}]$ by

[TABLE]

In addition, let $\mathbb{P}$ denote the conditional probability on $(\Omega,\mathcal{F},\mathbf{P})$ given $\Xi_{n}$. Thus the conditional expectation $\mathbb{E}[Z]=\mathbf{E}[Z|\Xi_{n}(s);s\in[0,L_{n}]]$. Denote $\gamma_{n}:=\mathbf{E}[\sup_{t\in[0,L_{n}]}\hat{f}_{t}]$.

The key ingredient to find probabilistic bounds for $\hat{f}$ as alluded to above is the Borell-TIS inequality included below for completeness. As a first step we find the conditional expectation of $\sup\hat{f}$ given $\Xi_{n}$. This is obtained in the following lemma.

Lemma 4.10**.**

Let Assumptions 4.1, 4.2 and 4.9 hold. Then there exist constants $M$ and $\tilde{C}$ such that

[TABLE]

Proof.

The canonical metric for $\hat{f}$ in $(\hat{\Omega},\hat{\mathcal{F}},\mathbb{P})$ is given by

[TABLE]

Let $\hat{D}$ denote the diameter of $[0,L_{n}]$ with respect to the canonical metric, i.e.,

[TABLE]

Let $\hat{N}(\varepsilon)$ be the metric entropy defined by the smallest number of balls of diameter $\varepsilon$ (with respect to the canonical metric $\hat{d}$) that cover $[0,L_{n}]$. Then from [adler, Theorem 1.3.3] there exists a universal constant $M$ such that

[TABLE]

It can be easily shown that the canonical metric $\hat{d}$ as defined in (33) satisfies:

[TABLE]

If the numbers $\xi_{n}(s),\Xi_{n}(s),\xi_{n}(t),\Xi_{n}(t)$ be arranged in ascending order form the vector $(d_{1},d_{2},d_{3},d_{4})$, then it can be shown that

[TABLE]

This implies, for $s$ and $t$ such that $\Xi_{n}(s)=\Xi_{n}(t)$,

[TABLE]

In addition, note that ${\hat{D}}/2=\sup_{s,t\in[0,L_{n}]}\hat{d}(s,t)/2\leq(\sup_{t\in[0,L_{n}]}|\Xi_{n}(t)-\xi_{n}(t)|)^{1/2}/\sqrt{2}=\alpha_{n}$.

In order to obtain an upper bound to $\hat{N}(\varepsilon)$, recall as mentioned earlier, $\hat{d}(s,t)\leq\sqrt{c_{\xi_{n}}}\sqrt{|s-t|}$ whenever $\Xi_{n}(s)=\Xi_{n}(t)$. Since $\Xi_{n}$ has at most $kn^{m}$ points of discontinuity, there are at most $(kn^{m}+1)$ intervals where $\Xi_{n}$ is constant. Let these intervals be $R_{0},\ldots,R_{kn^{m}}$. Then $\hat{N}(\varepsilon)$ can be bounded above as follows:

[TABLE]

Thus, using (36) we get from (34):

[TABLE]

Observe that $\log(x+y)\leq\log(x+1)+\log(y)$ for $x\geq 0$ and $y\geq 1$. Consequently we obtain

[TABLE]

where in the last step we have used the fact that $\sqrt{x+y}\leq\sqrt{x}+\sqrt{y}$. From (37) and (4) we obtain that there exists a constant $\tilde{C}$ such that:

[TABLE]

∎

Having obtained the conditional expectation of $\sup\hat{f}$, we next obtain the unconditional expectation of $\sup\hat{f}$. This is achieved in the following result and would be a crucial ingredient in the proof of the forthcoming Proposition 4.12.

Lemma 4.11**.**

Let Assumptions 4.1, 4.2, 4.3 and 4.9 hold. Then there exists constant $C$ such that

[TABLE]

Proof.

Using integration by parts, and denoting $c_{\xi_{n}}L_{n}$ by $L_{n}^{\prime}$, the first term in (32) yields

[TABLE]

where the last step is obtained by a change of variable ($\sqrt{\log(L_{n}^{\prime}/\varepsilon^{2}+1)}\mapsto t$). It is readily checked that

[TABLE]

Consequently we obtain from (4):

[TABLE]

Now, (40) can be represented using the standard normal distribution function as follows:

[TABLE]

Thus, from (32) and (41) we obtain:

[TABLE]

Having completed the first step in our attempt to bound $\gamma_{n}$, we now proceed to obtain an upper bound to the right hand side of (42). The expectation of the first term:

[TABLE]

where ${g}$ is the function $\tilde{g}(y)=y\sqrt{\log(L_{n}^{\prime}/y^{2}+1)}$ restricted to the domain $[0,s_{g}]$ and $s_{g}$ is the point of global maxima of $\tilde{g}$, thereby making ${g}$ invertible. In addition it is readily checked by comparing values of $\tilde{g}^{\prime}$ that $s_{g}<\sqrt{L_{n}^{\prime}/(e-1)}$. Following the above steps we thus obtain:

[TABLE]

The first term on the right hand side in (4) can be bounded above using (41) as follows:

[TABLE]

It is readily checked by using the standard upper bound for normal tail probability that there exists a constant $C^{\prime}$ such that the right hand side is bounded above by $C_{1}^{\prime}\frac{\sqrt{\log(L_{n}^{\prime}\vee n)}}{n^{\gamma/4}}$ and thus we have:

[TABLE]

In order to bound the second term in (4), let us perform a change of variable manipulation, namely replace $n^{\gamma}t^{4}$ by $e^{-z}$. We now obtain from Assumption 4.3:

[TABLE]

Consequently we have:

[TABLE]

It is readily checked that the integral on the right is finite and thus we have:

[TABLE]

for some generic constant $C_{2}^{\prime}$. Using (44) and (46) in (4) we now obtain:

[TABLE]

where $C_{3}^{\prime}=C_{1}^{\prime}+C_{2}^{\prime}$. For the second term in (42), we use the bound on the normal tail probability, namely, $1-\Phi(t)\leq\frac{e^{-t^{2}/2}}{t\sqrt{2\pi}}$. Thus, we have

[TABLE]

The right hand side can be bounded above by using the fact that $\frac{1}{\sqrt{\log(L_{n}^{\prime}+\alpha_{n}^{2})}}\leq\frac{1}{\sqrt{\log T_{1}^{\prime}}}$, and the bound for $\mathbf{E}[\alpha_{n}]$ achieved in Lemma 4.8. Consequently, combining (30), (42) and (47) we have thus obtained

[TABLE]

for some constant $C$. ∎

We finally arrive at the main result of this section, namely, in Proposition 4.12 we state a general non-asymptotic probabilistic bound on the difference between a time-changed Brownian motion evaluated on a stochastic jump process and its fluid limit.

Proposition 4.12**.**

Let Assumptions 4.1, 4.2, 4.3 and 4.9 hold. Then there exist a constant $C$ such that for all $n\geq 1$ and $x>0$:

[TABLE]

where $\upsilon_{n}^{2}=\sup_{s\in[0,L_{n}]}\mathbf{E}[|\Xi_{n}(s)-\xi_{n}(s)|]$.

Proof.

Observe that we have:

[TABLE]

This implies

[TABLE]

Recall $\gamma_{n}:=\mathbf{E}\left[\sup_{t\in[0,1]}\hat{f_{t}}\right]$. Then by the Borell-TIS inequality (see e.g. [adler, Theorem 2.1.1]) we have:

[TABLE]

Now invoking Lemma 4.11 we obtain our desired result (48). ∎

Finally, we will require the following proposition for proving strong embeddings under Assumption 2.5. Consider the following regularity condition.

Assumption 4.13**.**

For every $n\geq 1$, let $\xi^{(n)},\xi:[0,T]\mapsto\mathbb{R}$ satisfy:

[TABLE]

In addition, let $\xi^{(n)}$ and $\xi$ both be Lipschitz continuous with the Lipschitz coefficient of $\xi^{(n)}$ growing at most polynomially in $n$.

Proposition 4.14**.**

Let Assumption 4.13 hold and $B^{n},n\geq 1$ be any sequence of Brownian motions defined on a probability space $({\Omega},{\mathcal{F}},\mathbf{P})$. Then there exists constants $C$, $K$ and $\lambda$ such that:

[TABLE]

Proof.

The key ingredient of the proof is again the Borel-TIS inequality. First, let us reuse the same notations as before; viz., let $\hat{f_{s}}$ be the Gaussian process defined on $[0,T]$ by

[TABLE]

As before, the canonical metric for $\hat{f}$ in $({\Omega},{\mathcal{F}},\mathbf{P})$ is given by (33) and let $\hat{D}$ denote the diameter of $[0,T]$ with respect to the canonical metric, i.e.,

[TABLE]

Let $\hat{N}(\varepsilon)$ be the metric entropy defined by the smallest number of balls of diameter $\varepsilon$ (with respect to the canonical metric $\hat{d}$) that cover $[0,T]$. Then from [adler, Theorem 1.3.3] there exists a universal constant $M$ such that

[TABLE]

It can be easily shown that the canonical metric $\hat{d}$ as defined in (33) satisfies:

[TABLE]

If the numbers $\xi(s),\xi^{(n)}(s),\xi(t),\xi^{(n)}(t)$ be arranged in ascending order form the vector $(d_{1},d_{2},d_{3},d_{4})$, then it can be shown that

[TABLE]

This implies that we have

[TABLE]

Consequently from Assumption 4.13 we have

[TABLE]

In addition, (54) also implies:

[TABLE]

Using the Lipschitz continuity of $\xi^{(n)}$ and $\xi$, we obtain that

[TABLE]

where $l_{n}$ grows polynomially in $n$. This implies $\hat{d}(s,t)\leq l_{n}\sqrt{|s-t|}$ for all $s,t\in[0,T]$. Thus $\{s\in[0,T]:|s-x_{0}|\leq\frac{\varepsilon^{2}}{l_{n}^{2}}\}$ is contained in the $\varepsilon-$ball around $x_{0}$. The length of this ball is thus at least $\frac{2\varepsilon^{2}}{l_{n}^{2}}$. Hence the number of $\varepsilon-$balls which cover $[0,T]$ is atmost $\frac{Tl_{n}^{2}}{2\varepsilon^{2}}$. This leads to an upper bound for $\hat{N}(\varepsilon)$, namely

[TABLE]

Thus, using (57) we get from (53):

[TABLE]

It is readily checked using integration by parts that

[TABLE]

A change of variable $\sqrt{\log\frac{Tl_{n}^{2}}{\varepsilon^{2}}}\mapsto t$ in the integral on the right hand side yields

[TABLE]

The standard upper bound to the normal tail probability now gives:

[TABLE]

Observe that since $l_{n}\geq 1$ and $\hat{D}\leq 2$ we have:

[TABLE]

In addition, since $x\log x\geq-1$ for all $x>0$:

[TABLE]

Using the inequalities (60) and (61) in the right hand side of (59) we get:

[TABLE]

Finally using (55) and (56) in (62) we obtain:

[TABLE]

Observe that

[TABLE]

Then by the Borell-TIS inequality (see e.g. [adler, Theorem 2.1.1]) there exists constants $C$, $K$ and $\lambda$ such that:

[TABLE]

∎

5. A Strong Embedding for the Arrival Process

In this section we derive a strong embedding for the arrival process. The following proposition is an extension of Theorem 3.1 when the length of the random walk is provided by a time-varying not necessarily determinstic function.

Proposition 5.1**.**

Let $X_{1},\ldots,X_{n}$ be iid samples from a distribution which admits existence of a moment generating function in a neighborhood of zero. Let $\mu$ and $\sigma$ denote the mean and standard deviation respectively of this distrbution. Let $J_{n}:[0,\infty)\mapsto\{1,\ldots,n\}$ be any process. Then there exists a standard Brownian motion $B$ and a version of $X_{1},\ldots,X_{n}$, along with constants $C_{1}$, $K_{1}$ and $\lambda_{1}$ such that for all $x>0$ we have:

[TABLE]

Proof.

From Theorem 3.1 we have that there exists a standard Brownian motion $\tilde{B}$, a version of $X_{1},\ldots,X_{n}$, along with constants $C$, $K$ and $\lambda$ (depending on the distribution of $V$) such that for all $x>0$ we have:

[TABLE]

where $\sum_{i=1}^{k}$ is defined to be the null sum for $k=0$. Since $A_{n}$ takes values in $\{1,\ldots,n\}$, we may replace the supremum in the left hand side of (63) by a supremum over $k$ taking values in $\{A_{n}(t),t\in[0,\infty)\}$. Consider another version $B$ of the standard Brownian motion $\tilde{B}$ such that

[TABLE]

Then from (63) there exist constants $C_{1}$, $K_{1}$ and $\lambda_{1}$ such that for all $x>0$ the desired strong-embedding holds:

[TABLE]

∎

The following proposition is a consequence of Theorem 3.2 and holds for any general distribution as opposed to the uniform distributional assumption made in Theorem 3.2.

Proposition 5.2**.**

Let Assumption 2.1 hold with $p=1$, that is, let us consider the $\Delta_{(i)}/G/1$ model as explained in Remark 2.2. Then for every $n\geq 1$ there exists a Brownian bridge $\{{B^{\textnormal{br},n}};t\in[0,1]\}$ and a version of $T_{1},\ldots,T_{n}$ along with constants $C_{2}$, $K_{2}$ and $\lambda_{2}$ such that for all $x>0$ we have:

[TABLE]

Proof.

We will first consider the random variables $\{G(T_{i}):i=1,\ldots,n\}$. Observe that the $G(T_{i})$’s are independent and identically distribued as $U[0,1]$ random variables. Consider the corresponding empirical distribution function $F_{n}$ given by:

[TABLE]

Then by a little modification of (15) there exist a Brownian bridge $B^{\textnormal{br},n}$ (observe by Remark 3.4, the Brownian bridge under consideration depends on $n$), constants $C_{2}$, $K_{2}$ and $\lambda_{2}$ such that:

[TABLE]

Let the inverse distribution function $G^{-1}$ be defined as:

[TABLE]

Observe that owing to our definition of $G^{-1}$ we have:

[TABLE]

Applying relation (68) in (65), and using (1) we obtain:

[TABLE]

Inserting (69) in (66) yields:

[TABLE]

In addition, observe that for any $s_{1}<s_{2}$ such that $G(s_{1})=G(s_{2})$ we have for all $i=1,\ldots,n$:

[TABLE]

This implies even though $G^{-1}(G(s))\geq s$, we still have:

[TABLE]

Consequently we obtain:

[TABLE]

We utilize this property in (70). Notice that

[TABLE]

Thus we have the following inequality between the suprema of the same function over these two sets:

[TABLE]

Inserting (71) in (72) we thus get

[TABLE]

Looking at the complement probability in the left hand side of (70) we obtain as a result of (73):

[TABLE]

This yields our desired result (64). ∎

Remark 5.3*.*

Observe that the constants $C_{2}$, $K_{2}$ and $\lambda_{2}$ in (66) do not depend on $G$ and that the same constants satisfy (64). Thus owing to Remark 3.3 we have that $C_{2}=100$, $K_{2}=10$ and $\lambda_{2}=1/50$ satisfy (64).

We now adapt the statement of Proposition 5.2 under Assumption 2.5.

Corollary 5.4**.**

Let Assumption 2.5 hold. Then for every $n\geq 1$ there exists a Brownian bridge $\{{\tilde{B}^{\textnormal{br},n}};t\in[0,1]\}$ and a version of $T_{1},\ldots,T_{n}$ such that for all $x>0$, the same constants $C_{2}$, $K_{2}$ and $\lambda_{2}$ as in Proposition 5.2 satisfy:

[TABLE]

Proof.

Observe from Remark 5.3 for every $k\geq 1$, there exists a Brownian bridge $B^{\textnormal{br},k,n}$ such that for all $x>0$, the same constants $C_{2}$, $K_{2}$ and $\lambda_{2}$ as in Proposition 5.2 satisfy:

[TABLE]

In particular, for $k=n$ and writing $B^{\textnormal{br},n,n}$ as $\tilde{B}_{G^{(n)}(t)}^{\textnormal{br},n}$ we have:

[TABLE]

∎

In the sequel, we will require control over Brownian motion evaluated at the fluid-scaled arrival process $A_{n}/n$ and the corresponding fluid limit. This is achieved for $U[0,1]$ distributed time epochs in the following proposition.

Proposition 5.5**.**

Let $T_{1},\ldots,T_{n}$ be iid samples from the $U[0,1]$ distribution. Let $\hat{A}_{n}$ be the arrival process with dropouts given by

[TABLE]

where $F_{n}$ is the empirical distribution function corresponding to the sample $T_{1},\ldots T_{n}$, and $\zeta_{i}$ are iid $\textnormal{Ber}(p)$. Let $B$ be a Brownian motion independent of $T_{1},\ldots,T_{n}$ and $\zeta_{1},\ldots,\zeta_{n}$. for each $n\geq 1$. Then, there exist constants $C_{3}$, $K_{3}$ and $\lambda_{3}$ such that

[TABLE]

Proof.

The proof follows from the DKW-type inequality established for the Brownian motion in Proposition 4.12, the conditions for which are satisfied in Proposition 4.6. Consequently we obtain there exists constants $C_{3}$, $K_{3}$ and $\lambda_{3}$ such that:

[TABLE]

where

[TABLE]

In order to bound $\sigma_{n}^{2}$, we first apply Cauchy-Schwarz inequality to get:

[TABLE]

Observe $\hat{A}_{n}(s)$ has the ${\rm Bin}(n,ps)$ distribution, which implies:

[TABLE]

Combining (76) and (77) we obtain from (75):

[TABLE]

Using inequality (78) in (74) we now have our desired strong embedding result:

[TABLE]

which holds for all $n\geq 1$ and $x>0$. ∎

We now extend our result in Proposition 5.5 to generally distributed time epochs. This is the subject of the following corollary.

Corollary 5.6**.**

Let Assumption 2.1 holds. Let $B$ be a Brownian motion independent of $T_{1},\ldots,T_{n}$ and $\zeta_{1},\ldots,\zeta_{n}$ for each $n\geq 1$. Then for every $n\geq 1$ and $x>0$, the same constants $C_{3}$, $K_{3}$ and $\lambda_{3}$ as in Proposition 5.5 satisfy:

[TABLE]

Proof.

The proof of this result is similar to the reasonings we adopted in Proposition 5.2. Consequently, let us again consider the random variables $\{G(T_{i}):i=1,\ldots,n\}$. Observe that the $G(T_{i})$’s are independent and identically distributed as $U[0,1]$ random variables. Consider the corresponding distribution function $F_{n}$ given by:

[TABLE]

Let $\hat{A}_{n}(t)=\sum_{i=1}^{nF_{n}(t)}\zeta_{i}$. In addition, recall the definition of $G^{-1}(t)$ in (67). Using (69) we obtain:

[TABLE]

In a spirit similar to what is used to obtain (73) we have from the analogue to (71):

[TABLE]

Our desired result now follows from Proposition 5.5. Observe Proposition 5.5 guarantees existence of constants $C_{3}$, $K_{3}$ and $\lambda_{3}$ such that, for all $n\geq 1$ and $x>0$ we have:

[TABLE]

We now complete our proof by looking at the complement probability in (81) and using (79) we have:

[TABLE]

Now using (80) in (5) we obtain:

[TABLE]

which yields our desired result. ∎

Remark 5.7*.*

Observe that the constants $C_{3}$, $K_{3}$ and $\lambda_{3}$ in Corollary 5.6 do not depend on the particular distribution $G$.

Corollaries 5.4 and 5.6 provide approximations in terms of $G^{(n)}$. In order to further simplify our approximation processes we require the following lemma.

Lemma 5.8**.**

Let Assumption 2.5 hold. Then for any $q>0$ and any sequence of Brownian motions $B^{n}$, there exist constants $C_{4}$, $K_{4}$ and $\lambda_{4}$ such that

[TABLE]

Proof.

We apply Proposition 4.14. Take $T=1$, $\xi^{(n)}(s)=qG^{(n)}(G^{-1}(s))$, $\xi(s)=qs$ and observe that continuity of $G$ implies

[TABLE]

Notice $G(G^{-1}(G(s)))=G(s)$ for all $s$. Thus we have:

[TABLE]

From (3) we now obtain that

[TABLE]

The Lipschitz continuity of $\xi$ is obvious, while Lipschitz continuity of $\xi^{(n)}$ follows from that of $G^{(n)}$ and $G^{-1}$. In addition, from a similar property for $G^{(n)}$, the Lipschitz coefficient of $\xi^{(n)}$ grows at most polynomially in $n$. Thus we have constants $C_{4}$, $K_{4}$ and $\lambda_{4}$ such that:

[TABLE]

Since $G$ is strictly increasing in $[0,\infty)$, we have $G^{-1}(G(s))=s$ for $s\in[0,\infty)$. Consequently we have:

[TABLE]

This provides our desired result from (84). ∎

We can now extend Corollary 5.6 courtesy Lemma 5.8.

Corollary 5.9**.**

Suppose Assumption 2.5 holds. Let $B$ be a Brownian motion independent of $T_{1},\ldots,T_{n}$ and $\zeta_{1},\ldots,\zeta_{n}$ for each $n\geq 1$. Then there exist constants $C_{5}$, $K_{5}$ and $\lambda_{5}$ such that for all $n\geq 1$ and $x>0$ we have:

[TABLE]

Proof.

From Corollary 5.6 and observing Remark 5.7 we have

[TABLE]

From Lemma 5.8 with $q=p$ we have constants $C_{4}^{\prime}$, $K_{4}^{\prime}$ and $\lambda_{4}^{\prime}$ such that

[TABLE]

Combining the above we have our desired result. ∎

Corollary 5.4 can be extended using the following result which again is a consequence of Lemma 5.8.

Corollary 5.10**.**

Let Assumption 2.5 hold. Then we have:

[TABLE]

Proof.

Using the fact that the Brownian bridge $\tilde{B}^{\textnormal{br},n}$ can be represented as:

[TABLE]

for a Brownian motion $B^{n}$ we have that

[TABLE]

From Lemma 5.8 with $q=1$, there exist constants $C_{4}^{\prime\prime}$, $K_{4}^{\prime\prime}$ and $\lambda_{4}^{\prime\prime}$ such that for $z=(C_{4}^{\prime\prime}\frac{\sqrt{\log n}}{n^{1/4}}+x)$ we have:

[TABLE]

Recall notation $r_{n}(G)$ introduced in (3). Since $z\geq x$ we have

[TABLE]

Since $r_{n}(G)=O(\frac{1}{\sqrt{n}})$ and $\mathbf{P}(|B_{1}^{n}|>u)\leq 2e^{-u^{2}/2}$, we now have:

[TABLE]

for some constant $\lambda_{4}^{\prime\prime\prime}$. Using (86) and (87) in (5) we obtain that there exist constants $C_{6}$, $K_{6}$ and $\lambda_{6}$ such that:

[TABLE]

∎

We now arrive at our main result for this section, namely a strong embedding for the arrival process $A_{n}$.

Proposition 5.11**.**

Let Assumption 2.1 or 2.5 hold. Then there exists a Brownian motion $\hat{B}$, a Brownian bridge $B^{\textnormal{br},n}$ such that if $\hat{H}_{n}$ be defined as:

[TABLE]

then there exists a version of $T_{1},\ldots,T_{n}$, a version of $\zeta_{1},\ldots,\zeta_{n}$, along with constants $C_{7}$, $K_{7}$ and $\lambda_{7}$ such that

[TABLE]

Proof.

Step 1: Assumption 2.1: From Proposition 5.1, there exists a Brownian motion $\hat{B}$ along with constants $\hat{C}_{1}$, $\hat{K}_{1}$ and $\hat{\lambda}_{1}$ such that for all $n\geq 1$ and $x>0$ we have:

[TABLE]

From Proposition 5.2, there exists a Brownian bridge $B^{\textnormal{br},n}$ along with constants $\hat{C}_{2}$, $\hat{K}_{2}$ and $\hat{\lambda}_{2}$ such that for all $n\geq 1$ and $x>0$ we have:

[TABLE]

Observe that the arrival process $A_{n}$ given by (2) may be decomposed as follows:

[TABLE]

where

[TABLE]

and

[TABLE]

From Corollary 5.6, there exist constants $\hat{C}_{3}$, $\hat{K}_{3}$ and $\hat{\lambda}_{3}$ such that for all $n\geq 1$ and $x>0$:

[TABLE]

Using the bounds (88), (89) and (91) in the decomposition (90), we now obtain existence of constants $C_{7}$, $K_{7}$ and $\lambda_{7}$ such that for all $n\geq 1$ and $x>0$:

[TABLE]

Step 2: Assumption 2.5: From Proposition 5.1, there exists a Brownian motion $\hat{B}$ along with constants $\hat{C}_{1}$, $\hat{K}_{1}$ and $\hat{\lambda}_{1}$ such that for all $n\geq 1$ and $x>0$ we have:

[TABLE]

From Corollary 5.4, there exists a Brownian bridge $B^{\textnormal{br},n}$ along with constants $\hat{C}_{2}$, $\hat{K}_{2}$ and $\hat{\lambda}_{2}$ such that for all $n\geq 1$ and $x>0$ we have:

[TABLE]

Observe that the arrival process $A_{n}$ given by (2) may be decomposed as follows:

[TABLE]

where

[TABLE]

and

[TABLE]

From Corollary 5.6 and Lemma 5.8, there exist constants $\hat{C}_{3}$, $\hat{K}_{3}$ and $\hat{\lambda}_{3}$ such that for all $n\geq 1$ and $x>0$:

[TABLE]

Next from Corollary 5.10, there exist constants $\hat{C}_{4}$, $\hat{K}_{4}$ and $\hat{\lambda}_{4}$ such that for all $n\geq 1$ and $x>0$ we have:

[TABLE]

Ultimately note that according to Assumption 2.5 we have:

[TABLE]

Using the bounds (92), (93), (95), (96) and (97) in the decomposition (94), we now obtain existence of constants $C_{7}$, $K_{7}$ and $\lambda_{7}$ such that for all $n\geq 1$ and $x>0$:

[TABLE]

∎

6. A Strong Embedding for the Workload Process

In this section we derive a strong embedding for the workload process, as well as the total remaining workload.

Proposition 6.1**.**

Let Assumptions 2.1 or 2.5 and 2.7 hold. Then there exist Brownian motions $B$ and $\hat{B}$ and a Brownian bridge $B^{\textnormal{br}}$ such that if $\hat{R}_{n}$ be defined as:

[TABLE]

where $\hat{H}_{n}$ has been defined in Proposition 5.11, then there exists a version of $T_{1},\ldots,T_{n}$, a version of $V_{1},\ldots,V_{n}$, a version of $\zeta_{i},\ldots,\zeta_{n}$ along with constants $C_{8}$, $K_{8}$ and $\lambda_{8}$ such that for all $n\geq 1$ and $x>0$:

[TABLE]

Proof.

From Proposition 5.1, there exists a Brownian motion ${B}$ along with constants $\hat{C}_{1}$, $\hat{K}_{1}$ and $\hat{\lambda}_{1}$ such that for all $n\geq 1$ and $x>0$ we have:

[TABLE]

Observe that the workload process ${W_{n}}$ given by (4) may be decomposed as follows:

[TABLE]

where

[TABLE]

and

[TABLE]

From Proposition 5.11, there exist constants $\hat{C}_{1}$, $\hat{K}_{1}$ and $\hat{\lambda}_{1}$ such that for all $n\geq 1$ and $x>0$:

[TABLE]

From Corollary 5.6, there exists constants $\hat{C}_{3}$, $\hat{K}_{3}$ and $\hat{\lambda}_{3}$ such that for all $n\geq 1$ and $x>0$:

[TABLE]

Using the bounds (98), (100) and (101) in the decomposition (99), we now obtain existence of constants $C_{8}$, $K_{8}$ and $\lambda_{8}$ such that for all $n\geq 1$ and $x>0$:

[TABLE]

∎

The following result helps in extending strong embedding of processes to strong embedding of their Skorohod reflections defined in the sequel.

Lemma 6.2**.**

Let $f$ and $g$ be real-valued functions defined on $[0,\infty)$. Assume $f$ and $g$ satisfy the following property:

[TABLE]

Then we have the following:

[TABLE]

Proof.

Observe that (103) is immediate if the following holds for every $\varepsilon>0$:

[TABLE]

Let us prove (104) through contradiction. First, assume the contrary, namely, there exists $t_{0}\in[0,\infty)$ such that

[TABLE]

Consequently, assume without loss of generality:

[TABLE]

Observe there exist points $t_{1}^{\varepsilon},t_{2}^{\varepsilon}$ such that

[TABLE]

From (102) and (106) we obtain:

[TABLE]

Finally, from (105) and (107) we obtain

[TABLE]

which contradicts the definition of infimum and is not true. Hence our assumption (105) is wrong, and we must have

[TABLE]

Interchanging $f$ and $g$ in (108) allows us to conclude (104), and hence (103) as well. ∎

We now arrive at a strong embedding of the total remaining workload.

Proposition 6.3**.**

Let $\phi$ be the reflection map functional given by:

[TABLE]

Then, under the same assumptions and notations as in Proposition 6.1 there exist constants $C_{9}$, $K_{9}$ and $\lambda_{9}$ such that

[TABLE]

where $\textnormal{id}:x\to x$ denotes the identity function, and $c_{n}$ is a positive constant denoting the server efficiency rate.

Proof.

Observe from Proposition 6.1, we have:

[TABLE]

From Lemma 6.2 we have

[TABLE]

Combining equations (109) and (6), and recalling the definition of $\phi$ we obtain:

[TABLE]

for some constants $C_{9}$, $K_{9}$ and $\lambda_{9}$. ∎

7. A Strong Embedding for the Queue Length Process

In this section we obtain a strong approximation to the queue length process. Control of the truncated renewal process $M_{n}$ (recall relation 5) would lead to a strong approximation of the queue length.

Lemma 7.1**.**

Let $Z_{n}(t)$ be given by:

[TABLE]

Then there exist constants $C_{10}$, $K_{10}$ and $\lambda_{10}$ such that

[TABLE]

Proof.

Notice $Z_{n}(t)$ can be further decomposed as:

[TABLE]

Henceforth, the two terms in (112) will be approximated. Utilizing the definition of $M_{n}(t)$, the second term in (112) can be bounded as follows:

[TABLE]

Hence for any constant $C>0$ we obtain:

[TABLE]

Using Chebyshev’s inequality we obtain:

[TABLE]

for some constants $K$ and $\lambda$, where the last step is obtained by choosing $\delta$ sufficiently small.

In order to approximate the first term in (112), observe that from Proposition 5.1 there exist constants $C_{1}$, $K_{1}$ and $\lambda_{1}$ such that:

[TABLE]

which implies

[TABLE]

Using (113) and (114), our desired inequality (111) is obtained. ∎

We now approximate the Brownian motion evaluated at $M_{n}$ appearing in Lemma 7.1

Lemma 7.2**.**

There exist constants $C_{11}$, $K_{11}$ and $\lambda_{11}$ such that

[TABLE]

Proof.

It suffices to check the conditions in Proposition 4.12, which are satisfied by Proposition 4.7. Combining equations (26), (27) and (28) we obtain:

[TABLE]

Using the fact that $\mathbf{E}(\sup_{0\leq k\leq n}|\tilde{S}_{k}-k|)\leq C\mathbf{E}|\tilde{S}_{n}-n|$, we obtain:

[TABLE]

This yields our desired result. ∎

Remark 7.3*.*

Notice from definition, $M_{n}(t)$ equals $n$ for all $t>\frac{S_{n}}{c_{n}}$. Hence, control of $B_{\frac{M_{n}(t)}{n}}$ for $t\in[0,\infty)$ reduces to a control of $B_{\frac{M_{n}(t)}{n}}$ for $t\in[0,\frac{S_{n}}{c_{n}}]$. However Lemma 7.2 leads a control of $B_{\frac{M_{n}(t)}{n}}$ over $t\in[0,L_{n}]$ for a predetermined and fixed sequence $L_{n}$. Hence we need $\frac{S_{n}}{c_{n}}$ to be in an interval $[0,L_{n}]$ with exponentially high probability (i.e. the complement event has exponentially decreasing probability). This is achieved in the following lemma.

Lemma 7.4**.**

Let $S_{n}=V_{1}+\ldots+V_{n}$. Then for every $\eta>0$, there exists $\delta>0$ such that

[TABLE]

Proof.

[TABLE]

for some $\delta$ small enough such that $\mathbf{E}e^{\delta V_{1}}<e^{\eta}$. ∎

From Lemma 7.2 and Remark 7.3 we obtain the following result, which controls $B_{\frac{M_{n}(t)}{n}}$ for all $t$ positive.

Lemma 7.5**.**

[TABLE]

Finally we arrive at the main result for this section, namely a strong embedding for the queue length $Q_{n}$.

Proposition 7.6**.**

Let Assumptions 2.1 or 2.5 and 2.7 hold. Let $\tilde{E}_{n}$ be given by:

[TABLE]

Recall $\hat{H}_{n}$ defined in Proposition 5.11 and define the process $\hat{Y}_{n}$ as follows:

[TABLE]

Then there exists a version of $T_{1},\ldots,T_{n}$, a version of $V_{1},\ldots,V_{n}$, a version of $\zeta_{i},\ldots,\zeta_{n}$ along with constants $C_{12},K_{12}$, $\lambda_{12}$ and $\zeta$ independent of $n$, such that for all $n\geq 1$ and $x>0$:

[TABLE]

if $c_{n}=O(n^{m})$ for some $m>0$ and $\liminf_{n}c_{n}>0$, else:

[TABLE]

Proof.

Observe that from (6) the diffusion-scaled queue length $\frac{Q_{n}}{\sqrt{n}}$ can be further decomposed as:

[TABLE]

where the idle time process $I_{n}$ has been defined in (7), $Y_{n}$ is given by:

[TABLE]

and $\hat{H}_{n}$ has been defined in Proposition 5.11. By the Skorohod reflection theorem we have that

[TABLE]

and the busy time process is given by:

[TABLE]

The diffusion-scaled queue length process is now given by the Skorohod reflection of $Y_{n}$:

[TABLE]

Thus a strong embedding of $\frac{Q_{n}}{\sqrt{n}}$ would follow from a strong embedding of $Y_{n}$. Notice that we already have a strong embedding of the arrival process courtesy Proposition 5.11, namely there exist constants $C_{7}$, $K_{7}$ and $\lambda_{7}$ such that for all $n\geq 1$ and $x>0$:

[TABLE]

Thus to complete the strong embedding of $Y_{n}$, we need to approximate $\tilde{Z}_{n}$ given by:

[TABLE]

Observe that the busy time process $D_{n}(t)$ is non-decreasing and takes values in $[0,\frac{S_{n}}{c_{n}}]$ where $S_{n}=\sum_{i=1}^{n}V_{i}$. Consequently from Lemma 7.1, there exist constants $\hat{C}_{1}$, $\hat{K}_{1}$ and $\hat{\lambda}_{1}$ such that for all $n\geq 1$ and $x>0$:

[TABLE]

Let

[TABLE]

Using (119) and (120) in (116) we thus have constants $\hat{C}_{2}$, $\hat{K}_{2}$ and $\lambda_{2}$ such that for all $n\geq 1$ and $x>0$:

[TABLE]

From the converse of Lemma 6.2, there exist constants $\hat{C}_{3}$, $\hat{K}_{3}$ and $\hat{\lambda}_{3}$ such that for all $n\geq 1$ and $x>0$:

[TABLE]

Recalling the expression (7) and using (117) we have for every $n\geq 1$ and $x>0$:

[TABLE]

Consequently for every $n\geq 1$ and $x>0$:

[TABLE]

Recall $\tilde{E}_{n}$ given by (115). From the expression of $\tilde{Y}_{n}$ in (121) and recalling $\hat{H}_{n}$ from Proposition 5.11 we have:

[TABLE]

Observe that both $G(t)$ and $\frac{M_{n}(D_{n}(t))}{n}$ are less than $1$. Using the tail probability for the supremum of the standard Brownian bridge on $[0,1]$ we have:

[TABLE]

Using the tail probability for the supremum of the Brownian motion on $[0,1]$ we have

[TABLE]

and

[TABLE]

Using (125), (126) and (127) in (7), we have that there exist constants $k_{1}$ and $k_{2}$ such that:

[TABLE]

From the converse of Lemma 6.2, there exist constants $\hat{k}_{1}$ and $\hat{k}_{2}$ such that:

[TABLE]

Now, using (123) and (128) we have:

[TABLE]

Since $|x\wedge 1-y\wedge 1|\leq|x-y|$ we have from (129) constants $\hat{k}_{3}$, $\hat{k}_{4}$ and $\hat{k}_{5}$ such that:

[TABLE]

From Proposition 4.7 we have constants $\hat{k}_{6}$ and $\hat{k}_{7}$ such that:

[TABLE]

Combining (130) and (131) we have constants $\tilde{k}_{0}$, $\tilde{k}_{1}$, $\tilde{k}_{2}$ and $\tilde{k}_{3}$ such that:

[TABLE]

Now, from Proposition 4.12 we obtain constants $\hat{C}_{4}$, $\hat{K}_{4}$ and $\hat{\lambda}_{4}$ such that:

[TABLE]

Let $\hat{Y}_{n}$ be given by

[TABLE]

We now obtain from (121) and (132) existence of constants $\hat{C}_{5}$, $\hat{K}_{5}$ and $\hat{\lambda}_{5}$ such that :

[TABLE]

Hence from (122) we have constants $\hat{C}_{6}$, $\hat{K}_{6}$ and $\hat{\lambda}_{6}$ such that

[TABLE]

Recalling (118), we now obtain that there exist constants $\hat{C}_{7}$, $\hat{K}_{7}$ and $\hat{\lambda}_{7}$ such that

[TABLE]

Observe that for any two sets $A$ and $B$, we have that $P(A)\leq P(A|B)+P(B^{c})$. Thus using Lemma 7.4 we have:

[TABLE]

Finally choosing $L_{n}=n$ yields the desired result. ∎

8. Commentary and Conclusions

From a philosophy of science perspective, one can consider the bulk of the nonstationary queueing model literature as phenomenological ([FrHa2006]) in nature – that is, accurately reflecting empirical evidence, but not necessarily first principles. For instance, as noted in Remark 2.7 (see [Wh2018] as well) a widely used nonstationary traffic model uses a composition construction, where a cumulative intensity function that captures the time-varying effects is posited. However, these models are not necessarily a first principles explanation of how customers choose to arrive at a service system. This distinction between phenomenological and mechanistic modeling is not crucial from a performance analysis/prediction perspective, but it can be important from a system design or optimization and control perspective. For instance, in [ArAtHo2017] the authors study the design of an optimal appointment schedule to a nonstationary single server queueing system. In this instance, the composition constructed time inhomogeneous queueing models [Wh2018] are not appropriate and the authors used a transitory queueing model. On the other hand, stationary queueing models111we define this as the content in either [gross] or [kleinrock], on the other hand, are also mechanistic models [FrHa2006] in nature whereby the arrival and service models are directly related to the behavior and choices of individual customers. The  $\textnormal{RS}(G,p)/G/1$ model provides a mechanistic description of queueing behavior and can be used in optimization, control and game theoretic models of queueing behavior with much flexibility.

The  $\textnormal{RS}(G,p)/G/1$ model generalizes the  $\Delta_{(i)}/G/1$ model that has been studied in the literature. However, computing performance metrics for the discrete event  $\textnormal{RS}(G,p)/G/1$ queueing model is quite difficult, owing to the complicated time-dependencies in the model. The strong embeddings (and FSATs) proved in this paper provide error bounds from tractable diffusion approximations. Furthermore, these results can be specialized to yield prior diffusion limits obtained via weak convergence in [HoJaWa2016a, HoJaWa2016b, BeHoLe2019]. We anticipate that our FSAT results to be immensely useful for optimization and control problems involving queueing systems.

There are several avenues for further exploration. First, our most general conditions on the traffic model in Assumption 2.5 allows the arrival epoch distribution to depend on the population size – allowing for the possibility that an increase in the population will change customer behavior since there is an increase in demand for services. We currently assume that the “dropout” probability is stationary. A more general model would allow for time-dependent/non-stationary dropout probabilities. It seems possible to extend the current FSETs and FSATs to this setting. More complicated is establishing analogous results for a multi-server queue. In [HoJaWa2016a], diffusion limits were established for a fixed multi-server queue in the large population asymptotic limit, relying on the fact that in the large sample limit the regulator process is identical to the single server case. In our current paper, however, the FSETs (which are for finite $n$) are much harder to prove, since we can no longer use the asymptotic simplification. This issue is compounded when the servers are not identical, and we will investigate these results in future papers. A further avenue for investigation is how to prove FSETs in a scaling regime that is analogous to the many-server heavy-traffic (MSHT) scaling. In this case, we anticipate that the diffusion approximation should be some type of a non-stationary Halfin-Whitt diffusion process, but we have not been able to prove the FSET, and it appears we might require some new mathematical innovations to achieve this result.

Appendix

Theorem 8.1**.**

Let $F$ be a distribution function with mean [math] and variance $1$. In addition, suppose the moment generating function corresponding to $F$, $R(t)=\mathbb{E}(e^{tX})$, $X\sim F$, exists in a neighborhood of [math]. Then, given a Brownian motion $B$, and using it, one can construct a sequence of random variables $X_{1},X_{2},\ldots$ which are independent and identically distributed to $F$. Furthermore, the partial sums of $X_{i}$’s are strongly coupled to the Brownian motion $B$ in the following sense. For every $n\in\mathbb{N}$ and $x>0$:

[TABLE]

where $C$, $K$ and $\lambda$ are positive constants depending only on $F$.

Proof ideas and construction:.

Let us first give a sketch of the construction that yields the random variables $\mathbf{X}:=\{X_{i},i\in\mathbb{N}\}$. Instead of fixing the Brownian motion it would suffice to fix an infinite sequence of independent standard normals $\mathbf{Y}:=\{Y_{i},i\in\mathbb{N}\}$. Let $G$ be the distribution function for a generic element from $\mathbf{Y}$ and recall $F$ is the distribution function for a generic element from $\mathbf{X}$.

Observe that $F(X)$ and $G(Y)$ are both $U[0,1]$ distributed random variables, and a natural suggestion to generate $X$ given $Y$ (or vice versa) is to equate them:

[TABLE]

assuming $F^{-1}$ is properly defined. Equation (134) is the quantile transform which forms the backbone of the construction, details of which are to follow.

The above technique turns out to be useful especially if one is interested in bounding $|X-Y|$ as well (which is our aim too), e.g., when $F$ is the distribution of sum of $n$ iid mean [math] random variables with finite moment generating function, and $G$ is the distribution of appropriately normalized normal random variable, then

[TABLE]

Note the error of approximation does not increase with the number of components in $X$.

If someone wants to approximate $S_{n}=X_{1}+\cdots+X_{n}$ by $T_{n}=Y_{1}+\cdots+Y_{n}$ (we are interested in $Y_{i}$’s distributed as normals), such that the components are independent, then we may divide $S_{n}$ into blocks, in particular if:

[TABLE]

we may approximate the block $X_{n_{i}}+\cdots+X_{n_{i+1}}$ using $Y_{n_{i}}+\cdots Y_{n_{i+1}}$ whereby this approximation is first achieved by quantile transformation of $Y_{n_{i}}+\cdots+Y_{n_{i+1}}$ and then finding the individual $X_{i}$’s by conditioning on the sum $X_{n_{i}}+\cdots+X_{n_{i+1}}$ (details to follow). This is the key idea in the KMT construction.

Suppose $(X_{m+1}+\cdots+X_{m+2n})$ has already been obtained as quantile transform of $(Y_{m+1}+\cdots+Y_{m+2n})$. For normally distributed $Y_{i}$’s we have $[(Y_{m+1}+\cdots+Y_{m+n})+(Y_{m+n+1}+\cdots+Y_{m+2n})]$ independent of $[(Y_{m+1}+\cdots+Y_{m+n})-(Y_{m+n+1}+\cdots+Y_{m+2n})]$. Thus one may expect:

[TABLE]

to be approximately independent of

[TABLE]

Also let

[TABLE]

and

[TABLE]

As a consequence, one can obtain two blocks out of $U_{2n}$, having already obtained $U_{2n}$ from quantile transform of $V_{2n}$:

1. (i)

Consider $F(x|y)=P(\tilde{U}_{2n}<x|U_{2n}=y)$. 2. (ii)

Transform $G(V_{2n}/\sqrt{2n})$ which is $U[0,1]$ distributed by the application of the quantile transform, i.e.,

[TABLE] 3. (iii)

Now, let

[TABLE]

and

[TABLE]

Observe that $(U_{2n}+\tilde{U}_{2n})/2$ and $(U_{2n}-\tilde{U}_{2n})/2$ are independent.

In the construction these steps are applied iteratively, until the individual random variables get constructed. To this end, the following structure/blocking is considered in the beginning: the set of standard normals $\mathbf{Y}$ is broken down into: $\{Y_{1},Y_{2}\}\cup_{j=1}^{\infty}\{Y_{2^{j}+1},\ldots,Y_{2^{j+1}}\}$. From each of these blocks we construct $\{X_{1},X_{2}\}\cup_{j=1}^{\infty}\{X_{2^{j}+1},\ldots,X_{2^{j+1}}\}$ and then apply the three steps outlined above repeatedly until we arrive at individual random elements. This is possible because the block sizes are powers of $2$.

The random variables $X_{i}$’s thus generated are clearly independent of each other. To see this, observe that the blocks $\{Y_{2^{j}+1},\ldots,Y_{2^{j+1}}\}$ are independent and hence the blocks $\{X_{2^{j}+1},\ldots,X_{2^{j+1}}\}$ are independent as well. Each division of the block also generate independent random variables. Thus the newly constructed $\mathbf{X}$ is composed of independent entries. Furthermore each of them is distributed with distributed with distribution function $F$.

The only remaining argument is to show that the $X_{i}$’s thus generated are close to the $Y_{i}$’s. For some $m\in\mathbb{N}$, let $j=\max\{i:2^{i}|m\}$. Then $m=(2k+1)2^{j}$ for some $k\in\mathbb{N}$. One can write $S_{m}$ as follows:

[TABLE]

or written according to the notations in [KMT]:

[TABLE]

Applying this recursion repeatedly for $2^{n}<m\leq 2^{n+1}$ one arrives at:

[TABLE]

where $\tilde{S}_{m}$ is the linear interpolation between $S_{2^{n}}$ and $S_{2^{n+1}}$ given by:

[TABLE]

and $c(i)$ and indices $k(i)$ depend on $m$, $0\leq c(i)\leq 1$, $k(i)=\lceil\frac{m}{2^{i}}\rceil$.

A similar representation is possible for $T_{n}=Y_{1}+\cdots+Y_{n}$, $n\in\mathbb{N}$. One obtains:

[TABLE]

Lemma 1 in [KMT] gives bounds on these differences. These are used in a string of inequalities to bound the probability:

[TABLE]

Existence of moment generating function is required in order to apply Chernoff inequality to obtain the best possible bound. Note that non-existence of moment generating function will not provide an exponentially decreasing upper bound to (135). ∎

Theorem 8.2**.**

There exists a probability space with independent $U[0,1]$ random variables $U_{1},U_{2},\ldots$ and a sequence of Brownian bridges $B_{1}^{\textnormal{br}},B_{2}^{\textnormal{br}},\ldots$ such that for all $n\geq 1$ and $x\in\mathbb{R}$:

[TABLE]

where the empirical process $\alpha_{n}$ is given by

[TABLE]

and

[TABLE]

Proof ideas and construction:.

The following has been borrowed extensively from [mattner] and the reader is encouraged to look at Sections 25.1-25.3 for a clearer exposition. Here as well, we would construct the random variables $U_{1},\ldots,U_{n}$ (for a fixed $n$) from the Brownian bridge $B_{n}^{\textnormal{br}}$. In fact, a similar dyadic scheme is used to construct the required variables.

The first key idea in the construction is the fact that every real number in $[0,1]$ has a binary representation (terminating for rationals, while non-terminating for irrationals).

1. (i)

Let $\mathscr{Z}=\{Z\}\cup_{j=1}^{\infty}\{Z_{i}:i\in{\{0,1\}}^{j}\}$ be an indexed set of independent standard normal random variables. 2. (ii)

Let $H_{m}$ be the inverse distribution function for $\textnormal{Bin}(m,\frac{1}{2})$. Define $H_{0}=0$. 3. (iii)

Let

[TABLE]

Next for each $i\in{\{0,1\}}^{j}$, $j\geq 1$:

[TABLE]

Observe that $n=N_{0}+N_{1}$, and for each $i\in{\{0,1\}}^{j}$, $j\geq 1$ we have $N_{i}=N_{i,0}+N_{i,1}$. This implies for each $j\geq 1$, $N_{i},i\in{\{0,1\}}^{j}$ is $\textnormal{Multinomial}(n,\frac{1}{2^{j}}\cdots\frac{1}{2^{j}})$.

We have obtained $\{N_{i},i\in{\{0,1\}}^{j},j\geq 1\}$ from standard normals. Using these $N_{i}$’s we will now construct a uniform empirical distribution function $G_{n}$ as follows: For $j\geq 1$ and $k=1,\ldots,2^{j}$ let

[TABLE]

where

[TABLE]

With increasing $j$, the $G_{n}^{j}$’s become finer and finer. Thus by taking limits (as rationals are dense in [0,1]) we can construct a uniform empirical distribution function defined on the entire interval $[0,1]$:

[TABLE]

Inverting $G_{n}$ one obtains the order statistics:

[TABLE]

of $n$ independent $U[0,1]$ random variables. A random permutation of these order statistics provides a sample of $n$ independent $U[0,1]$ random variables.

Our aim is to construct $\alpha_{n}(s)=\sqrt{n}(G_{n}(s)-s)$, such that it is close to the Brownian bridge, $B^{\textnormal{br}}$. This means that the set of standard normals we considered previously should be obtained from the Brownian bridge in a nice manner. The following is the nice manner (here we have denoted the $n^{\text{th}}$ Brownian bridge as $B$):

1. (i)

$Z=2B\left(\frac{1}{2}\right)$. 2. (ii)

For $i\in{\{0,1\}}^{j}$, $j\geq 1$:

[TABLE]

where $k$ is obtained from $\frac{k}{2^{j}}=\sum_{s=1}^{j}\frac{i_{s}}{2^{s}}$ $(i=(i_{1},\ldots,i_{j}))$. 3. (iii)

It can be checked that these $Z_{i}$’s are standard normals and independent of each other.

Next one needs to argue why this construction of the uniform random variables produces an empirical process which is close to the Brownian bridge in the sense of (16). Some sketch of the proof can be found in [Mason-notes]. ∎

Notice the following remark already alluded to in the main text.

Remark 8.3*.*

The above construction of uniform random variables is for a fixed $n$. Having obtained $\{U_{1},\ldots,U_{n}\}$, one is unable to obtain another $U_{n+1}$ such that the new set $\{U_{1},\ldots,U_{n+1}\}$ satisfies (16) with the same Brownian bridge. Instead we have to redo our construction. This necessitates the need for a different Brownian bridge $B^{\textnormal{br},n}$ for every $n$. However the following statement is true:

There exists a Brownian bridge $B$ such that for each $n$, there exists $n$ independent uniforms $U_{1},\ldots,U_{n}$ whose empirical process $\alpha_{n}$ satisfies:*

[TABLE]