Uncertainty in modeling overdispersed queues
-An analysis of scaling limits for infinite-server queues in a random environment. Mariska Heemskerk
Supervisors: prof. dr. Michel Mandjes & prof. dr. Johan van Leeuwaarden June 25, 2015
Korteweg-de Vries Institute for Mathematics Faculty of Science
Contents
1 Introduction 4
2 Overdispersion in an infinite-server context 9
2.1 Pre-limit results . . . 9 2.2 Limit results . . . 13
3 Model with correlated arrival streams: functional central limit theorem 17 3.1 Qualitative behavior of first two moments . . . 18 3.2 Functional central limit theorem . . . 20
4 Large deviations 26
4.1 Univariate large deviations . . . 26 4.2 Large deviations for the coupled model . . . 32
5 Conclusions and outlook 33
5.1 Alternative models using mixed Poisson arrivals . . . 33 5.2 Correlated queues . . . 34 5.3 General service times . . . 34
Abstract
Large-scale service systems are typically modeled using a multi-server system with a nonho-mogeneous Poisson arrival process (NHPP). Empirical studies provide statistical evidence that actual arrival processes in call centers and healthcare often fluctuate more than the NHPP-assumption can capture. This calls for arrival processes and hence queueing systems that model overdispersion. To this end, we introduce a mixed Poisson arrival process in a random environment. Studying this arrival process in the context of an infinite-server queue, we establish multiple stochastic-process limits through asymptotic analysis. This gives rise to a trichotomy in system behavior, in which different regimes arise depending on a delicate interplay between the overdispersion caused by the mixed Poisson distribution and the rate at which the random environment changes.
1
Introduction
Queueing theory was developed in order to solve capacity-sizing problems for service systems [7], a type of stochastic networks consisting of a combination of installed work-processing resources and typically a large number, say N , of work-generating sources. Each source generates arrivals at a certain arrival rate. Examples of service system usage include modeling service operations in call centers and healthcare. The capacity-sizing problem is to determine the optimal capacity to guarantee a certain desired system performance (here the capacity is defined as the number of servers). The required capacity typically increases with the level of uncertainty in the system. In service system modeling, as in stochastic modeling in general, one can distinguish three types of uncertainty [14]. First, there is uncertainty about the choice of the probabilistic model that is used to approximate reality. This is referred to as model uncertainty. Given that the model is appropriate, there can still be uncertainty concerning the parameters of the model: parameter uncertainty. Finally, of course, once the model and parameters are set, intrinsic stochastic variability provides for uncertainty of the outcomes of a process.
Taking into account the economies-of-scale principle, service systems are typically modeled using many-server queues, where the canonical models are M/M/s and M/M/∞ systems with Poisson arrivals, exponentially distributed service times and s or infinitely many servers, respectively. In this thesis, we consider infinite-server queues, as this will make the analysis of our rather complicated setting tractable. As observed by Whitt [13], the infinite-server assumption is natural when the goal to achieve is to immediately serve all jobs. Besides, using techniques as PSA and MOL, first-order insights in system behavior in a finite-server setting can be obtained afterwards [15].
In the majority of studies, it is assumed that the arrival rate of the combined work-generating sources, the so-called total arrival ratePN
i=1λi= N λ, is a known parameter that can be extracted
from the data as the mean number of arrivals per time unit [2]. Consequently, the arrivals of the N sources are considered as one single process, the arrival process, the arrival rate of which equals the total arrival rate. Given an M/M/∞ system, the assumption of independent sources (the individual arrival processes are independent processes) results in a Poisson arrival process with mean and variance both equal to N λ. However, recent studies suggest that this independence assumption seems invalid: there is growing evidence that in many practical cases sources are highly correlated [2]. For instance, it is observed that the number of arrivals fluctuate more than we would expect from a Poisson process [12, 10, 4]. Given a non-negative random variable X, the index of dispersion or the variance-to-mean ratio is
VMRX = Var X
EX .
Under the standard Poisson assumption the index of dispersion of the number of arrivals in a time unit equals 1, whereas in other models it can be larger (smaller) than 1; this phenomenon is referred to as overdispersion (underdispersion).
Overdispersion was observed in all kinds of data sets, for instance in call center data [4] and hospital arrivals [9]. As this violates the Poisson assumption, it renders results as the widely used square-root safety capacity prescription deficient for implementation. Indeed, Zeevi et al. [2] showed that this capacity rule performs poorly when handling overdispersed data. However, they observe that in practice as well as in literature, forecast errors due to inadequate capacity prescriptions are often ascribed to stochastic variability. Of course, in the event of significant overdispersion, such malfunctioning should be explained by another type of uncertainty instead. Only recently, researchers have begun to recognize the urge to add extra built-in uncertainty to service system models. As first suggested by Whitt [13], overdispersion can be explained by parameter uncertainty, which actually could be viewed (but not necessarily) as a manifestation of correlation between sources. The simplest way to implement parameter uncertainty is using mixed Poisson arrivals; for related work, see e.g. [8, 10, 11] and references therein.
B Parameter uncertainty as a manifestation of correlated sources
When modeling an arrival process with overdispersed arrivals from a many-sources perspective while holding on to the paradigm of an M/M/∞ system, it can be concluded from the above that we should implement correlation between sources. This can be done via a random arrival rate: instead of λ, the sources all have random arrival rates sampled from the distribution of some random variable Λ, such that an underlying random mechanism determines the rates for all sources. In more general language: parameter uncertainty causes extra variability in the arrival process, rendering it a useful implementation when modeling dispersed data.
To understand the relation between parameter uncertainty and correlated sources, note that correlation can be interpreted as the effect of an underlying random mechanism that influences multiple processes. As before, suppose that there is one arrival rate for all sources (λi = λ for
i ∈ {1, . . . , N }), which is called the state of the total arrival process. The idea behind a random arrival rate is that λ is replaced by a random variable Λ such that EΛ = λ: an underlying random mechanism determines the state.
In previous work, overdispersion is often modeled using so-called mixed Poisson arrivals. As we will show, this construction leads to a level of overdispersion proportional to the size of the system and the variance-to-mean ratio of Λ. Consider the following extension of a Poisson process.
Definition 1.1. Let Λ : Ω1 → R a non-negative random variable, absolutely continuous with
density fΛ(λ). Then XΛ = (XtΛ)t∈T is a mixed Poisson process if conditional on the value of
Λ, XΛ is a Poisson process. Then XtΛ ∼ mPoisson(Λt), i.e. the distribution of XΛ
t is mixed
Poisson and defined via:
P(XtΛ= n) = EΛ (Λt)n n! e −Λt =Z ∞ 0 (λt)n n! e −λt fΛ(λ) dλ.
Let PiΛ denote the mixed Poisson arrival process for source i, with random parameter Λ. When the PiΛ all respond to the same realization of Λ, the total arrival process P :=PN
i=1PiΛis mixed
Poisson with parameter N Λ, hence the process is dispersed. Note that the PiΛ are correlated through their parameter
From the system’s perspective, the mean number of mixed Poisson arrivals with random rate N Λ in 1 time unit is N λ and its variance is N λ + N2Var(Λ). Hence for large values of N the variance is huge compared the mean, and the VMR equals
1 + NVar(Λ)
λ . (1)
This VMR is tunable via Var(Λ), which is an important first step in fitting a model to actual overdispersed data. We next enrich this mixed Poisson setting with a random environment in order to establish a more flexible, and more realistic formalism for actual arrival processes.
B Channeling uncertainty via random environment
Mixed Poisson arrivals are overdispersed simply due to the arrival rate being randomly drawn from the distribution of Λ, thus considerably generalizing the case where the rate is a deterministic value EΛ = λ. In order to channel parameter uncertainty, we introduce a sampling frequency: not only is the arrival rate a random unknown value, it also keeps changing over time at a certain pace. The classical mixed Poisson setting is then enriched by extending the random rate to a random environment, where a new arrival rate is independently drawn from Λ with frequency ∆1, for some ∆ > 0. One can imagine that in a rapidly changing random environment, the randomness of Λ does not play a role anymore, whereas in a slowly changing random environment, the variability of Λ has a lot of influence on the variance of the arrival process. This trade-off between the influence of the sampling frequency and parameter uncertainty is a convenient and interesting device to have in a model. It could be used to tune the model in such a way that the resulting level of overdispersion fits the data.
Number of arrivals ( ) Time ( ) Δ 2Δ 3Δ 4Δ 0
Figure 1 shows a realization of a mixed Poisson arrival process in a random environment. The number of arrivals per time unit is locally Poisson and fluctuates proportionally to its local parameter, independently drawn from the distribution of the random parameter Λ with EΛ = λ. As the paramater changes every ∆ time units in this random environment, the number of arrivals per time unit fluctuates more severely than normal Poisson data would (e.g. the sample variance is larger than the sample mean, which is approximately EΛ = λ).
In order to explain how we extend the standard mixed Poisson construction with a random environment, let us now define the precise arrival process to be used in our model. Consider a Markovian arrival process with a random arrival rate that is newly sampled after each ∆ time units, i.e. at deterministic sampling times reoccuring at frequency ∆1. Assuming N sources, the arrival process is an instance of a Cox process [5] with time-dependent arrival rate N Λ∆(s), defined as
Λ∆(s) :=X
j
Λj1[j∆,(j+1)∆)(s), (2)
with i.i.d. Λjdistributed as Λ, absolutely continuous with density fΛ(λ) and finite moments. The
number of arrivals after t time units, denoted by Pt, is mixed Poisson with random parameter
NR0tΛ∆(s) ds, hence if R0tΛ∆(s) ds =Pt/∆ j=1λj∆ (suppose t/∆ ∈ N), then P(Pt= n) = Z ∞ 0 · · · Z ∞ 0 (NR0tΛ∆(s) ds)n n! e −NRt 0Λ ∆(s) ds t/∆ Y j=1 fΛ(λj) dλ1. . . dλt/∆.
Note that on itself, introducing a sampling frequency does not lead to interesting results. For instance, the mean number of arrivals per time slot is still N λ and the variance is again N λ + N2Var(Λ); even when ∆ < 1 (suppose ∆1 ∈ N),
Var(P1) = 1 ∆ · Var(P∆) + 1 ∆( 1 ∆− 1)Cov(P 1 ∆, P∆2) = 1 ∆(N λ∆ + N 2∆2Var(Λ)) + 1 ∆( 1 ∆− 1)N 2∆2Var(Λ)
= N λ + N2∆Var(Λ) + N2Var(Λ) − N2∆Var(Λ) = N λ + N2Var(Λ),
where P∆1 and P∆2 are i.i.d. mixed Poisson random variables with random parameter N Λ∆. Therefore, the VMR of the number of arrivals per time slot is the same as in (1). Nevertheless, imposing a scaling on the sampling frequency, the asymptotic behavior of the resulting system is interesting.
B Asymptotic analysis results
We have thus arrived at a single-class Markovian infinite-server queue with mixed Poisson arrival governed by a random environment, so that the random arrival rate is refreshed at fixed points in time. Surprisingly, despite the rather complicated set-up, it turns out to be possible to
compute the moments and the probability generating function of the queue length process. Also, asymptotic analysis of the queue length under space and time scaling is tractable and will expose powerful results, e.g. (functional) central limit results and results concerning rare events. The asymptotic analysis is not only tractable, it can also be motivated as follows. Following the economies-of-scale principle, we consider large systems (N large), rendering the sampling frequency ∆1 small compared to N Var(Λ), the variance of the random arrival rate. The battle of the sampling frequency against the variance of Λ is won before it’s fought; the variance will have a lot of influence on the system’s behavior. In the limit, as N tends to infinity, the frequency will get infinitely small compared to the system size, which will of course not lead to interesting behavior. On the contrary, in the event of a relatively high sampling frequency the variance of Λ has no impact on the system’s behavior. Therefore, we want to control the rate at which the frequency increases as the system size grows. This is accomplished by scaling ∆1 to N∆α for some α > 0, where for α < 1 the variance eventually has the upper hand, whereas for α > 1, the sampling frequency grows faster than the variance of Λ. Choosing α = 1, hence letting the system size and the sampling frequency grow at the same rate, we find a third regime that combines the two others. This trichotomy plays an important role in the analysis and will continue to pop up throughout the different sections.
B Related work
This work is inspired by analyses of standard mixed Poisson arrival processes used to model overdispersion in e.g. [2, 8, 9, 10, 11]. Here, it should be noted that in these studies, dispersed data over disjoint time slots are modeled with individual, independent mixed Poisson arrival processes, whereas we use a continuous dispersed arrival process to model at once data from, say, a day.
In a way, the random environment adaptation of the standard mixed Poisson process resembles Markov-modulation in arrival processes, where the process has rate λi whenever the background
process, a continuous-time Markov chain, is in state i. In this thesis, we will follow joint work of Anderson, Blom and de Turck on Markov-modulated infinite-server queues [1, 3, 6]. Whereas we scale the sampling frequency, in this setting the jump rate corresponding to the background Markov process is scaled (sub/super-)linearly as the system size tends to infinity. This leads to interesting asymptotic behavior; see [1] for the derivation of a functional central limit theorem related to a generalized OU process, and [3, 6] for a rare-event analysis.
2
Overdispersion in an infinite-server context
In this section we present a transient and stationary analysis of a single-class Markovian infinite-server queue with random environment of sample frequency ∆1 given by Λ(s) = Λ∆(s) as in (2). Remember that we assumed that the sampled arrival rates are i.i.d. and distributed as a random variable Λ > 0 with finite first two moments and density fΛ(·). The corresponding service times
are assumed i.i.d. (and in addition independent of the arrival process), exponentially distributed random variables with mean 1/µ; some results carry over to the case of general i.i.d. service times, as we indicate below. We particularly focus on asymptotic normality under a specific parameter scaling.
2.1 Pre-limit results
In this first subsection we analyze the queue’s stationary behavior, in terms of the Laplace trans-form and the corresponding moments. We then extend these findings to the associated transient behavior. This subsection presents ‘pre-limit results’; later we impose parameter scalings which facilitate the derivation of a central limit theorem.
B Transform of stationary behavior
Let M be the random variable associated with the stationary number of customers (also some-times referred to as ‘jobs’) in the system. We describe this random variable in terms of its pgf φ(z) := EzM.
In the sequel we write pt := e−µt for the probability that a job present at kt is still present at
(k + 1)t and qt:= (1 − e−µt)/(µt) for the probability that a job arriving at a uniform epoch in
[kt, (k + 1)t) is still present at (k + 1)t. Denote ¯pt:= 1 − pt and ¯qt:= 1 − qt.
It is not hard to see that, owing to the fact that we observe the system in stationarity,
φ(z) = ∞ X k=0 P(M = k) k X `=0 z`k ` p`∆(1 − p∆)k−` ! × ∞ X k=0 Z ∞ 0 fΛ(λ)e−λ∆ (λ∆)k k! k X `=0 z`k ` q∆` (1 − q∆)k−`dλ ! = φ(zp∆+ ¯p∆)gΛ,∆(z), (3)
where gΛ,t(z) is defined as follows, with rt:= tqt:
gΛ,t(z) :=
Z ∞
0
exp (−λrt(1 − z)) fΛ(λ) dλ = E exp (−Λrt(1 − z)) .
Observe that gΛ,t(z) is a pgf of ‘mixed Poisson’ type: conditional on Λ = λ the pgf corresponds
An iteration argument yields φ(z) = ∞ Y k=0 gΛ,∆ zpk∆+ k−1 X j=0 ¯ p∆pj∆ = ∞ Y k=0 gΛ,∆ 1 − (1 − z)pk∆. (4)
This expression has an insightful interpretation. Note that the stationary number in the system, M , can be written as the sum of M0, M1, M2, . . ., where Mk represents the number of jobs
that arrived in [−(k + 1)∆, −k∆) that are still present at time 0. Furthermore, observe that the Mk are independent. It is easy to see that the probability of an arbitrary arriving job in
[−(k + 1)∆, −k∆) (i.e., having arrived at a uniform epoch in this interval) is still in the system at time 0 with probability
Z ∆ 0 1 ∆e −µ(k∆+s)ds = 1 µ∆ e−µk∆− e−µ(k+1)∆= p k ∆− p k+1 ∆ µ∆ = pk∆p¯∆ µ∆ . (5) As a consequence, EzMk = ∞ X `=0 Z ∞ 0 fΛ(λ)e−λ∆ (λ∆)` `! ` X m=0 zm ` m pk ∆p¯∆ µ∆ m 1 −p k ∆p¯∆ µ∆ `−m dλ = Z ∞ 0 fΛ(λ) exp −λ µp k ∆p¯∆(1 − z) dλ,
which, after an easy verification, turns out to equal φk(z) := gΛ,∆(1 − (1 − z)pk∆). The factor
φk(z) can thus be seen as the contribution of those jobs present at time 0 that arrived in the
interval [−(k + 1)∆, −k∆), for k = 0, 1, . . .. Therefore, the pgf φk(z) corresponds to that of a
mixed Poisson random variable, i.e., a Poisson random variable with random parameter
Vk :=
Λk
µ p
k
∆p¯∆= Λkpk∆r∆,
with Λk the value of the arrival rate in the interval [−(k + 1)∆, −k∆). We conclude that M is
equal in distribution to the sum of a countably infinite number of independent mixed Poisson random variables with random parameters Vk (and is consequently mixed Poisson itself).
Interestingly, the above argumentation carries over to the infinite-server model with i.i.d. general service times; say that the jobs have the same distribution as a random variable B. Write Be
for the random variable that corresponds to the residual lifetime of a job in the system in stationarity from an arbitrary time in equilibrium. Then it is readily verified that we have the following generalization of the probability in (5):
Z ∆ 0 1 ∆P(B > k∆ + s)ds = EB ∆ ξ(k, ∆) with ξ(k, ∆) := P(Be∈ [k∆, (k + 1)∆]).
From this, (5) can be recovered for the case that jobs are exponentially distributed (then Be
and B are identically distributed). The probability generating function of Mkfor general service
times is
EzMk = E exp (−Λk· EB · ξ(k, ∆) (1 − z)) ,
hence Mk is mixed Poisson with random parameter Vk= Λk· EB · ξ(k, ∆), the expected fraction
of jobs with a residual time long enough to reach zero starting from a uniform epoch in the interval [−(k + 1)∆, −k∆), for k ∈ {0, 1, . . .}. We thus obtain that M is mixed Poisson, too:
EzM = E exp − ∞ X k=0 Λk· EB · ξ(k, ∆)(1 − z) ! .
Hence M is Poisson with random parameter P∞
k=0Vk =
P∞
k=0Λk · EB · ξ(k, ∆), with i.i.d.
Λk. Considering the special case Λk = λ a.s., we see (as expected) that M is Poisson with
(deterministic) parameter λEB.
B First two moments
We now evaluate the first two moments of M . This is an interesting computation in its own right, but it also provides useful results that we can use when considering this system under a specific parameter scaling (as we do in the next subsection).
Differentiating (3) and letting z ↑ 1, we obtain a fixed-point equation,
φ0(1) = φ0(1)e−µ∆+ gΛ,∆0 (1) = φ0(1)e−µ∆+ r∆EΛ.
Hence EM = φ0(1) = r∆EΛ/(1 − e−µ∆) = EΛ/µ. This mean could have been computed more directly as well, using a standard identity for conditional means:
EM = ∞ X k=0 EMk= ∞ X k=0 E [E[Mk| Λk]] . (6)
Then observe that (Mk| Λk) is Poisson, and hence its mean equals Vk, its parameter. In the
sequel, we denote vk:= EVk. As a result, (6) equals
EM = ∞ X k=0 vk= EΛ µ ∞ X k=0 pk∆p¯∆= EΛ µ .
Note that for general service times, this generalizes to EM = EΛ · EB. For the variance we use that
φ00(1) = φ00(1)p2∆+ 2φ0(1)p∆g0Λ,∆(1) + g 00 Λ,∆(1),
and hence φ00(1) = 2φ 0(1) p ∆g0Λ,∆(1) 1 − p2∆ + g00Λ,∆(1) 1 − p2∆ = 2 EΛ µ 2 p∆ 1 + p∆ + EΛ 2 µ2 1 − p∆ 1 + p∆ .
It thus follows that, after some algebra,
Var M = φ00(1) + φ0(1) − (φ0(1))2 = EΛ µ + C Var Λ µ2 , (7) where C := (1 − p∆)/(1 + p∆).
Alternatively, one could use the ‘law of total variance’ to identify Var M :
Var M = ∞ X k=0 Var Mk= ∞ X k=0 E[Var(Mk| Λk)] + ∞ X k=0 Var(E[Mk| Λk]). (8)
Observe that, because of the ‘mixed Poisson property’, E[Var(Mk| Λk)] = EVk =: vk, and as a
result the first term at the right-hand side of (8) equals EM. The second term, which is inherently non-negative, gives rise to ‘overdispersion’, i.e., the effect that the variance of the stationary number in the system exceeds the corresponding mean. This is a distinguishing feature compared to the analogous system in which the Poissonian arrival rate is deterministic: the stationary number of jobs in an M/G/∞ queue is Poisson, and cannot accommodate any overdispersion. In order to evaluate the second term in the right-hand side of (8), we note that Var(E(Mk| Λk)) =
Var Vk, which equals
Var Vk=
1 µ2p
2k
∆ p¯2∆· Var Λ. (9)
Substituting (9) in the second term in the right-hand side of (8), we find that Var M equals (7), as desired.
From the above reasoning we conclude that formula (7) lends itself to a nice interpretation. The first term, EΛ/µ that is, is the contribution to the variance that one would have if the arrival rate would have had the deterministic value EΛ, whereas the second term needs to be added to deal with the intrinsic variability of the arrival rate.
For the model with general service times the counterpart of (7) is
Var M = EB EΛ +C (EB)¯ 2Var Λ, with C :=¯
∞
X
k=0
ξ(k, ∆)2.
B Transient behavior
As the system’s transient behavior can be dealt with in the same way, we restrict ourselves to a short account of this. We let the system start empty (for ease of presentation; a non-empty
initial condition can be analyzed without any additional difficulty). Denote by M (t) the number of jobs present at time t. Then it follows that, for n ∈ N,
M (n∆) = n X j=1 ¯ Mj,
where (for given n) ¯Mj represents the number of jobs that have arrived between (j − 1)∆ and
j∆ and are still around at time n∆, with pgf
EzM¯j = ∞ X k=0 Z ∞ 0 fΛ(λ)e−λ∆ (λ∆)k k! k X `=0 z`k ` pn−j∆ p¯∆ µ∆ !` 1 −p n−j ∆ p¯∆ µ∆ !k−` dλ = Z ∞ 0 fΛ(λ) exp −λ µp n−j ∆ p¯∆(1 − z) dλ,
as before. As the individual random variables ¯M1, ¯M2, . . . are independent, M (n∆) is mixed
Poisson: EzM (n∆)= n Y j=1 E exp −Λ µp n−j ∆ p¯∆(1 − z) .
Again the model with general service times can be dealt with as well. Following the steps that are standard by now, we obtain that M (n∆) has a mixed Poisson distribution, as
EzM¯j = Z ∞
0
fΛ(λ) exp (−λ EB ξ(j, ∆)(1 − z)) dλ = E exp (−Λ EB ξ(j, ∆)(1 − z)) .
2.2 Limit results
This section focuses on the central limit regime that results from simultaneously scaling both the sampling frequency ∆1 and the arrival rate Λ to infinity, in a controlled way. More specifically, we consider the following regime:
Λ 7→ N Λ ∆1 7→ N∆α, (10)
where the parameter N ∈ N is sent to ∞. The scaled counterpart of Λ(s) is N · Λ(N )(s), where
Λ(N )(s) :=
∞
X
j=0
Λj1[j∆N−α,(j+1)∆N−α)(s).
In words: the sampling frequency and the arrival rates are both inflated, but, importantly, this occurs at rates that are not necessarily identical. An important conclusion from the computations in the remainder of this subsection, is that the behavior of the number of jobs in the system for α > 1 (relatively high rate-renewal frequency) crucially differs from that for α < 1 (relatively low rate-renewal frequency); the case α = 1 will be treated as a special case.
We consider a sequence of systems indexed by N , where the N -th system uses a sampling fre-quency of N∆α and i.i.d. mixed Poisson arrival rates with random parameter N Λ. Let M(N ) denote the stationary number of jobs in the system. We start our exposition by a preliminary calculation, in which we compute the mean and variance of M(N ) and study their behavior for large N , and indeed observe the above-mentioned trichotomy. Then, after centering and nor-malizing M(N ), we derive a central limit theorem. This result will serve as an illustration for the reader, to create intuition as for why the scaled stationary number of customers is asymp-totically normal. However, the goal of the next section is to derive a functional central limit theorem for the (scaled) transient process M(N )(·) (for a system that is more general than the one under consideration now). Given the behavior of the transient system, asymptotic normality in stationarity could be derived by sending t → ∞.
B Qualitative behavior of first two moments: trichotomy in the variance
First, we identify the steady-state mean and variance in our scaling regime, using (6) and (7). We find that EM(N )= NEΛ µ (11) Var M(N )= NEΛ µ + N 21 − e−µ∆N −α 1 + e−µ∆N−α Var Λ µ2 ∼ N EΛ µ + N 2−α∆Var Λ 2µ , (12)
where the ‘∼’ in the last display indicates that the ratio of the left-hand side and right-hand side converges to 1 as N → ∞. We thus observe that
Var M(N ) ∼ NEΛ µ if α > 1; N2−α∆Var Λ 2µ if α < 1; N EΛ µ + ∆Var Λ 2µ if α = 1. (13)
These three cases allow an appealing interpretation, that was already mentioned in the intro-duction. For α > 1, the sampling frequency will overrule the variability of Λ and the variance of M(N ) is linear in N . Consequently, the model behaves essentially as an M/M/∞. For α < 1, we find a superlinear relation between N and Var Λ, and both the sampling frequency and the variance of Λ play a role. We indeed find that the asymptotic variance grows faster than the asymptotic mean for α < 1; in this regime the system is overdispersed. For α = 1, the variance is ‘slightly larger’ than for α > 1, but it is still linear in N . In this case the sampling frequency and the variance of Λ diverge at the same rate, resulting in a balanced variance for M(N ) that combines the two former cases.
Due to the above computation, Var M(N ) is essentially proportional to Nγ, with γ := max{1, 2 − α}. As a consequence, one may expect that, under (10),
ˇ
M(N ):= M
(N )− E M(N )
Nγ/2 (14)
converges to a (zero-mean) normally distributed random variable. It is this property that we verify now.
B Asymptotic normality
We show how to establish asymptotic normality for ˇM(N ), i.e., the centered and normalized version of the random variable M(N ). We use a direct proof, in which we evaluate the Laplace transform of ˇM(N ) and show that it converges to that of a (zero-mean) Gaussian random variable
as N → ∞; appealing to L´evy’s convergence theorem, we establish the desired convergence in distribution.
The proof is included to provide further motivation and intuition for the scaling proposed. In the next section, the setting is considerably generalized and additionally, a functional result is derived: we prove weak convergence of (a centered and normalized version of) the process M(N )(·) to a generalized Ornstein-Uhlenbeck (OU) process.
Theorem 2.1 (clt). As N → ∞, ˇM(N ) converges to a zero-mean normally distributed random variable with variance
σ2:= EΛ
µ 1{α≥1}+
∆Var Λ
2µ 1{α≤1}.
Proof. Let φ(N )(z) be the counterpart of (4) under scaling as in (10) and g(N )Λ,∆(z) the counterpart of gΛ,∆. Then φ(N )(z) = ∞ Y k=0 g(N )Λ,∆ 1 − (1 − z)e−µkN−α∆ .
We are interested in the behavior of M(N ) in the central limit regime, hence we need to analyze the limiting distribution (as N → ∞, that is) of (14). To this end, we evaluate its Laplace transform: E exp −sM (N )− EM(N ) Nγ/2 ! = exp sN1−γ/2EΛ µ φ(N ) e−s/Nγ/2 . (15)
One readily obtains that
log φ(N ) e−s/Nγ/2 = ∞ X k=0 log E exp −N Λ µ 1 − e−µN−α∆ 1 − e−s/Nγ/2 e−µkN−α∆ . (16)
The right-hand side of (16) is a sum of cumulants (which is a cumulant itself), all related to the random variable Λ. Let m` denote the `-th cumulant moment of Λ (for ` ∈ N); in particular
m1 = EΛ and m2 = Var Λ. In addition, we define
ζk(N )(s) := −N µ
1 − e−µN−α∆ 1 − e−s/Nγ/2e−µkN−α∆.
Then it follows that
log φ(N )e−s/Nγ/2= ∞ X `=1 m` `! ∞ X k=0 ζk(N )(s)`. (17)
Let us first consider the contribution of the term corresponding to ` = 1. Observe that
m1 ∞ X k=0 ζk(N )(s) = −N EΛ µ 1 − e−s/Nγ/2 ∼ −sN1−γ/2EΛ µ + s 2N1−γEΛ 2µ (18)
as N → ∞. Note that the first term in the right-hand side of (18) cancels the first term in the right-hand side of (15), so that we are left with the second term, i.e.
s2N1−γEΛ
2µ. (19)
The second term in (17), corresponding to ` = 2, gives
N2Var Λ 2µ2 1 − e−µN−α∆ 2 1 − e−2µN−α∆ 1 − e−s/Nγ/2 2 ∼ 1 2s 2N2−α−γ ∆ Var Λ 2µ . (20)
Now compare the asymptotic expansion identified in (19) and (20). In case α > 1, we have that γ = 1, so that (19) equals 12s2EΛ/µ, whereas (20) converges to zero. On the other hand, for α < 1 we have γ = 2 − α, and hence (19) converges to zero, whereas (20) behaves as
1 4s
2∆ Var Λ/µ. Finally, if α = 1, we indeed find that both terms converge to a finite positive
limit.
We now check that the terms in (17) for ` ≥ 3 vanish as N → ∞. It is readily verified that
∞ X k=0 ζk(N )(s) ` ∼ N` 1 − e−µN−α∆ ` 1 − e−`µN−α∆ 1 − e−s/Nγ/2 ` ∼ N`·µ `N−α`∆` `µN−α∆ s` Nγ`/2,
which behaves like Nδ, with δ = δ(α) := `(1 − α − γ/2) + α. In case α ≥ 1, we have γ = 1 and hence δ = `(12 − α) + α = 12` + α(1 − `) ≤ 12` + 1 − ` = 1 −2`; in case α < 1, we have γ = 2 − α and hence δ = −`α/2 + α = α(1 − 2`). Hence δ < 0 for ` ≥ 3 and the corresponding terms in (17) can indeed be neglected. We have therefore established that, as N → ∞,
log E exp −sM (N )− EM(N ) Nγ/2 ! → 1 2σ 2s2, as claimed.
3
Model with correlated arrival streams: functional central limit
theorem
In this section we generalize the central limit result of Thm. 2.1 in two ways. Firstly, we establish its functional version: the centered and normalized transient queue length process converges to a limiting process of Ornstein-Uhlenbeck type with parameters that depend on the regime. Ad-ditionally, we consider a multidimensional setting with correlated arrivals: every arrival triggers jobs in multiple queues. The resulting functional central limit theorem is a multidimensional Gaussian process. We explicitly identify the corresponding correlation structure.
Let us start by describing the mechanics of this setting. The system works as a fork-join model in which d queues are fed by a single arrival process that is constructed as the one featuring in the previous section. The service times in queue j are i.i.d. exponential random variables with mean µ−1j ; the service processes of the individual queues are independent, and also independent of the arrival process.
Clearly, both the pre-limit results and the asymptotic normality that were established in the former section can be extended to this more general setting. Moreover, we perform the same scaling as before: the sampling frequency is sped up by a factor Nα, while the (random) arrival rate is blown up by a factor N . Let Mj(N )(t) denote the number of customers at time t in the j-th queue of the scaled system.
We now present an alternative way of writing Mj(N )(t), which facilitates the use of the martingale central limit theorem. We introduce the functional
Ψ[X](t) := Z t
0
X(s) ds,
mapping the stochastic process {X(s) : s ∈ [0, t]} to a real number; observe that µjΨ[Mj(N )](t)
can be interpreted as the ‘cumulative service capacity’ in queue j in the interval [0, t]. In addition, we write ¯Λ(t) for 1tRt
0 Λ(s) ds, such that t · ¯Λ(t) is the total input rate up to time t. The scaling
counterpart of ¯Λ(t) is N ¯Λ(N )(t), where ¯ Λ(N )(t) := 1 t Z t 0 Λ(N )(s) ds.
By the law of large numbers, ¯Λ(t) and ¯Λ(N )(t) converge a.s. to EΛ (as t → ∞ and N → ∞, respectively).
With Y0(·) and Yj(·) denoting independent unit-rate Poisson processes, we can characterize
Mj(N )(t) as follows: Mj(N )(t)= Md j(N )(0) + Y0 N t · ¯Λ(N )(t) − Yj µjΨ[Mj(N )](t) . (21)
The objective in this section is to derive a d-dimensional functional central limit theorem (fclt) for the Mj(N )(·). This result characterizes the time-dependent number of customers in the scaled system and makes explicit how the correlated arrivals lead to correlation between the individual queue length processes. It will be stated and proved in the second subsection; first we consider the queues’ stationary behavior by presenting the corresponding first two moments of the system in stationarity (including the covariance of individual queues).
3.1 Qualitative behavior of first two moments
Note that the individual queue lengths are only coupled through the arrival process, so under scaling (10), the mean and variance of Mj(N ) are as in (11) and (12):
EMj(N )= N EΛ µj , Var Mj(N )= N EΛ µj + N2Var Λ µ2j 1 − pj(∆N−α) 1 + pj(∆N−α) ∼ NEΛ µj + N2−α∆Var Λ 2µj .
Hence, we find the same behavior as in (13), and interestingly, the same trichotomy is observed for the covariance between Mi(N ) and Mj(N ), for i, j ∈ {1, . . . , d}, as stated in the next lemma.
Lemma 3.1 (Covariance in M(N )). For i, j ∈ {1, . . . , d} with i 6= j, and for large N ,
Cov(Mi(N ), M (N ) j ) ∼ N EΛ µi+ µj if α > 1; N2−α∆Var(Λ) µi+ µj if α < 1; N EΛ µi+ µj +∆Var(Λ) µi+ µj if α = 1. (22)
Proof. Without loss of generality, we take i = 1 and j = 2. We first consider the non-scaled model, by studying the joint probability generating function,
Ez1M1(k∆)z M2(k∆) 2 = k Y j=1 ξjk(z1, z2)
where ξjk(z1, z2) is the contribution due to the slot between (j − 1)∆ and j∆; as z1 and z2 are
held fixed for the moment, we suppress them. Now introduce the functions
fi(r, k) := e−µi(k∆−r), gj(µ, k) :=
1 − e−µ∆
µ∆ e
where it is noted that gj(µ, k) looks like e−µ(k−j)∆ for ∆ small. In addition, we define the quantities ζjk++ := Z j∆ (j−1)∆ 1 ∆f1(r, k)f2(r, k)dr = gj(µ1+ µ2, k), ζjk+− := Z j∆ (j−1)∆ 1 ∆f2(r, k)(1 − f1(r, k))dr = gj(µ1, k) − gj(µ1+ µ2, k), ζjk−+ := Z j∆ (j−1)∆ 1 ∆(1 − f1(r, k))f2(r, k)dr = gj(µ2, k) − gj(µ1+ µ2, k), ζjk−− := Z j∆ (j−1)∆ 1 ∆(1 − f1(r, k))(1 − f2(r, k))dr = 1 − gj(µ1, k) − gj(µ2, k) + gj(µ1+ µ2, k), Using the arguments that are standard by now,
ξjk : = E ∞ X `=0 e−λ∆(λ∆) ` `! ζjk++z1z2+ ζjk+−z1+ ζjk−+z2+ ζjk−− ` ! = E expΛ∆ζjk++z1z2+ ζjk+−z1+ ζjk−+z2+ ζjk−−− 1 ∼ E exp Λ∆ 2 Y i=1 (zi− 1)e−µi(k−j)∆+ 1 − 1 !! for small ∆.
Because the contributions to the Mi(k∆) resulting from different time intervals are independent,
we obtain that Cov M1(k∆), M2(k∆) = k X j=1 ∂2 ∂z1∂z2 ξjk(z1, z2) − ∂ ∂z1 ξjk(z1, z2) ∂ ∂z2 ξjk(z1, z2) z1↑1,z2↑1 .
Now imposing scaling (10) and considering the stationary behavior by letting k tend to infinity, it is readily derived that
E zM (N ) 1 1 z M2(N ) 2 ∼ ∞ Y j=1 E exp Λ∆N1−α 2 Y i=1 (zi− 1)e−µij∆/N α + 1− 1 !! ;
observe that, for reasons of symmetry, we replaced k − j by j in the definition of the gj(µ, k).
We thus arrive at Cov M1(N ), M (N ) 2 ∼ ∞ X j=1 (EΛ2)∆2N2−2αe−(µ1+µ2)j∆/Nα+ EΛ ∆N1−αe−(µ1+µ2)j∆/Nα − ∞ X j=1 2 Y i=1 EΛ ∆N1−αe−µij∆/N α = EΛ ∆N 1−α+ Var(Λ) ∆2N2−2α 1 − e−(µ1+µ2)∆N α ,
which behaves for N large as (22).
Recall that γ = max{1, 2 − α}; due to the above computation we have that the covariance matrix of the d-dimensional vector M(N ) is essentially proportional to Nγ. Therefore, we expect that
the centered and normalized version of the joint stationary number of customers converges to a (zero-mean) d-dimensional normal random vector with covariance matrix C such that
Cij = 1{α≤1}EΛ µi + 1{α≥1}∆Var(Λ) 2µi if i = j, 1{α≤1} EΛ µi+ µj + 1{α≥1}∆Var(Λ) µi+ µj if i 6= j.
3.2 Functional central limit theorem
The objective of this section is to derive a functional limit theorem for M(N )(t), the vector describing the numbers of customers in the scaled system at time t. To this end, we consider the process ˜Mj(N )(·) := Mj(N )(·)/N , for which we have
˜ Mj(N )(t) = ˜Mj(N )(0) + N−1Y0 N t · ¯Λ(N )(t)− N−1Yj N µjΨ[ ˜M(N )](t) .
Applying the law of large numbers to Poisson processes, the following useful lemma is obtained; see e.g. [1].
Lemma 3.2. Let Y be a unit-rate Poisson process. Then for any U > 0, almost surely
lim N →∞0≤u≤Usup Y (N u) N − u = 0.
The uniform convergence in Lemma 3.2 entails that ˜Mj(N )(t) converges almost surely to the solution of the functional equation
%j(t) = %j,0+ t EΛ − µjΨ[%j](t), (23)
as N → ∞, under the proviso that ˜Mj(N )(0) converges a.s. to some value %j,0 for j = 1, . . . , d.
The solution is given by a convex mixture of the initial position %j(0) = %j,0 and the limiting
value EΛ/µj:
%j(t) = %j,0e−µjt+EΛ
µj
1 − e−µjt .
Having identified this ‘fluid limit’, the next objective is to establish an fclt for the centered and normalized process
Uj(N )(t) := Nβ ˜Mj(N )(t) − %j(t)
,
with β := 1 − γ/2 = min{12,α2}. Here we closely follow the proof in [1], where the idea is to use a martingale central limit theorem (mclt) , so as to obtain weak convergence to a (generalized) Ornstein-Uhlenbeck process. The version that we need in our setting is stated below.
Theorem 3.3 (mclt, [1]). Let {R(N )}N ∈N be a sequence of Rd-valued martingales. Assume that
the following condition on the jump sizes is met:
lim N →∞E sup s≤t R (N )(s) − R(N )(s− ) = 0; (24)
in addition, assume that, as N → ∞, h
R(N )i , R(N )j i
t
→ Cij(t)
for a deterministic function Cij(t), continuous in t for all t > 0 and for i, j = 1, . . . , d.
Then the process R(N ) weakly converges to a centered Gaussian process W with independent increments whose covariance matrix is characterized by
EWi(t) · Wj(t)T = Cij(t).
Introducing compensated unit-rate Poisson processes ˜Yj(t) := Yj(t) − t, we define
ˇ Y(N )0 (t) := Nβ−1 ˜ Y0 N t · ¯Λ(N )(t) .. . ˜ Y0 N t · ¯Λ(N )(t) , ˇ Y(N )(t) := Nβ−1 ˜ Y1 N µ1Ψ[ ˜M (N ) 1 ](t) .. . ˜ Yd N µdΨ[ ˜Md(N )](t) (25)
both define sequences of d-dimensional processes amenable for this mclt.
Lemma 3.4. Consider the d-dimensional processes ˇY(N )0 (·) and ˇY(N )(·). If α ≥ 1, then as N → ∞ these processes converge weakly to d-dimensional zero-mean Brownian motions with covariance matrices K0(t) := (t EΛ)11T and K(t) := diag{µ1Ψ[%1](t), . . . , µdΨ[%d](t)}; if α < 1
the limiting covariance matrices equal 0.
This result is in correspondence with the fact that the first martingale consists of d identical processes, whereas the other consists of d independent processes.
Proof. We start by checking the conditions of Thm. 3.3. First, observe that for each N , ˇY(N )0 (·) and ˇY(N )(·) are d-dimensional real-valued martingales; scaled and centered Poisson processes satisfy the martingale property. Also, condition (24) is met, as both for R(N ) = ˇY(N )0 and R(N ) = ˇY(N ), lim N →∞E sup s≤t R (N )− R(N )(s− ) < ∞,
Let us determine the entries of the corresponding quadratic covariation matrices: h Nβ−1Y˜0 N t · ¯Λ(N )(t)i t= N 2β−2Y 0 N t · ¯Λ(N )(t), h Nβ−1Y˜j N µjΨ[ ˜Mj(N )](t) i t= N 2β−2Y j N µjΨ[ ˜M(N )](t) , whereas, for i 6= j, h Nβ−1Y˜i N µ1Ψ[ ˜Mi(N )](t) , Nβ−1Y˜j N µjΨ[ ˜Mj(N )](t) i t= 0
(using that, for i 6= j, ˜Yi(·) and ˜Yj(·) are independent, so that their quadratic covariation is zero).
Note that β = min{12,α2}, so that 2β − 2 = min{−1, α − 2}. Now observe that for α ≥ 1 (and hence 2β − 2 = −1), by virtue of Lemma 3.2 and the fluid limit derived in (23),
lim N →∞N 2β−2Y 0 N t · ¯Λ(N )(t) = t EΛ, lim N →∞N 2β−2Y j N µjΨ[ ˜Mj(N )](t) = µjΨ[%j](t).
Thm. 3.3 yields that the processes converge weakly to d-dimensional Brownian motions with covariance matrices K0(t) and K(t), respectively.
For α < 1, we have 2β − 2 = α − 2 < −1, and hence the entries of the covariance matrices all vanish as N → ∞. Therefore, both ˇY(N )0 (·) and ˇY(N )(·) converge to a process identical to 0.
Stated below is the main theorem of this section: an fclt for the process U(N )(·). In line with earlier findings, three regimes need to be distinguished: α > 1 (the fast regime), α < 1 (the slow regime), and α = 1 (the intermediate regime).
Theorem 3.5 (fclt). As N → ∞, U(N )(·) converges weakly to a zero-mean d-dimensional Gaussian process with covariance matrix
Cij(t) := 1{α≥1} EΛ µi + %i,0e−µit (1 − e−µit) + 1 {α≤1}∆Var Λ 2µi (1 − e−2µit) if i = j, 1{α≥1} EΛ µi+ µj + 1{α≤1}∆Var Λ µi+ µj · (1 − e−(µi+µj)t) if i 6= j. (26)
Proof. We first write, for j = 1, . . . , d, using the definition of Uj(N )(t),
Uj(N )(t) = Nβ M˜j(N )(0) + Y0 N t · ¯Λ (N )(t) N − Yj N µjΨ[ ˜Mj(N )](t)) N − %j(t) ! .
Adding and subtracting %j,0, this expression is equivalent to
Nβ M˜j(N )(0) − %j,0 − Nβ %j(t) − %j,0 + Nβ Y0 N t · ¯Λ(N )(t) N − Yj N µjΨ[ ˜M (N ) j ](t) N ! ,
which, by filling out the explicit form of %j(t), simplifies to Uj(N )(t) = U1,j(N )(t) + U2,j(N )(t) + U3,j(N )(t), with U1,j(N )(t) := Uj(N )(0) − µjΨ[Uj(N )](t), U2,j(N )(t) := Nβt · ¯Λ(N )(t) − t EΛ, U3,j(N )(t) := Nβ−1 ˜Y0 N t · ¯Λ(N )(t) − ˜Yj N µjΨ[ ˜Mj(N )](t) .
We consider the three individual components separately.
◦ Component U1,j(N )(t) consists of the starting value of the process, which is assumed to converge to some value Uj(0), minus a reverting term. It is now straightforward that, as
N → ∞, U1,j(N )(·) converges to the solution of U1,j(t) = Uj(0) − µjΨ[Uj](t).
◦ Then consider U2,j(N )(·). For α ≤ 1 (and hence β = α2), the standard fclt for partial sums of i.i.d. random variables entails that, as N → ∞,
U2,j(N )(·) →√∆ Var Λ · V (·),
with V (·) a standard Brownian motion. On the other hand, the limiting process equals 0 for α > 1, as we then have β = 12 < α2.
◦ Finally, from Lemma 3.4, we conclude that U(N )3 (·) weakly converges to a d-dimensional zero-mean Brownian motion with covariance matrix K0(t) + K(t) if α ≥ 1, and to 0 else.
Using the above observations, we can now complete the proof. Each of the three regimes will be considered separately.
1. Fast regime (α > 1). We have obtained above that U(N )(·) converges weakly to the solution U (·) of the d-dimensional stochastic integral equation given by
Uj(t) = Uj(0) − µjΨ[Uj](t) + W0(t EΛ) + Wj(µjΨ[%j](t)) for j = 1, . . . , d
with W0(·), W1(·), . . . , Wd(·) independent standard Brownian motions, or, equivalently
Uj(t) = Uj(0) − µjΨ[Uj](t) + ˜Wj(t EΛ + µjΨ[%j](t))
with ˜W1(·), . . . , ˜Wd(·) standard Brownian motions (but not independent). It takes a routine
calculation to find that
Uj(t) = e−µjt Uj(0) + Z t 0 q EΛ + µj%j(s) eµjsd ˜Wj(s) .
All linear combinations of the Uj(t) are normally distributed, and we conclude that this
d-dimensional process is multivariate normal. It is readily seen that the expectations equal Uj(0)e−µjt. For the variances we find, after an elementary computation,
Var Uj(t) = e−2µjt Z t 0 (EΛ + µ j%j(s)) e2µjsds = EΛ µj + %j,0e−µjt (1 − e−µjt).
Likewise, for the covariance, with
Ui(t) := √ EΛ Z t 0 eµisdW 0(s) + Z t 0 p µi%i(s)eµisdWi(s),
it follows that, for i 6= j,
Cov(Ui(t), Uj(t)) = e−(µi+µj)tEUi(t)Uj(t) = e−(µi+µj)t EΛ · E Z t 0 e−µisdW 0(s) · Z t 0 e−µjsdW 0(s) = e−(µi+µj)t EΛ Z t 0 e(µi+µj)sds = EΛ µi+ µj (1 − e−(µi+µj)t).
This shows (26) for α > 1.
2. Slow regime (α < 1). In the slow regime, we have convergence to the solution of
Uj(t) = Uj(0) − µjΨ[Uj](t) +
√
∆ Var ΛV (t), which can be written as
dUj(t) = −µjUj(t) dt +
√
∆ Var Λ dV (t).
Therefore the Uj(·) are OU-processes, with marginal solutions, for j = 1, . . . , d:
Uj(t) = e−µjt Uj(0) + Z t 0 p ∆Var(Λ) eµjsdV (s) .
As before, we can conclude that this d-dimensional process is multivariate normal with expecta-tions Uj(0)e−µjt. Computations as above reveal that for α < 1, as claimed in (26),
Cov(Ui(t), Uj(t)) = ∆Var(Λ)
µi+ µj
(1 − e−(µi+µj)t).
3. Intermediate regime (α = 1). In this regime, we have a combination of the processes from the other cases: dUj(t) = −µjUj(t) dt + √ EΛ dW0(t) + q µj%j(t) dWj(t) + p ∆Var(Λ) dV (t).
The marginal solutions are, for j = 1, . . . , d, equal to Uj(t) = e−µjt Uj(0) + Z t 0 p EΛ + ∆Var(Λ) eµjsdW0(s) + Z t 0 q µj%j(s) eµjsdWj(s) .
Again, we conclude that this d-dimensional process is multivariate normal with expectations Uj(0)e−µjt; routine computations yield the desired covariance matrix, as given in (26). This
completes the proof.
It is interesting to study the impact of the scaling parameter α on the correlation between the individual queue lengths. Remarkably, it turns out that for α 6= 1 this correlation depends on the service rates only, whereas for α = 1 also the first and second moment of Λ play a role; see the following corollary for a result on the stationary regime.
Corollary 3.6 (Correlation coefficients). In stationarity, the correlation coefficient for i 6= j satisfies lim N →∞Corr(M (N ) i , M (N ) j ) = cij(α) · √ µiµj µi+ µj ∈ [0, 1], (27)
for some constant cij(α) ∈ [1, 2]; in addition, cij(α) = 1 for α > 1 and cij(α) = 2 for α < 1.
Proof. From Thm. 3.5, as t → ∞, Cij(t) → 1{α≥1} EΛ µi+ µj1{i6=j} + 1{α≤1}∆Var(Λ) µi+ µj .
We observe that (27) holds, with
cij(α) =
EΛ 1{α≥1}+ ∆Var(Λ) 1{α≤1}
EΛ 1{α≥1}+12∆Var(Λ) 1{α≤1}
,
which equals 1 for α > 1, and 2 for α < 1.
For α = 1, we find interesting boundary behavior: when ∆ → 0 the random arrival rate Λ becomes increasingly deterministic (i.e., Var Λ ↓ 0), cij(α) converges to 1, as was the case for
α > 1. If, on the other hand, EΛ ↓ 0, we find that cij(α) converges to 2, which equals the result
for α < 1. It is also observed that √ µiµj
µi+ µj
≤ 1 2, due to the inequality of arithmetic and geometric means.
4
Large deviations
Whereas the previous section focused on the behavior of the process M(N )(·) around the mean (by establishing an fclt), this section analyzes rare events. As it will turn out, the observed trichotomy applies here again. The results again relate to the setting with d coupled queues; for ease we first present (and prove) the results for d = 1, to return to the coupled model at the end of the section.
4.1 Univariate large deviations
As before, we study the system resulting from an arrival process with time-dependent scaled arrival rate N · Λ(N )(s), where
Λ(N )(s) = Λi1[i∆N−α,(i+1)∆N−α)(s)
with i.i.d. Λi distributed as Λ > 0. We are interested in the tail asymptotics of M(N )(t); more
specifically, our objective is to evaluate the decay rate
lim N →∞ 1 N log P MN(t) N ≥ a .
We start this section with a number of preliminary results.
◦ The following law of large numbers holds: almost surely, as N → ∞, for any t > 0,
¯ Λ(N )(t) := 1 t Z t 0 Λ(N )(s) ds → EΛ.
Hence, for all δ > 0, withSδ := [EΛ − δ, EΛ + δ],
P ¯Λ(N )t ∈Sδ → 1 as t → ∞. ◦ Furthermore, define κt Λ(N ):= Z t 0 Λ(N )(s)e−µ(t−s)ds.
Then, as observed earlier, M(N )(t) is a mixed Poisson random variable, with random parameter N · κt(Λ(N )). It is verified that (given that the queue starts empty), almost
surely, as N → ∞, for any t > 0, M(N )(t) N → %(t) := EΛ µ (1 − e −µt) =Z t 0 EΛ e−µ(t−s)ds.
The following theorem enables us to compute the large deviations for our model, using the moment generating function (mgf) of M(N )(t).
Theorem 4.1 (G¨artner-Ellis). Let (ZN)N ≥0 be a sequence of random variables such that the
mgf EeθZN exists and is finite for θ > 0 and let a > lim
N →∞EZNN. Define
I(θ) := θa − lim
N →∞ 1 N log E e θZN. Then lim N →∞ 1 N log P ZN N ≥ a = − sup θ I(θ).
Note that the maximizing value for θ is strictly positive, as we have that I(0) = 0 and I(θ) < 0 for θ < 0. Namely, the slope of limN →∞N1 log E eθZN is smaller than a, as limN →∞N1 log E eθZN
is convex in θ and its slope equals limN →∞EZNN < a in zero. Therefore, from now on we will
only consider strictly positive values for θ, which implies for instance that eθ− 1 > 0.
In order to apply G¨artner-Ellis, we first compute the mgf of M(N )(t). To this end, we use the known expression for the pgf of a mixed Poisson random variable, and then plug in z = eθ:
EeθM (N )(t) = E exp N κt(Λ(N ))(eθ− 1) = btN α ∆ c−1 Y i=0 E exp N Λi(eθ− 1) Z (i+1)N α∆ iN α∆ e−µ(t−s)ds ! × E exp N ΛbtN α ∆ c (eθ− 1) Z t btN α ∆ c ∆ N α e−µ(t−s)ds ! .
After taking the logarithm and dividing by N , we observe that for large values of N it is justified to replace t by (btN∆αc + 1) ∆ Nα, as an immediate consequence of (btN α ∆ c + 1) ∆ Nα → t. We rely on
this simplification in the rest of our analysis, i.e. we work with, as N → ∞,
1 N log btN α ∆ c Y i=0 E exp N Λ(eθ− 1) Z (i+1)N α∆ iN α∆ e−µ(t−s)ds ! , which simplifies to 1 N btN α ∆ c X i=0
log E exp N Λ(eθ− 1)e−µt e
µ(i+1)∆/Nα
− eµi∆/Nα
µ
!! .
This in turn behaves as (N large)
1 N btN α ∆ c X i=0
which, for reasons of symmetry, can be replaced by 1 N btN α ∆ c X i=0 log E exp
∆N1−αΛ(eθ− 1)e−µi∆/Nα. (28)
As before, we distinguish between the cases α > 1, α = 1, and α < 1. Informally, in the first regime, the sampling frequency is so high that the Λi’s can replaced by EΛ; the system’s
M/M/∞ behavior is also inherited when considering large deviations. In the second regime, on the contrary, rather detailed information on the distribution of Λ is required. The intermediate regime is the least insightful, as it involves a rather complicated integral.
1. Fast regime (α > 1). Our analysis heavily relies on the following proposition. Proposition 4.2. If α > 1, then lim N →∞ 1 N log Ee θM(N )(t) = %(t)(eθ− 1).
Proof. The proof consists of a lower bound (which is relatively straightforward) and an upper bound (which is more involved).
Lower bound. Due to Jensen’s inequality, (28) majorizes
1 N btN α ∆ c−1 X i=0
∆N1−αEΛ (eθ− 1)e−µi∆/N
α .
The stated lower bound follows by recognizing a Riemann sum, and by letting N go to ∞, thus obtaining EΛ Z t 0 e−µsds(eθ− 1) = EΛ µ (1 − e −µt)(eθ− 1).
Upper bound. We fix an ` ∈ (0, α). Let
KN := ∆ Nα · bN `c ∼ ∆N`−α, K¯ N := (btN α ∆ c + 1) · bN `c−1 ∼ t ∆N α−`.
Clearly, ¯KN ∼ t/KN, and bearing in mind that α > `, KN → 0 and ¯KN → ∞.
The first step is to bound (28) from above by N−1PK¯N−1
m=0 Im, where for independent Λn∼ Λ
Im :=
bN`c−1
X
n=0
log E exp∆N1−αΛ(eθ− 1)e−µ(mbN`c+n)∆/Nα
= log E exp ∆N1−α bN`c−1 X n=0 Λne−µ(mbN `c+n)∆/Nα (eθ− 1) ≤ log E exp ∆N1−α bN`c−1 X n=0 Λne−µmKN(eθ− 1) .
The main idea is to intersect the expectation with the event EN := 1 bN`c bN`c X n=1 Λn< EΛ + δ
and its complement, for an arbitrarily chosen δ > 0. To this end, we introduce JN,m :=
∆N1−α exp(−µmKN)(eθ − 1); realize that JN,m → 0 as N → ∞ (as a consequence of α > 1).
We thus have Im≤ log E exp JN,m bN`c−1 X n=0 Λn EN + E exp JN,m bN`c−1 X n=0 Λn1{Ec N} .
The idea of the remainder of the proof is to show that
lim sup N →∞ 1 N ¯ KN−1 X m=0 log E exp JN,m bN`c−1 X n=0 Λn EN = %(t)(eθ− 1) > 0, (29) whereas lim sup N →∞ 1 N ¯ KN−1 X m=0 log E exp JN,m bN`c−1 X n=0 Λn1{Ec N} , ≤ 0, (30)
which implies the claim.
Let us start by establishing (29). Because of the conditioning, the expression in the left-hand side of (29) is bounded from above by
lim sup N →∞ 1 N ¯ KN−1 X m=0 log E exp ∆N1−αbN`c EΛ + δe−µmKN(eθ− 1) = EΛ + δ KN ¯ KN−1 X m=0 e−µmKN(eθ− 1).
For ` ∈ (0, α), this Riemann sum converges to
(EΛ + δ) Z t 0 e−µsds(eθ− 1) = EΛ + δ µ (1 − e −µt)(eθ− 1) as N → ∞. Letting δ ↓ 0, we derive (29).
We continue by proving (30) by a change-of-measure argument. We define the (exponentially twisted) measure Q ≡ Qδ such that EQΛ = EΛ + δ, with mgf
EQ[e
θΛ] = E[e(θ+ˆθ)Λ]
E[eθΛˆ ] ,
where ˆθ ≡ ˆθδ > 0 maximizes θ(EΛ + δ) − log EeθΛ. The change-of-measure machinery then yields E exp JN,m bN`c−1 X n=0 Λn 1{Ec N} = EQ exp JN,m− ˆθ bN`c−1 X n=0 Λn 1{Ec N} .
Recalling that the Λn are nonnegative random variables, and realizing that JN,m− ˆθ → −ˆθ < 0
as α > 1, we conclude that for N large the expression in the previous display is majorized by Q(ENc), which converges to 12 as N → ∞ due to the central limit theorem. As a consequence, the left-hand side of (30) is bounded from above by
lim sup N →∞ 1 N ¯ KN X m=0 log1
2 = − log 2 · lim supN →∞
¯ KN
N ≤ 0.
This proves the claim.
An application of ‘G¨artner-Ellis’ then yields the following result. Proposition 4.3. If α > 1, lim N →∞ 1 N log P M(N )(t) N ≥ a ! = − sup θ
θa − %(t)(eθ− 1) = a log %(t) a
− %(t) + a.
2. Slow regime (α < 1). Here two cases need to be distinguished: Λ having either a bounded or an unbounded support.
Let us start with the unbounded case. It means that for all x > 0 there is an ε = εx > 0 such
that P(Λ > x) > ε. As a consequence, (28) is bounded from below by
1 N btN α∆ c X i=0 log E exp ∆N1−αxe−µi∆/Nα(eθ− 1)+ 1 N btN α∆ c X i=0 log ε.
The first of these two terms converges to x µ(1 − e
−µt
)(eθ− 1),
wheres the second converges to 0 (due to α < 1). Letting x → ∞ leads to
lim N →∞ 1 N log Ee θM(N )(t) = ∞.
It means that the supremum over θ ≥ 0 in ‘G¨artner-Ellis’ is attained at θ = 0, and hence we find decay rate 0. Intuitively, Λ is re-sampled relatively infrequently, and as a consequence there is a relatively high likelihood to attain a large value.
Now consider the bounded case. For ease we assume that Λ attains value in a discrete set of positive values, of which x is the largest (with probability p > 0), and y < x the one-but-largest. As above, lim N →∞ 1 N log Ee θM(N )(t) ≥ x(1 − e −µt)(eθ− 1) µ . (31) In addition, (28) is majorized by 1 N log btN α ∆ c X i=0 p exp
∆N1−αx(eθ− 1)e−µi∆/Nα+ (1 − p) exp
∆N1−αy(eθ− 1)e−µi∆/Nα,
which converges to (use x > y and α < 1) to the right-hand side of (31). Applying ‘G¨artner-Ellis’, we thus find that
lim N →∞ 1 N log P M(N )(t) N ≥ a ! = − sup θ θa − x(1 − e −µt)(eθ− 1) µ = a log x(1 − e −µt) µa + a −x(1 − e −µt) µ .
3. Intermediate regime (α = 1). In this case it is directly derived that
lim N →∞ 1 N log Ee θM(N )(t) = 1 ∆ Z t 0
log E exp [Λ∆e−µs(eθ− 1)] ds,
and hence lim N →∞ 1 N log P M(N )(t) N ≥ a ! = − sup θ I(θ) = − sup θ θa − 1 ∆ Z t 0
log E exp [Λ∆e−µs(eθ− 1)] ds
.
In the deterministic case, Λ = λ that is, this is the same result as for α > 1, but in general we cannot compute the integral. We will now consider a special case in which computation is possible.
Suppose that Λ ∼ exp(ν). In order to find the value for θ maximizing I(θ), we differentiate with respect to θ and solve
∂ ∂θ 1 ∆ Z t 0
log E exp [Λ∆e−µs(eθ− 1)] ds = Z t
0
e−µsEΛ exp [Λ∆e
−µs(eθ− 1)]
E exp [Λ∆e−µs(eθ− 1)] ds = a, (32) where, with z = ∆(eθ− 1),
E exp [Λ∆e−µs(eθ− 1)] = Z ∞ 0 ν exp [λ ∆e−µs(eθ− 1) − ν] dλ = ν ν − ze−µs
EΛ exp [Λ∆e−µs(eθ− 1)] = ν (ν − ze−µs)2
Therefore, (32) equals a = Z t 0 e−µs ν − ze−µs ds = 1 µν Z 1 e−µt 1 1 −zνxdx = log (1 − z νe −µt) − log (1 − z ν) µz = 1 µz · log ν − ze−µt ν − z
Which comes down to solving µaze−µt1−e−eµazµaz = µaν. Therefore, define f (x) := xC1
−ex
1−ex − C2 and
note that there are certain conditions on C1 and C2under which f (x) has multiple zeroes. These
could be computed numerically. Let ˆx a zero of f (x) such that ˆθ with ˆx = µa∆(eθˆ−1) maximizes I(θ). Then lim N →∞ 1 N log P M(N )(t) N ≥ a ! = −I(ˆθ).
4.2 Large deviations for the coupled model
We conclude the section by pointing out how the large devations for the coupled model (where each arrival generates work in d queues) can be dealt with. The multivariate version of the G¨artner-Ellis theorem entails that, modulo the validity of mild regularity conditions to be im-posed on the set A ⊂ Rd+,
lim N →∞ 1 N log P M(N )(t) N ∈ A ! = − inf a∈Asupθ d X i=1 θiai− lim N →∞ 1 N log E exp d X i=1 θiMi(N )(t) ! .
The problem therefore reduces to characterizing the limiting log moment generating function. It takes standard computations to verify that for α > 1, relying on specific intermediate results in the proof of Lemma 3.1,
lim N →∞ 1 N log E exp d X i=1 θiMi(N )(t) = EΛ Z t 0 d Y i=1 e−µis(eθi− 1) + 1 ds − t ! ,
whereas for α = 1 it turns out to equal
1 ∆ Z t 0 log E exp Λ∆ d Y i=1 e−µis(eθi− 1) + 1− 1 !! ds.
For α < 1 we again find decay rate 0 in the case Λ is unbounded, while in the bounded case the probability of interest is (as in the single-dimensional case) essentially determined by the largest value it can attain. Note that these results indeed reduce to the univariate case if we plug in d = 1.
5
Conclusions and outlook
In this thesis, we have introduced a new framework to model systems with overdispersed arrivals in the form of an infinite-server system with mixed Poisson arrivals in a random environment. For this model, we were able to derive stochastic-process limits and large-deviations results, and in both cases, we revealed a trichotomy in system behavior depending on the imposed scaling on system size and sampling frequency. In the process of finding these results, we encountered various possible directions for future research. Having said this, we would like finish off by exploring some extensions of the model, considering general service times and correlated queues. In order to place the model in a broader perspective, we will start this final section with a discussion on some alternative models for overdispersion.
5.1 Alternative models using mixed Poisson arrivals
The three alternative models that we will discuss have been introduced in order to deal with overdispersion. Although they all involve mixed Poisson arrival processes, it is noted that they also show crucial differences.
Jongbloed and Koole [8] considered the setting of an M/M/s system, but instead of Poisson arrivals assume mixed Poisson arrivals, so that conditional on the outcome of some random variable Λ, say Λ = λ, the number of arrivals in a unit length interval is Poisson with rate λ. In [8] they show how the additional uncertainty in Λ leads to more uncertainty in the system’s behavior and resulting dimensioning rules.
Maman [10] again considers an M/M/s-setting, yet making a stronger assumption than Jongbloed and Koole by restricting the random arrival rate to be of the form Λ = λ + λcX for some zero-mean random variable X with finite variance. Overdispersion due to parameter uncertainty thus increases with c. Maman established scaling limits for the stationary queue length behavior in the large-λ limit and a specific coupling between λ and s of the form (using unit mean exponential service times) s ≈ λ + βλc, where β is a Quality-of-Service parameter. Although the model and scalings differ from our model, Maman also found a trichotomy in scaling limits, depending on the value of c. Consequently, she identifies three levels of arrival rate variability: conventional (c ≤ 12), moderate (c ∈ (12, 1)) and extreme (c = 1). Comparing these levels with our regime, we observe striking similarities. For c ≤ 12, the situation is indeed conventional and the square-root rule applies. This corresponds to the case α > 1 in our model, where the system indeed behaves essentially like M/M/∞. Note that in our model, we considered the case where α = 1 (c = 12) separately; although there is no (persistent) overdispersion in this regime, we observed slightly different behaviour. For c ∈ (12, 1), Maman proposes the so-called qed-c staffing rule to replace the classical square-root rule; additional capacity is needed to hedge for substantial parameter uncertainty. This situation corresponds to our slow regime (0 < α < 1), where overdispersion is observed. For c = 1, Maman introduces a new regime in which qed performances can not
be achieved: the ed regime. The standard deviation of M is the same size as its mean. This situation corresponds to setting α = 0 in our model; the sampling frequency is fixed while the system size grows to ∞ and there is no scaling limit.
Mathijsen et al. [11] considered a discrete-time queue, with arrivals per slot following some mixed Poisson distribution, and a capacity of serving a deterministic number of s jobs per slot. Like in Maman [10], large-scale system limits are obtained, where A, the generic random variable describing the arrivals per slot, and s are coupled: s = EA + βσA. Again, this gives rise to a
trichotomy, depending on how the variance of A grows with the system size.
The three models in [8, 10, 11] provided a source of inspiration for this work, and share at some level the feature of the trichotomy due to the level or growth of overdispersion. The details of the models being rather different, we will not pursue a further comparison among these models here. However, this is an interesting topic for future research.
5.2 Correlated queues
Considering d queues using a single arrival process, we can think of an extension where each queue has an additional, individual arrival process (all independent). Additionally, the arrivals from the joint arrival process could be distributed over the queues according to a certain probability distribution, in two ways. Denote such arrival processes as follows:
Pi(t) = ˜Pi(t) + Pthinned by i(t),
with ˜Pi(t) for i ∈ {1, . . . , N } i.i.d. Poisson processes, and P (t) yet another independent Poisson
process that generates jobs for all queues. Using the pgf, it can be shown that the resulting processes are again Poisson. The joint pgf for the Pthinned by i(t) is
E Y i∈I zPthinned by i(t) i = e −λt· eλt·QN i=1(pi(zi−1)+1),
which equals the pgf for PN
i=1Pthinned by i(t).
Obviously, the queue length processes of the d queues are positively correlated when the jobs of P enter the queue i with probability pi, indepedently of what happens at the other sources. Defining
the Pi in such a way that the jobs of P are sent to exactly one of the i queues would provide
us with negatively correlated sources. Note that from a call center point of view, we would not expect that arrivals enter at exactly the same time (only with probability zero), whereas in a hospital setting arriving patients might generate jobs in multiple queues (e.g. doctors or machines).
5.3 General service times
We mentioned several times that results could be generalized for the case of general service times. In order to see that, we consider the stationary queue length process M with general
service times. The corresponding pgf is φ(z) = ∞ Y k=0 E exp (−ΛkEBP(Be∈ [k, k + 1])(1 − z)) = E exp (− ∞ X k=0 ΛkEBP(Be ∈ [k∆, (k + 1)∆])(1 − z)).
Imposing the usual scaling,
Λ 7→ N Λ, 1 ∆ 7→
Nα ∆ , the pgf of the scaled stationary queue length M(N ) will look like
φ(N )(z) = E exp (−
∞
X
k=0
N ΛkEB · P(Be∈ [k∆N−α, (k + 1)∆N−α))(1 − z)).
We can now determine the first two moments of M(N ):
EM(N ) = N EΛEB.
Var M(N ) = N EΛEB + N2C(N )(EB)2Var Λ.
where C(N ) =P∞
k=0P(Be ∈ [kN
−α, (k + 1)N−α))2. As N → ∞, we find that C(N ) = O(N−α)
under the assumption of exponential service times. Namely, we see that
C(N )= ∞ X k=0 e−(k+1)N−α− e−kN−α2 = (e −N−α − 1)2 1 − e−2N−α = 1 − e−N−α 1 + e−N−α = N−α+ O(N−2α) 2 + O(N−α) = 1 2N −α.
This raises the question: What is the order of C(N ) for general service times? If we find that it is of order N−α, all results in this thesis could be generalized.
Let B denote a deterministic service time. We need to impose some extra conditions to obtain that C(N )= ∞ X k=0 P(Be∈ [kN−α, (k + 1)N−α))2 = O(N−α).
