1. science
2. mathematics
3. statistics

A
How Hard are Steady-State Queueing Simulations?
ERIC CAO NI and SHANE G. HENDERSON, Cornell University
Some queueing systems require tremendously long simulation runlengths to obtain accurate estimators of
certain steady-state performance measures when the servers are heavily utilized. However, this is not uniformly the case. We analyze a number of single-station Markovian queueing models, demonstrating that
several steady-state performance measures can be accurately estimated with modest runlengths. Our analysis reinforces the meta result that if the queue is “well dimensioned,” then simulation runlengths will be
modest. Queueing systems can be well dimensioned because customers abandon if they are forced to wait in
line too long, or because the queue is operated in the “quality and efficiency driven regime” where servers
are heavily utilized but wait times are short. The results are based on computing or bounding the asymptotic
variance and bias for several standard single-station queueing models and performance measures.
Categories and Subject Descriptors: G.3 [Probability and Statistics]: Markov Processes, Queueing Theory; I.6.6 [Simulation and Modeling]: Output Analysis
General Terms: Design, Performance, Theory
Additional Key Words and Phrases: Diffusion approximations, Markovian queues, asymptotic variance
ACM Reference Format:
Eric C. Ni and Shane G. Henderson. 2013. How hard are steady-state queueing simulations? ACM Trans.
Model. Comput. Simul. V, N, Article A (January YYYY), 20 pages.
DOI:http://dx.doi.org/10.1145/0000000.0000000
1. INTRODUCTION
There is a widely held perception that using simulation to estimate steady-state performance measures for queueing systems with heavily utilized servers is hard. By
“heavily utilized” servers we mean that the fraction of time that the servers are busy
is close to 1. By “hard” we mean that the runlengths needed to obtain narrow confidence intervals with the desired coverage level are very large. On the contrary, we
will argue that for “well-dimensioned” single-station queueing systems, the simulation
runlengths needed to obtain accurate estimates are often modest. Queueing systems
can be well dimensioned because customers abandon if they are forced to wait in line
too long, or because the queue is operated in the “quality and efficiency driven regime”
where servers are heavily utilized but wait times are short. Our argument is based
on extending existing results [Whitt 2006] that support this view to additional singlestation queueing models with infinite waiting room and first-come-first served service
discipline. See Srikant and Whitt [1996] for closely related results for loss-systems,
which we do not explore.
To make this discussion more precise, let X = (X(t) : t ≥ 0) be a stochastic process
representing the number of customers or jobs in a queueing system as a function of
time, and suppose that X possesses a steady-state, i.e., there exists a random variThis work is partially supported by the National Science Foundation, under grant CMMI-1200315.
Authors’ address: School of Operations Research and Information Engineering, Cornell University
Permission to make digital or hard copies of part or all of this work for personal or classroom use is granted
without fee provided that copies are not made or distributed for profit or commercial advantage and that
copies show this notice on the first page or initial screen of a display along with the full citation. Copyrights
for components of this work owned by others than ACM must be honored. Abstracting with credit is permitted. To copy otherwise, to republish, to post on servers, to redistribute to lists, or to use any component
of this work in other works requires prior specific permission and/or a fee. Permissions may be requested
from Publications Dept., ACM, Inc., 2 Penn Plaza, Suite 701, New York, NY 10121-0701 USA, fax +1 (212)
869-0481, or permissions@acm.org.
c YYYY ACM 1049-3301/YYYY/01-ARTA $15.00
DOI:http://dx.doi.org/10.1145/0000000.0000000
ACM Transactions on Modeling and Computer Simulation, Vol. V, No. N, Article A, Publication date: January YYYY.
A:2
Ni and Henderson
able X(∞) say, for which X(t) ⇒ X(∞) as t → ∞, where ⇒ denotes convergence in
distribution. Furthermore, let f : {0, 1, 2, . . .} → R be a real-valued cost function and
suppose we wish to estimate the steady-state performance measure α = Ef (X(∞)).
For example, if f (x) = x, then our goal is to estimate the mean steady-state number of
customers in the system.
A natural estimator of α is
Z
1 t
α(t) =
f (X(s)) ds.
t 0
For a wide class of queueing systems and cost functions, it is known that as t → ∞,
√
t(α(t) − α) ⇒ σN (0, 1),
where N (0, 1) denotes a (standard) normal random variable with mean 0 and variance
1, and σ 2 is the asymptotic variance, which is also called the time-average variance.
Accordingly, an asymptotic 100(1 − κ)% confidence interval for α is α(t) ± zκ/2 σt−1/2 ,
where zκ/2 is the 1−κ/2 quantile of a standard normal random variable. The confidence
interval halfwidth is zκ/2 σt−1/2 , which is proportional to σ. Accordingly, the asymptotic
variance σ 2 , or the standard deviation σ, is an indicator of the absolute accuracy of the
estimator α(t). Similarly, relative error, which is perhaps preferable to absolute error,
is indicated through the ratio σ 2 /α2 or instead σ/α.
Whitt [1989] and Asmussen [1992] explore the magnitude of σ for a range of queueing systems and performance measures. The most important of their results shows
that for certain queueing systems with a fixed number of servers in which the servers
are utilized for a large fraction ρ < 1 of the time, when estimating the mean steadystate number of customers in the system, σ is typically of order (1 − ρ)−2 while α is of
order (1−ρ)−1 . Accordingly, when ρ is close to 1, σ/α is of order (1−ρ)−1 , and hence the
runlengths required to obtain estimators of α with small relative error are very large.
This observation has been exploited within the simulation community in “stress testing” of output-analysis algorithms. Indeed, the heavily loaded M/M/1 queue is a standard test problem for batching algorithms; see, e.g., Steiger et al. [2005].
It is now well understood that heavily loaded queueing systems as described above
require large simulation runlengths to obtain accurate estimators of steady-state performance measures, at least for steady-state moments of queue size and waiting time.
But what of other performance measures, such as the steady-state probability of delay, i.e., that a customer will have to wait for service? Perhaps more importantly, such
heavily loaded queues do not necessarily reflect “real” queueing systems. In reality,
customers will not queue forever; a common feature in queueing systems is customer
abandonment, where customers leave without receiving service if they have to wait
too long. Furthermore, it is usually the case that the number of servers in a queueing
system is chosen to ensure good customer service in the sense of short waiting times.
We call queueing systems in which customers may abandon, and/or where the number of servers is chosen to deliver short waiting times “well dimensioned” queueing
systems. (The notion of “dimensioning” queueing systems is not ours, although our use
of the term “well dimensioned” is specific to this paper; see Borst et al. 2004.) The key
question considered herein is how hard it is to accurately estimate various steady-state
performance measures associated with well-dimensioned queueing systems.
To answer this question we compute the asymptotic variance in a range of Markovian queueing models, including the M/M/∞, M/M/c, and M/M/c + M models, and
also several diffusion models. We confine our attention to these tractable models, even
though performance measures can be computed directly so that simulation is unnecessary, specifically because they are tractable. This allows us to perform the needed calACM Transactions on Modeling and Computer Simulation, Vol. V, No. N, Article A, Publication date: January YYYY.
How Hard are Steady-State Queueing Simulations?
A:3
culations. We believe that similar conclusions will hold for many less-tractable queueing systems, partly because many of our results are obtained for diffusion models that
are known to accurately approximate more general queues.
Assuming a Poisson arrival process is often appropriate, as justified by the PalmKhintchine theorem; see Karlin and Taylor [1975, p. 221], Cinlar [1972], and Nelson
[2013, p. 107]. Exponential service times are sometimes reasonable, but usually some
other distribution is more appropriate. Finally, assuming exponential customer patience times, (the “+M ” in the M/M/c + M queue) is not ideal, although the results of
Zeltyn and Mandelbaum [2005] suggest that in many queueing systems the value of
the density of the patience time distribution at 0 is the key quantity, in which case the
full distribution is unimportant and assuming an exponential distribution with rate
equal to this density value is an accurate approximation. (See the excellent surveys
Dai and He [2011; 2012] for much more on queueing systems with abandonment.) In
any case, our goal is to obtain the right “order of magnitude” of the asymptotic variance, so as long as our results are interpreted as applying to queueing systems that
are robust in this sense, we believe that confining attention to Markovian systems
is reasonable. So, for example, one should not attempt to extend our conclusions to
queueing systems with heavy-tailed interarrival and/or service time distributions, nor
to systems in which the sequences of these quantities exhibit long-range dependence.
(See Whitt 2002 for much more on such queues.)
In addition to considering the asymptotic variance of estimators, we also consider
their asymptotic bias. It turns out that either variance or bias can be more problematic in terms of delivering narrow confidence intervals that have the desired coverage,
depending on the performance measure and queueing system. In fact, bias is often the
more important property, at least in certain asymptotic regimes, as previously noted in
Srikant and Whitt [1996] for a variety of estimators of loss probabilities in loss models.
Our overall approach and philosophy is mostly adopted from Whitt [2006], who analyzed Markovian single-server queues and infinite server queues in some detail, along
with some results for multiserver queues. Indeed, on p. 411, Whitt stated that “Other
important classes of stochastic models should be analyzed in the same way.” We actually work with the same stochastic processes that Whitt did, except that we emphasize
the phenomenon of customer abandonment, we work with a greater variety of performance measures, we consider what happens in queues when the number of servers is
chosen so as to ensure that large backups do not arise, and we use slightly different
technical tools, especially for diffusion models. The paper Srikant and Whitt [1996]
mentioned above is also relevant. In that paper, asymptotic approximations are derived for the asymptotic variance and bias for four loss-probability estimators in loss
systems. Similar calculations to those we employ for diffusion models are used in Wang
and Glynn [2014], where the properties of a certain bias reduction scheme are studied.
Our primary contribution is to reinforce the meta result that for well-dimensioned
queueing systems, estimating steady-state performance measures using simulation is
not hard.
This meta-result, supported by analysis in Srikant and Whitt [1996] and Whitt
[2006] and reinforced here, shows that not only does abandonment or appropriate sizing of server pools relieve congestion (as is well understood), but the benefits extend
to simulation models in the sense that the runlengths required to obtain high-quality
confidence intervals for a number of steady-state performance measures are modest.
The remainder of this paper is structured as follows. In Section 2 we explain the
mathematical tools used to obtain the asymptotic variance and bias, and the interpretation of those quantities in choosing runlengths that deliver high-quality confidence intervals. Then, in Section 3 we review the so-called “efficiency-driven” regime,
which is the source of the common view that simulating heavily loaded queues is hard.
ACM Transactions on Modeling and Computer Simulation, Vol. V, No. N, Article A, Publication date: January YYYY.
A:4
Ni and Henderson
We then present some results in Section 4 that show that the presence of customer
abandonment changes the situation dramatically. In Section 5 we turn to the so-called
“quality and efficiency driven” regime associated with queueing systems with many
heavily-loaded servers, but where customer wait times are also modest. Finally, in
Section 6 we discuss and compare our results.
2. PRELIMINARIES
The primary queueing models we consider in this paper are the M/M/c model with
arrival rate λ, service rate µ and c servers with λ < cµ, and the M/M/c + M model
where a patience time is associated with each customer, and each customer is willing
to wait in queue only up to its patience time, at which point it abandons, i.e., leaves
without receiving service. In these systems, the sequences of customer interarrival
times, service times and patience times are mutually independent iid sequences of
exponential random variables. If X(t) gives the number of customers in the system
(for either queueing model) at time t ≥ 0, then X := {X(t) : t ≥ 0} is an irreducible,
positive-recurrent continuous-time Markov chain on the state space S = {0, 1, 2, . . .}.
Let π be the unique stationary distribution, and with an abuse of
let π(k) =
Pnotation,
∞
πk = π({k}), k ≥ 0. Let f : S 7→ R+ be a cost function and let α := k=0 f (k)π(k) be the
desired performance measure, namely the expected steady-state cost. We approximate
α by
Z t
−1
α(t) := t
f (X(s)) ds,
(1)
0
the average cost over [0, t]. The regenerative strong law of large numbers ensures that
α(t) → α as t → ∞ almost surely. See, e.g., Resnick [1992, p. 123, p. 396], and Crane
and Iglehart [1974a; 1974b; 1975] for an introduction to the regenerative method for
steady-state simulation output analysis.
Let A be the rate matrix for X, and define the function V : S → [1, ∞) by V (k) = aebk
for k = 0, 1, 2, . . ., where we leave a, b > 0 unspecified. It is straightforward to show, for
each of our queueing models, that there exist strictly positive constants a, b, β, δ such
that
AV (k) ≤ −βV (k) + δ,
(2)
for all k ∈ S, which is known as a Lyapunov drift criterion. (In this expression we take
V to be a column vector with kth component V (k), k ≥ 0, so that AV is a matrix-vector
product.) This condition implies that the chain X is “V-uniformly ergodic,” which allows
us to make a number of conclusions below; see Meyn and Tweedie [1993b, Theorem 7.1]
and also Down et al. [1995]. It turns out that one can also apply this same theory to
the diffusion models we consider in this paper to ensure that the same results apply to
those models, with only modest modifications, e.g., the rate matrix A in (2) is replaced
by the so-called generator of the diffusion process.
The Lyapunov drift criterion (2) implies [Glynn and Meyn 1996, Theorem 4.3] the
central limit theorem (CLT)
√
t(α(t) − α) ⇒ N (0, σ 2 ),
(3)
as t → ∞, where ⇒ denotes convergence in distribution, provided that for some γ > 0,
f 2 (k) ≤ γV (k) for all k. (An expression for the asymptotic variance constant σ 2 is given
below.) The functions f we consider grow at most linearly, so this condition is assured
and the CLT indeed holds. The CLT establishes that an asymptotic confidence interval
for α is given by
zσ
α(t) ± √ ,
(4)
t
ACM Transactions on Modeling and Computer Simulation, Vol. V, No. N, Article A, Publication date: January YYYY.
How Hard are Steady-State Queueing Simulations?
A:5
where z is an appropriate quantile of the standard normal distribution. (In practice we
must replace σ with an estimator thereof, but that is not important for our present purpose.) If we want the half-width of this confidence interval to be smaller than some absolute error tolerance > 0, then we require that the simulation runlength t ≥ z 2 σ 2 /2 .
Hence, the (asymptotic) variance constant σ 2 provides information on the accuracy of
the estimator α(t) in terms of the amount of simulated time that is required to obtain a
narrow confidence interval. This remains true if we want to assure that the half-width
of the confidence interval, relative to the true performance measure, is smaller than
some relative error tolerance > 0. In that case we require
t ≥ z 2 σ 2 /(2 α2 ).
(5)
(While relative error is typically the more relevant quantity, there is no additional
work to obtaining absolute error as well, so we discuss both measures.)
The estimator α(t) is almost always biased, owing to the fact that X(0) cannot usually be generated from the stationary distribution π. This bias can deteriorate the
coverage probability of the confidence interval (4). Suppose that we initiate the chain
X in some fixed state x, and let Ex and Px denote the corresponding expectation and
probability. As in, e.g., Proposition 2.1 of Awad and Glynn [2007], the bias is then
Z t
1
α(t) − α = Ex
[f (X(s)) − α] ds
t
0
Z t
1
=
[Ex f (X(s)) − α] ds
(6)
t 0
Z ∞
Z ∞
1
1
=
[Ex f (X(s)) − α] ds −
[Ex f (X(s)) − α] ds
t 0
t t
g(x)
− o(t−1 ),
(7)
=
t
where the bias constant
Z ∞
g(x) =
[Ex f (X(s)) − α] ds,
0
R∞
provided that the interchange (6) is valid, and that 0 |Ex f (X(s)) − α| ds < ∞. These
conditions are satisfied for our examples as assured by (2); see Down et al. [1995].
If we consider the (standard) regime where the desired confidence interval halfwidth
→ 0, then the required runlength according to the relative error criterion (5) is of the
order −2 . For such a runlength, the asymptotic bias is, according to (7), of the order 2 ,
which is asymptotically negligible compared to the confidence interval halfwidth.
We instead consider a different asymptotic regime, where is held fixed and the limiting behavior of σ and g(x) are considered as a function of some other quantity, such
as the arrival rate of customers and/or the number of servers, in order to understand
how desired runlengths scale with these quantities.
We adopt the philosophy that we want to ensure that confidence intervals are of a
desired width and the coverage of the confidence interval is not unduly affected by
bias. From this perspective, it is important that bias is small relative to the confidence
interval width. The confidence interval width is of the order σt−1/2 while the bias is
of the order g(x)/t. In order to achieve a narrow confidence interval, we must choose
t so that t1/2 is large relative to σ, i.e., t is large relative to σ 2 . Likewise, to ensure
that the bias g(x)/t is small, we must take t large relative to g(x). Relative to the
simulation runlength t then, the appropriate comparison is between variance σ 2 and
bias g(x). This may seem strange if one is used to measuring the quality of an estimator
ACM Transactions on Modeling and Computer Simulation, Vol. V, No. N, Article A, Publication date: January YYYY.
A:6
Ni and Henderson
through its mean-squared error, where variance and squared bias are often balanced.
The difference arises from our goal of having the bias be negligible relative to the
confidence interval width.
If we instead consider relative error, then the relative confidence interval width (relative to the performance measure α) is proportional to (σ/α)t−1/2 and the relative bias
is (g(x)/α)/t, so in terms of desired runlengths we then compare σ 2 /α2 with g(x)/α.
In the limiting regimes we consider, the confidence interval width criterion can dominate the bias criterion or vice versa, and there are also situations where neither criteria dominates. When the bias dominates the variance, or is of the same order as
the variance, then confidence interval coverage will be affected, and one might turn
to bias mitigation schemes such as initial transient deletion or careful choice of the
initial conditions.
But how can we compute the variance σ 2 and bias constant g(x) for a particular
model and choice of parameters?
It is known (see Meyn and Tweedie 1993a, Section 17.4 for the result for discretetime chains, and Steckley and Henderson 2006, Section 6 for a direct proof for the
continuous-time chains corresponding to our queueing models) that
−Ag = f˜ := f − α
(8)
where A is the rate matrix of the chain and f˜ is the “centered” cost function in the
sense that π > f˜ = π > f − α = 0. (Here > denotes the usual matrix transpose.) In fact, g
is the π-integrable solution of these equations that satisfies π > g = 0.
We can therefore compute g, and hence the bias constant g(x) for any initial condition X(0) = x, by identifying the π-integrable solution to Poisson’s equation (8) that
satisfies π > g = 0. It turns out that this also allows us to compute σ, because [Glynn
and Meyn 1996, Theorem 4.3]
∞
X
σ2 = 2
f˜(k)g(k)π(k).
(9)
k=0
Thus, in the sections to come, we will compute the stationary distribution π of the
appropriate Markov process, use this to compute α = π > f and hence f˜ = f − α where f
represents the performance measure in question, solve (8) for the π-integrable solution
g satisfying π > g = 0, and hence obtain the bias constant g(x) for any fixed initial condition (X(0) = x), and compute the variance constant using (9). The magnitude of these
quantities then tells us how “hard” it is to estimate certain steady-state performance
measures of Markovian queues using simulation. Whitt [2006] uses a very similar
approach for continuous time Markov chains, with the key differences being that we
emphasize the phenomenon of customer abandonment, we work with a greater variety of performance measures, we consider what happens in queues when the number
of servers is chosen so as to ensure that large backups do not arise, and we use a
slightly different version of Poisson’s equation. For birth-death processes the methodology above is the same as that of Whitt [1992], except that we use what Whitt calls
the “alternate form of Poisson’s equation.” For diffusions we work with the infinitesimal generator of the process, as employed in Glynn and Meyn [1996].
A similar agenda could be followed to analyse estimators other than those considered
here, provided that they can be represented as a time-average for a suitably defined
cost function f (·) as in (1).
3. THE EFFICIENCY-DRIVEN REGIME
Consider as performance measure the steady-state expected number of customers in
system (the expected occupancy), so we take f (k) = k. In this section we analyze this
ACM Transactions on Modeling and Computer Simulation, Vol. V, No. N, Article A, Publication date: January YYYY.
How Hard are Steady-State Queueing Simulations?
A:7
performance measure within what is known as the efficiency-driven regime, first looking at the M/M/c special case and then at general GI/GI/c queues. The results in
this section are known, but our derivations are included in an appendix because the
method of derivation is instructive of our general approach and is new in some cases
that we clarify there.
One might be tempted to apply a similar analysis to the steady-state delay probability, i.e., the steady-state probability that a customer will have to wait. In doing so,
one might exploit the “Poisson arrivals see time averages” property, e.g., Wolff [1989,
Section 5.16], taking f (k) = I(k ≥ c), i.e., f (k) equals 1 if k ≥ c and 0 otherwise.
Indeed, we were so tempted, but as pointed out by a referee, in the efficiency-driven
regime, this delay probability converges to 1, so there is (asymptotically) no value in
using simulation if the error precision remains fixed. Moreover, the neglected term in
the bias approximation (7) can in fact be non-negligible in the regime we consider, so
we do not attempt to analyze this performance measure in this section. More refined
tools are needed.
3.1. The M/M/c Queue
Suppose we initiate a simulation of the M/M/c queue with X(0) = 0, and consider the
efficiency-driven regime where we keep c and µ fixed while λ → cµ from below, i.e.,
ρ → 1 from below.
Rt
From (7) the bias in the estimator t−1 0 X(s) ds is asymptotically g(0)/t, which calculations in the appendix show is asymptotically −c−1 (1 − ρ)−3 t−1 (taking µ = 1). The
asymptotic variance is of the order 4c−1 (1 − ρ)−4 /t as ρ → 1. (These values agree with
the M/M/1 special case in Whitt [2006].) Recall from the discussion in Section 2 that
to obtain a desired absolute error (confidence interval halfwidth) of ±, the required
simulation runlength t is z 2 σ 2 /2 . For a 95% confidence interval, z ≈ 2, so if µ = 1,
then the desired runlength is 4σ 2 −2 ∼ 16c−1 (1 − ρ)−4 −2 as ρ → 1. To ensure the
asymptotic bias, g(0)/t, is smaller than , we require a runlength that is of the order
c−1 (1 − ρ)−3 −1 . Consequently, as ρ → 1, the variance is the dominant criterion. Considering relative error rather than absolute error, the simulation runlength needed is
asymptotically t = z 2 σ 2 /(2 α2 ) which is of the order 16c−1 (1 − ρ)−2 −2 . Also, to ensure
that the bias relative to α, g(0)/(tα), is smaller than requires a runlength of order
c−1 (1−ρ)−2 −1 , which is of the same order (in terms of ρ) as that required from the perspective of the confidence interval width. Nevertheless the constant multipliers ensure
that variance is the primary driver of runlengths. These conclusions reinforce similar
conclusions given in Whitt [2006] for the M/M/1 queue.
One way to potentially reduce bias is to choose the initial state to be “representative of steady-state conditions,” which one might interpret as meaning taking
X(0) = (1 − ρ)−1 , the approximate steady-state mean. In the appendix we compute
the exact solution to Poisson’s equation and then obtain its order as ρ → 1. This enables us to conclude that, when estimating the mean occupancy, the bias constant is
g((1 − ρ)−1 ) ∼ −(2cµ)−1 (1 − ρ)−3 , which is of the same order as g(0) so, at least in order,
bias is not reduced.
3.2. The GI/GI/c Queue
The results above shed light on what happens in heavily loaded Markovian queues.
The assumption that the arrival process is Poisson is often easily justified, owing to
the Palm-Khintchine theorem; see, e.g., Karlin and Taylor [1975, p. 221], Cinlar [1972],
and Nelson [2013, p. 107]. However, service times are often not well modeled as exponential random variables, with, e.g., the lognormal distribution often fitting empirical
data. We now review the GI/GI/c queue where the sequences of interarrival and serACM Transactions on Modeling and Computer Simulation, Vol. V, No. N, Article A, Publication date: January YYYY.
A:8
Ni and Henderson
vice times are independent and each consists of i.i.d. random variables. Such queues
defy exact analysis in general. We rely on a reflected Brownian motion approximation
for the queue-size process due to Iglehart and Whitt [1970a; 1970b]. See Whitt [2002,
Theorems 10.2.1 and 10.2.3] for a recent review. We develop similar results to those of
Whitt [1989] and Whitt [2006] using the tools sketched in Section 2. The derivations
are given in the appendix.
Let Xρ = (Xρ (s) : s ≥ 0) be the stochastic process giving the number of customers
in the system over time as a function of ρ, the utilization of the servers. Iglehart and
Whitt [1970a; 1970b] established that Xρ can be approximated by a reflected Brownian
motion (RBM) on [0, ∞) with drift −η and variance δ 2 , where η = cµ(1 − ρ) and δ 2 =
cµ((cµσU )2 + (µσV )2 ). (The exact sense in which this approximation is appropriate is
described in the appendix.) We take this approximation as exact in the sense that we
compute results (bias and variance constants) for the approximating RBM rather than
the original intractable queueing model, and use those to develop our conclusions.
Consider the steady-state mean occupancy. The bias constant when the simulation
is initiated at 0 is of the order −(1 − ρ)−3 as seen in our M/M/c results. The variance
σ 2 is of order (1 − ρ)−4 . Thus, exactly as with the M/M/c results, from the perspective
of absolute error the variance dominates, while from the perspective of relative error,
both variance and bias are of the same order, so that bias mitigation schemes should
be considered.
Accordingly, we come to the same conclusions for general GI/GI/c queues that we
did for the M/M/c queue in that the bias becomes important to consider as ρ → 1. We
might try to mitigate bias by initializing the simulation in the (deterministic) state
corresponding to the steady-state mean of the approximating RBM. In doing so, the
initial bias when estimating the mean occupancy remains of order (1 − ρ)−3 . Unfortunately, our tools are too crude to quantify the benefits from initiating a simulation of
the queue from the steady-state distribution of the diffusion (or an analog thereof in
the original queueing model), since we are confining our analysis to diffusion models
and for the diffusion the initial bias is then exactly 0.
3.3. The M/M/c + M Queue
Zeltyn and Mandelbaum [2005] defined an ED regime for queues with abandonment
in an asymptotic setting where the number of servers and the arrival rate both increase, while the patience time and service time parameters remain constant. They
assumed that c = c(λ) = (1 − γ)λ/µ where γ ∈ (0, 1) is fixed. Thus the queue has insufficient servers to meet demand. As a result, some fraction of customers must abandon to ensure stability, and this fraction approaches γ as λ → ∞. We do not analyze
this queueing system in this paper, because we believe that the quality and efficiency
driven regime that we analyze later is almost always more relevant in practice; see
Dai and He [2011; 2012] for more discussion about this regime.
4. THE IMPACT OF ABANDONMENT
In the models we considered in the previous section, customers are willing to wait indefinitely, and this leads to very large queue sizes and persistent periods of congestion
with the associated very large asymptotic variance constants. However, in almost all
true queueing systems, customers will not wait indefinitely, and this can lead to dramatic differences in performance. Consider the M/M/c + M (or Erlang-A) queue in
which customers are only willing to wait for an exponentially distributed patience time
with mean θ−1 ∈ (0, ∞). Patience times of successive customers are iid and independent of the sequences of interarrival and service times.
ACM Transactions on Modeling and Computer Simulation, Vol. V, No. N, Article A, Publication date: January YYYY.
How Hard are Steady-State Queueing Simulations?
ρ=1
4
2
0
1024
256
64
0
−2 −1
16
4
−4 −3
1 −6 −5
c
log (θ)
1
10
ρ = 1.02
10
5
0
1024
256
8
6
4
2
0
1024
256
log10(σ2)
log10(σ2)
log10(σ2)
ρ = 0.95
A:9
64
0
−2 −1
16
4
−4 −3
1 −6 −5 log (θ)
c
10
1
64
16
4
1 −3
c
−2
−1
0
1
log10(θ)
Fig. 1. Asymptotic variance σ 2 for the average number of jobs in the system under µ = 1
4.1. The M/M/∞ Queue
Suppose that θ = µ so that the mean patience time and mean service times are the
same. In this case, the queue-size stochastic process X = (X(t) : t ≥ 0) coincides with
that of the M/M/∞ queue. Even if θ 6= µ, the stochastic process X is stochastically
dominated by the queue size in an infinite-server queue with service rate min{µ, θ}.
Therefore, the M/M/∞ queue is an interesting first model to consider.
Let ρ = λ/µ. (We use this notation even though ρ no longer represents the server
utilization, which is 0.) Whitt [2006] showed that when estimating the mean steadystate number of customers in the system, the bias is −ρ/µ and the asymptotic variance
constant is 2ρ/µ. We conclude that in terms of absolute error, the asymptotic variance
and bias are both of the same order in the regime where λ → ∞ with µ held constant.
Consequently, to ensure satisfactory confidence interval coverage, bias reduction must
be explicitly considered. Interestingly, Whitt [2006] shows that when one considers
relative error in this same regime, then the bias becomes the dominant criterion. This
happens because the runlength required to achieve a given confidence interval width
relative to the true performance measure ρ is proportional to 1/ρ, while the bias relative to ρ remains constant.
4.2. The M/M/c + M Queue
In general, when 0 < θ 6= µ, the solution to Poisson’s equation can be computed but is
complicated, and we turn to numerical experimentation to illustrate the effect of abandonment. We report computational results for the asymptotic variance σ 2 and asymptotic bias under different levels of λ, c and θ, with µ = 1 held fixed, for the expected
steady-state number of customers in the system. Additional numerical results for the
performance measures steady-state probability of delay and steady-state probability
of abandonment are reported in Section 5.
Figure 1 shows that σ 2 decreases significantly in the presence of abandonment relative to the no-abandonment, ED-regime case. Inspecting the plots, we see that for
0 θ < µ, we have, approximately, that σ 2 ∝ θ−2 which is similar to the M/M/∞
case where σ 2 ∝ µ−2 , except that the abandonment rate θ replaces the service rate
µ. Recalling that σ 2 ∝ (1 − ρ)−4 in the ED regime for the M/M/c queue, this result
suggests that the reduction in asymptotic variance is of order θ2 (1 − ρ)−4 , even when
θ µ, i.e., the abandonment rate is a small fraction of the service rate. Furthermore,
the “plateau” we see in the plot of variance when ρ = 0.95 suggests that when ρ < 1, the
variance constant σ 2 is upper bounded by the M/M/c variance as θ → 0. Also, when
ρ ≥ 1, the queue without abandonment would be overloaded, but with abandonment
the results are very much like those for an M/M/∞ queue with service rate θ.
Recall that in the M/M/∞ queue, the bias constant differs from σ 2 by a multiplicative constant -2. We observe a similar scaling relationship in the M/M/c + M case in
the plots of Figure 2, which are approximately proportional to the plots in Figure 1.
ACM Transactions on Modeling and Computer Simulation, Vol. V, No. N, Article A, Publication date: January YYYY.
A:10
Ni and Henderson
3
2
1
0
1024
256
64
0
−2 −1
16
4
−4 −3
1 −6 −5
c
log (θ)
10
1
ρ = 1.02
log10(|bias|)
ρ=1
log10(|bias|)
log10(|bias|)
ρ = 0.95
5
0
1024
256
64
0
−2 −1
16
4
−4 −3
1 −6 −5
c
log (θ)
10
1
6
4
2
0
1024
256
64
16
4
1 −3
c
−2
−1
0
1
log10(θ)
Fig. 2. Absolute asymptotic bias for the average number of jobs in the system under µ = 1
5. THE QUALITY AND EFFICIENCY DRIVEN REGIME
In this section, we consider Markovian queues operating in the Halfin-Whitt regime
named in honor of Halfin and Whitt [1981], which is now also known as the “quality
and efficiency driven” regime, a name coined by Avi Mandelbaum, because not only
are customers served quickly, but the servers are also heavily utilized. This regime is
most relevant for systems with moderate to large numbers of servers, so we will be
interested in asymptotics where both the arrival rate λ and the number of servers c
increase, with the service rate µ held fixed. More precisely, we require that for some
finite constant β,
√
(1 − ρ) c → β
as √
c → ∞, where ρ = λ/(cµ). Hence, for a given value of c, the arrival rate is λ = cµ −
βµ c. When there is no abandonment (θ = 0) we must have β > 0 so that the system
is stable, but this restriction is not necessary when the abandonment rate θ > 0, since
abandonment stabilizes the system.
We continue to think of “hardness” in terms of the simulation runlength t needed
to obtain high-quality confidence intervals, although this is imperfect in the following
sense. The computational effort required to simulate to simulated-time t is proportional to the number of random variates generated over the interval [0, t], which is
proportional to λt. In the asymptotic regime considered here both c and λ increase
without bound. So the computational effort required to simulate to time t is better
represented by λt, than by t alone. In previous sections where λ was bounded, these
quantities are equivalent in order, but now that λ → ∞, they are not. Nevertheless,
we continue to estimate and report the asymptotic variance and bias constants, which
imply a desirable t, and which can in turn be scaled by λ (or cµ, since cµ and λ are of
the same order in the asymptotic regime we consider) to estimate the computational
effort required.
Exact calculations for the continuous-time Markov chain models can be performed,
but it appears to be difficult to extract insight from the results. Accordingly, we employ
a combination of analytical results from diffusion models and numerical results for
continuous-time Markov chain models.
5.1. The M/M/c Queue
Consider a sequence of M/M/c queueing systems, indexed by c = 1, 2, . . . All systems
have a fixed service rate µ and are assumed
to start out empty. The arrival rate in the
√
cth system is chosen to ensure that c(1 − ρ) is constant and equal to β > 0, where
ρ = λ/(cµ). Let Xc = (Xc (t) : t ≥ 0) be the stochastic process giving the number of
customers in the system over time in the cth system. Halfin and Whitt [1981] proved
that
Xc (·) − c
√
⇒ Y (·)
(10)
c
ACM Transactions on Modeling and Computer Simulation, Vol. V, No. N, Article A, Publication date: January YYYY.
How Hard are Steady-State Queueing Simulations?
A:11
as c → ∞, where Y is a diffusion on (−∞, ∞) with drift function
−βµ
y>0
,
µ(y) =
−µ(β + y) y ≤ 0
and infinitesimal variance 2µ, with Y (0) = 0. (In fact, Halfin and Whitt 1981 proved a
version of this result for a sequence of GI/M/c queues, but we restrict attention to a
Poisson arrival process.)
The convergence result (10) suggests the process approximation
√
Xc (·) ≈ c + cY (·).
(11)
We take this approximation as an equality, which then allows us to obtain a number of
insights that agree with our numerical results for exact calculation for M/M/c models.
In other words, we redefine Xc to be the right-hand side of (11), which is a diffusion,
and compute the order of magnitude of the variance and bias for our performance
measures for these diffusions that are indexed by c. The scaling relationship makes
this calculation quite tractable, but it is certainly not trivial, because the asymptotic
bias depends on growth rates in the solution to Poisson’s equation for the process Y (·),
which we therefore need to obtain.
To begin, consider the cost function f (x) = x, so that we wish to estimate the exRt
pected steady-state number of customers in the system, with estimator t−1 0 Xc (s) ds.
We can compute the asymptotic variance and bias of this estimator as in Section 3.2,
but to emphasize the role of the scaling we relate these quantities to similar ones associated with the process Y . Let gc be the desired solution to Poisson’s equation for Xc
and f (x) = x, and let gY be the solution to Poisson’s equation for Y and f (y) = y. Let
αc = Ef (Xc (∞))
√ be the expected steady-state cost for Xc , and define αY similarly, so
that αc = c + cαY . The functions gc and gY are related, since
Z ∞
gc (x) =
E[Xc (t) − αX |Xc (0) = x] dt
0
Z ∞
√
√
√
√
=
E[c + cY (t) − (c + cαY )|(Xc (0) − c)/ c = (x − c)/ c] dt
0
Z ∞
√
√
= c
E[Y (t) − αY |Y (0) = (x − c)/ c] dt
√
√ 0
= cgY ((x − c)/ c).
√
√
Thus, the asymptotic bias constant gc (0) = cgY (− c), and so we will need to compute
gY . Before doing so, consider the calculation of the asymptotic variance. Let σc2 be the
asymptotic variance of Xc and σY2 be the asymptotic variance of Y (for the cost function
f (x) = x). Let πc and πY be the stationary densities of Xc and Y respectively, and note
that πc (x) = c−1/2 πY (c−1/2 (x − c)). Hence,
Z ∞
2
σc = 2
(x − αX )gc (x)πc (x) dx
−∞
Z ∞
√ x−c
√
1
√ − αY
=2
c
cgY (c−1/2 (x − c)) √ πY (c−1/2 (x − c)) dx
c
c
−∞
Z ∞
= c2
(y − αY )gY (y)πY (y) dx
−∞
= cσY2 ,
so that σc2 grows linearly in c, with multiplicative constant σY2 .
ACM Transactions on Modeling and Computer Simulation, Vol. V, No. N, Article A, Publication date: January YYYY.
A:12
Ni and Henderson
√
√
So now we return to obtaining the asymptotic bias cgY (− c), for which we need to
compute gY , the πY -integrable solution to Poisson’s equation, with πY integral 0, that
satisfies the differential equation
µgY00 (y) + µ(y)gY0 (y) = −y + αY .
The solution for y ≤ 0 is
gY (y) = A1 +
y
β + αY
−
µ
µ
Z
0
y
Φ(y + β)
dy,
φ(y + β)
where the constant A1 does not depend on c and is not important for our purposes.
Now we use the fact that
yΦ(y + β)
lim
=1
y→−∞ φ(y + β)
so that
√
gY (− c) ∼ A1 −
√
√
c β + αY
− c
−
ln c ∼
µ
2µ
µ
(12)
for large c. We conclude that the asymptotic bias is of order −c/µ as c → ∞.
A cautionary note is necessary at this point. The diffusion approximation
(11) is most
√
appropriate for measuring fluctuations in the process of order c around the “central”
value c. In considering the bias starting
from initial state 0, we are considering a larger
√
fluctuation that is of order c c, so we are extrapolating past the usual range over
which we can expect the diffusion approximation to accurately match the dynamics of
the continuous-time
Markov chain it approximates. If we instead take as initial state
√
some a ≥ 0, then the diffusion
c − a c for √
√ approximation gives the asymptotic bias
constant as cgY (−a), which is of order
√ c. So we might expect that the bias starting
from initial state 0 is at least of order c, and furthermore that bias is reduced to order
√
c by choosing the initial state as c rather than 0. Our numerical experiments below
support the view that the bias starting from 0 is of order c. Furthermore, as pointed
out by a referee, the fluid model also suggests that the asymptotic bias starting from
that state is of order −c/µ; see Section 6.
Hence, the bias and variance are both of√the same order, being asymptotically linear
in c, and the bias can be reduced to order c by starting from initial state c.
Next, consider the cost function f (x) = I(x ≥ c), so that we wish to estimate
the Rsteady-state probability that an arriving customer must wait, with estimator
t
t−1 0 I(Xc (s) ≥ c) ds. Let us redefine gc to be the desired solution to Poisson’s equation
for Xc and f (x) = I(x ≥ c), and let gY be the solution to Poisson’s equation for Y and
f (y) = I(y ≥ 0). Redefine αc = P (Xc (∞) ≥ c) and αY = P (Y (∞) ≥ 0) similarly, so that
αc = αY .
Using the same arguments used for the cost function f (x) = x, we find that gc (x) =
gY (c−1/2 (x − c)), σc2 = σY2 is constant, and
Z 0
Φ(y + β)
αY
dy,
(13)
gY (y) = A2 −
µ y φ(y + β)
for y < 0. Hence the asymptotic bias constant when initiating in State 0 is
gc (0) ∼
−αY
ln c,
2µ
while the bias is reduced to constant order if we initiate in State c.
ACM Transactions on Modeling and Computer Simulation, Vol. V, No. N, Article A, Publication date: January YYYY.
How Hard are Steady-State Queueing Simulations?
A:13
(The same cautionary note above about the range of applicability of the diffusion
approximation also applies here.)
We conclude that when estimating the steady-state probability of delay, the bias is
of order ln c, while the variance is constant in c, suggesting that at least for large c,
the bias is the dominant criterion. However, given that ln c grows extremely slowly, it
is likely that both quantities are important to consider, and this remains true even if
we reduce bias by initiating in State c.
5.2. The M/M/c + M Queue
Now consider the case where θ > 0, so that customers abandon if their waiting times
are too long. Again consider a sequence of M/M/c + M queueing systems, indexed by
c = 1, 2, . . . All systems have a fixed service rate µ and are assumed
to start out empty.
√
The arrival rate in the cth system is chosen to ensure that c(1 − ρ) is constant and
equal to β, where ρ = λ/(cµ). Hence, we use exactly the same asymptotic regime as in
the previous section where customers did not abandon, except that we explicitly allow
β ≤ 0, since abandonment ensures that the systems are stable. Let Xc = (Xc (t) : t ≥ 0)
be the stochastic process giving the number of customers in the system over time in
the cth system. Garnett et al. [2002] proved that
Xc (·) − c
√
⇒ Y (·)
c
(14)
as c → ∞, where Y is a diffusion on (−∞, ∞) with drift function
−(βµ + θy) y > 0
µ(y) =
,
−µ(β + y) y ≤ 0
and infinitesimal variance 2µ, with Y (0) = 0. We see that abandonment modifies the
drift function for y > 0, but otherwise the diffusion is unchanged.
We again take the process approximation implied by (14) as exact, so that we redefine
√
Xc (·) = c + cY (·).
(15)
Consider the cost function f (x) = x, so that we wish to estimate the expected steadyRt
state number of customers in the system, with estimator t−1 0 Xc (s) ds. We can compute the asymptotic variance and bias of this estimator exactly as in the M/M/c case
above. Redefining
all the √quantities of √
interest for the case in point, we find that
√
αc = c + cαY , gc (x) = cgY ((x − c)/ c), σc2 = cσY2 , gY (y) is given, for y ≤ 0 by
(12) although with a different additive constant, and gc (0) ∼ −c/µ as c → ∞. Hence,
our conclusions for the M/M/c queue continue to hold in the case of abandonment,
although with a different variance constant σY2 . This is perhaps to be expected, since
the Halfin-Whitt regime corresponds to a situation where a nontrivial fraction of customers have to wait, but they have to wait for a vanishingly small amount of time as c
increases, and so abandonment has very little effect asymptotically.
The analysis for the cost function f (x) = I(x ≥ c) follows similar, albeit nontrivial,
lines, and we omit the details. The asymptotic bias is of order ln c while the asymptotic
variance does not depend on c.
There is an additional cost function we should consider for this model. Some managers use the steady-state probability of abandonment as a performance measure for
design, so it is worth understanding how this measure might be estimated, along with
the asymptotic bias and variance of the estimator. The discrete-time process consisting
of the indicators of whether successive customers abandon or not is not very tractable.
Fortunately, there is an alternative based on the system-size process [Garnett et al.
ACM Transactions on Modeling and Computer Simulation, Vol. V, No. N, Article A, Publication date: January YYYY.
A:14
Ni and Henderson
Probability of having to wait
Probability of abandonment
3
0.2
log10(σ2)
2
2.5
σ
log10(σ2)
Number of jobs in the system
2
1024
512
256
128
64
c 32
0
−1.5−1 −0.5
β
0.5 1
1.5 2
0
1024
512
256
128
64
c 32
0
−1.5 −1 −0.5
β
0.5 1
1.5 2
−2
−4
1024
512
256
128
c
64
1
0 0.5
32 −1.5−1 −0.5
β
1.5 2
Fig. 3. Asymptotic variance σ 2 for various performance measures for the M/M/c + M queue under the
QED regime with µ = θ = 1
2002]. When there are x customers in the system, the abandonment rate is [x − c]+ θ.
On the other hand, the long run abandonment rate is λαX , where αX is the steadystate probability that an arriving customer will abandon. Thus
θ
E[Xc (∞) − c]+ ,
λ
αX =
which can be estimated via
θ1
λt
Z
t
[Xc (s) − c]+ dt.
(16)
0
First consider the cost function f (x) = [x −√c]+ . Following our
√ now familiar argument, we redefine αX = E[Xc (∞) − c]+ = cE[Y (∞)]+ = cαY . Again, gc (x) =
√
cgY (c−1/2 (x − c)), and gY for this model is of the same form as (13) with different
constants, so the asymptotic bias is of the order c1/2 ln c and the asymptotic variance is
cσY2 . The bias can be reduced to order c1/2 by initiating the simulation in State x = c.
Of course, we are more interested in the cost function λθ [x − c]+ , and since λ ∼ cµ
as c → ∞, the asymptotic bias is of the order c−1/2 ln c and the asymptotic variance is
θ2 µ−2 σY2 /c. Hence, when estimating the probability of abandonment using (16), both
the bias and the variance decay as c grows, with the bias being asymptotically of larger
order.
5.3. Numerical Examples
We derived the results above assuming that the diffusion approximation was exact.
We now confirm the predictions of that approximation by numerically computing the
asymptotic variance σ 2 and bias for the M/M/c + M queue under the QED regime. We
present the results in Figures 3 and 4. In these plots,
√ we fix µ and θ, and for each value
of β and c considered we choose λ so that (1 − ρ) c = β. We then choose the scaling of
the c axis and vertical axis according to the predictions made by the diffusion models,
and find that both (scaled) σ 2 and bias on the vertical axis appear to be linear with
respect to (scaled) c. This suggests that the diffusion model estimates the true orders
of the variance and bias accurately.
5.4. The M/M/c + GI Queue
We conclude this section with a brief comment about M/M/c + GI queues. In these
queues, the patience times are still iid, but may not have an exponential distribution.
Zeltyn and Mandelbaum [2005] proved that (14) still holds for such queues, with the
proviso that the term θ in the drift function of the limiting diffusion Y is redefined to be
the value of the density of the patience time distribution at zero. To understand why,
note that in the QED regime, customer wait times become very small, being of order
c−1/2 as c → ∞. Hence, while a nontrivial fraction of customers have to wait, their
waiting times are almost all very small. Consequently, customers have very little time
ACM Transactions on Modeling and Computer Simulation, Vol. V, No. N, Article A, Publication date: January YYYY.
How Hard are Steady-State Queueing Simulations?
Probability of having to wait
Probability of abandonment
3
2
1024
512
256
128
64
c 32
0
−1.5 −1 −0.5
β
0.5 1
1.5 2
0.4
2
|bias|
3
|bias|
log10(|bias|)
Number of jobs in the system
A:15
1
1024
512
256
128
64
c 32
0
−1.5 −1 −0.5
β
0.5 1
1.5 2
0.2
c−1/2ln(c)
0
−1.5−1 −0.5
β
0.5 1
1.5 2
Fig. 4. Absolute asymptotic bias for various performance measures for the M/M/c + M queue under the
QED regime with µ = θ = 1
to abandon, and the patience time distribution is relevant only in terms of its behavior
near 0. Assuming the patience distribution has a positive continuous density at 0, our
conclusions about the order of the variance and bias for the performance measures we
analyzed for M/M/c + M queues remain valid for M/M/c + GI queues, assuming that
our approximation (15) does not introduce significant error.
6. DISCUSSION AND COMPARISONS
Table I summarizes our results. The values given represent the highest-order term
in the property (variance or bias accordingly) and do not include any multiplicative
constants. For example, when estimating the steady-state probability of delay in the
M/M/c queue operated in the QED regime, one can expect the asymptotic bias when
starting the simulation in State 0 to be O(ln c), while the asymptotic bias is O(1) when
starting the simulation in State c, which is more representative of steady-state conditions. These values are also proportional to the order of magnitude of the simulation
runlength t required to give a confidence interval of a fixed width in the case of variance, or to obtain a fixed bias respectively. In interpreting these results, recall that in
the QED regime, the arrival rate is of the same order as c, so that the computational
effort needed is of the order ct.
Table I. A summary of our results. Values represent the order of magnitude of the
variance or bias, ignoring multiplicative constants, for the stated steady-state performance measure and model in the stated regime. The three performance measures are
the mean number of customers in the system, the probability of delay and the probability of abandonment. The columns labelled Bias0 and Biasα respectively give the
order of the bias constant when initiating simulations with an empty system or when
initiating from an approximation to the steady-state mean α obtained from the diffusion
approximation.
Performance
Measure
EX
P (X ≥ c)
P (Ab)
Regime
Model
Variance
|Bias0 |
|Biasα |
ED
QED
QED
QED
QED
QED
M/M/c
M/M/c
M/M/c + M
M/M/c
M/M/c + M
M/M/c + M
(1 − ρ)−4
c
c
1
1
c−1
(1 − ρ)−3
c
c
ln c
ln c
c−1/2 ln c
(1 − ρ)−3
c1/2
c1/2
1
1
c−1/2
The values in Table I are appropriate when errors are measured in absolute terms. If
we instead measure errors relative to the true values of the performance measure, then
as discussed in Section 2 we must divide the variance by the square of the performance
measure, and the bias by the performance measure. Doing so yields Table II below.
The values in Table II are striking in the sense that the bias when initiating with an
empty system is of larger order than the variance in all cases, except for the M/M/c
ACM Transactions on Modeling and Computer Simulation, Vol. V, No. N, Article A, Publication date: January YYYY.
A:16
Ni and Henderson
Table II. Values are interpreted as in Table I above, except that variance is relative
to the square of the performance measure, while bias is relative to the performance
measure.
Performance
Measure
EX
P (X ≥ c)
P (Ab)
Regime
Model
Variance
|Bias0 |
|Biasα |
ED
QED
QED
QED
QED
QED
M/M/c
M/M/c
M/M/c + M
M/M/c
M/M/c + M
M/M/c + M
(1 − ρ)−2
c−1
c−1
1
1
1
(1 − ρ)−2
1
1
ln c
ln c
ln c
(1 − ρ)−2
c−1/2
c−1/2
1
1
1
queue in the ED regime, where the two properties are equal in magnitude. This suggests that bias should receive careful consideration in simulations of heavily-loaded
queues, in agreement with results for loss models in Srikant and Whitt [1996], and
results for the M/M/∞ queue in Whitt [2006]. To mitigate this bias, Whitt [2006]
suggested simulating starting from an initial state where all servers are busy, with
residual service times sampled from the equilibrium residual-life distribution, instead
of starting with an empty system. Our results suggest that this would substantially
reduce bias in the QED regime, as seen in the final columns of the tables above. For
example, in estimating the expected steady-state number of customers√in the M/M/c
queue in the QED regime, the absolute bias would then be of the order c rather than
c, and in estimating the probability of delay the bias would be of order 1 rather than
ln c. However, as seen in the ED results for EX, the order of the bias reduction may
depend on the performance measure; an order of magnitude in bias reduction is not
guaranteed.
Even more substantial bias reduction might result if the initial state of the simulation is randomly chosen with a distribution that is related to the stationary distribution of the heavy-traffic approximation. While we expect bias to be reduced, our
methods cannot shed light on the effect, because we analyze the bias reduction from
the perspective of the heavy-traffic approximation itself. Thus our prediction of the
resulting bias would be 0, and a deeper analysis is needed.
It is interesting that in Table II the asymptotic variances relative to the square of
the mean in estimating EX in the QED regime are of order c−1 , showing that the
simulation runlength needed to obtain confidence interval widths with given relative
error shrinks as c → ∞. It is worth keeping in mind that in the QED regime the arrival
rate is approximately proportional to c, so that the total number of customer arrivals
simulated is constant. This phenomenon was noted in Srikant and Whitt [1996] and
in Whitt [2006] for related performance measures and systems. This is a striking observation, especially when one compares it with the situation in the ED regime in
the absence of abandonment, where the number of customer arrivals that need to be
simulated is of the order (1 − ρ)−2 , which grows extremely rapidly as ρ → 1.
Although the relative bias is of equal or larger order than the relative variance in all
cases, it is important to keep in mind the discussion from Section 2 that in the usual
asymptotic setting where the desired confidence interval width → 0, the confidence
interval width will eventually dominate the bias. The comments above apply to the
setting where is fixed and ρ → 1 (in the case of ED) or c → ∞ (in the case of QED).
A referee suggested that fluid models underlie and explain the large difference in
bias results for the multi-server (c remaining bounded) and many-server (c → ∞)
regimes that we obtained through tractable diffusion models. This suggests that our
results, and others, might instead be obtained by studying the even-more tractable
fluid models associated with these processes. For example, for the M/M/c queue in the
QED regime of Section 5, the fluid model initiated in State 0 is x0 (t) = λ − µ min(c, x)
ACM Transactions on Modeling and Computer Simulation, Vol. V, No. N, Article A, Publication date: January YYYY.
How Hard are Steady-State Queueing Simulations?
A:17
with x(0) = 0. The solution when λ < cµ is
x(t) =
λ
(1 − e−µt ) t ≥ 0.
µ
The corresponding approximation for g(0) is
Z ∞
λ
λ
g(0) ≈
(x(t) − ) dt = − 2 .
µ
µ
0
This is asymptotically of order −c/µ since λ ∼ cµ in the QED regime, matching our
order of the bias computed using the diffusion model.
A. APPENDIX
Here we provide further details on the calculations in Section 3.
A.1. The M/M/c Queue
The stationary distribution π is given by
" c−1
#−1
X (cρ)k
(cρ)c 1
π0 =
+
, and
k!
c! 1 − ρ
k=0
(
k
π0 (cρ)
, 0 < k < c,
k!
πk =
ρk cc
π0 c! , k ≥ c.
(17)
(18)
With this stationary distribution and the cost function f (i) = i, the long-run average
cost π > f is cρ + Cρ/(1 − ρ), where the constant C is the delay probability
C=
∞
X
(cρ)c /[(1 − ρ)c!]
πk = Pc−1 (cρ)k
.
(cρ)c
k=c
k=0 k! + (1−ρ)c!
(19)
Poisson’s equation (8) can be solved directly to yield
g(0) = K1 − K2 ,
(
Pj Pi−1
K1 − K2 + µ1 i=1 k=0
g(j) =
aj 2 + bj − K2 ,
α−k (i−1)!
,
(cρ)i−k k!
0 < j < c,
j ≥ c,
where
1
,
2cµ(1 − ρ)
ρ
b=a 1+2
(1 − C) − cρ ,
1−ρ
a=
and the constant K1 is determined by
c
i−1
1 X X α − k (i − 1)!
K1 = ac + bc −
µ i=1
(cρ)i−k k!
2
(20)
k=0
We then select K2 so that π > g = 0, which gives
K2 = (1 − C)K1 +
j i−1
c−1
∞
1 X X X α − k (i − 1)!
cc X 2
πj
+
π
(aj + bj)ρj .
0
µ j=1 i=1
(cρ)i−k k!
c! j=c
k=0
ACM Transactions on Modeling and Computer Simulation, Vol. V, No. N, Article A, Publication date: January YYYY.
(21)
A:18
Ni and Henderson
The asymptotic bias when initiating with an empty system is g(0) = K1 − K2 . Thus
we need to understand the asymptotics of this quantity as ρ → 1, while we hold c and
µ fixed. First consider K1 as in (20), which, in turn, depends on C, a and b, all of which
depend on ρ. First a is of order (1 − ρ)−1 . Second, 1 − C can be seen to be of the order
(1 − ρ) as ρ → 1, and thus b is asymptotically (1 − ρ)−1 (3 − 2c)/(2cµ). Thus, the first two
terms in (20) are of order (1 − ρ)−1 . As to the last term, as ρ → 1,
α = α(ρ) ∼
∞
X
jπj = π0
j=c
∞
cc X j
jρ .
c! j=c
P∞
One can verify that j=c j r ρj ∼ (1 − ρ)−(r+1) r! for r = 1, 2, . . ., and that π0 cc /c! ∼ 1 − ρ.
Thus, α ∼ (1 − ρ)−1 as ρ → 1. Hence α − k in (20) is of order (1 − ρ)−1 , which, when
taken out as a common factor, leaves a quantity that is bounded in ρ as ρ → 1. We
conclude that K1 is at most of order (1 − ρ)−1 as ρ → 1.
Using similar reasoning, we see that in (21), the first two terms are asymptotically
of order 1 and the final term is asymptotically (1 − ρ)−3 /(cµ). Hence g(0)/t is asymptotically −(1 − ρ)−3 /(cµt) as ρ → 1, agreeing with the M/M/1 special case discussed in
Whitt [2006].
Turning to the asymptotic variance, σ 2 , substituting the expressions for π and g into
(9) gives
j X
c−1
i−1
X
2X
α − k (i − 1)!
σ = 2K1
πj (j − α) +
πj (j − α)
i−k
µ
(cρ)
k!
j=0
j=0
i=1
c−1
X
2
k=0
+ π0
∞
c X
c
c!
(j − α)(aj 2 + bj)ρj .
j=c
Using the same asymptotic-order calculations, the dominant term in this expression is the last one, and it is asymptotically 4(cµ)−1 (1 − ρ)−4 , again agreeing with the
M/M/1 special case in Whitt [2006].
This particular form of derivation of the asymptotic constants where we directly
compute the solution to Poisson’s equation and then estimate the asymptotic order of
the expressions is, to the best of our knowledge, new for the M/M/c queue, although
the order of the constants has been known for some time.
A.2. The GI/GI/c Queue
Consider a family of queueing systems all of which have c servers serving jobs in firstin-first-out order, indexed by ρ ∈ (0, 1), constructed as follows. Let U = (Ui : i ≥ 1)
denote an iid sequence of unscaled interarrival times, and let V = (Vi : i ≥ 1) denote
an iid sequence of service times. We assume that EV1 = µ−1 , EU1 = (cµ)−1 , and that
2
both U1 and V1 have finite variances σU
and σV2 respectively. In the ρth system, the
service time sequence is V , and the interarrival time sequence is ρ−1 U , so that the ith
interarrival time is Ui /ρ. All systems are initially empty at time 0, so that the first
customer arrives at time U1 /ρ. Let Xρ (s) be the number of customers in the system at
time s in the ρth system and let Xρ = (Xρ (s) : s ≥ 0) be the corresponding stochastic
process.
For constants a, b > 0, let aXρ (·/b) be the stochastic process taking the value aXρ (t/b)
at time t. Iglehart and Whitt [1970a; 1970b] proved that
·
⇒ R(·; −cµ, cµ((cµσU )2 + (µσV )2 ), 0),
(22)
(1 − ρ)Xρ
(1 − ρ)2
ACM Transactions on Modeling and Computer Simulation, Vol. V, No. N, Article A, Publication date: January YYYY.
How Hard are Steady-State Queueing Simulations?
A:19
where ⇒ denotes convergence in distribution of stochastic processes as in Billingsley
[1968], and R(·; r0 , r1 , r2 ) is a reflected Brownian motion (RBM) on [0, ∞) with drift r0 ,
variance r1 , and initial state r2 . Using scaling properties of RBM as in Whitt [2006],
(22) suggests the process approximation
Xρ (·) ≈ R(·; −η, δ 2 , 0),
2
2
(23)
2
where η = cµ(1 − ρ) and δ = cµ((cµσU ) + (µσV ) ). The stationary distribution of this
RBM is known [Harrison 1990, p. 94] to be exponential with mean δ 2 /(2η).
Consider the steady-state mean occupancy. We approximate all quantities for this
performance measure (bias, variance etc) by the corresponding values for the approximating RBM (23). Accordingly,
α=
δ2
,
2η
which simplifies to (1 − ρ)−1 for M/M/c queues where δ 2 = 2cµ, agreeing (in order
as ρ → 1) with the exact result. Poisson’s equation for the RBM takes the form of a
differential equation (see Mandl 1968, p. 39 and Karlin and Taylor 1981, p. 305), and
is
δ 2 00
g (x) − ηg 0 (x) = −x + α
g 0 (0) = 0.
2
The solution we seek (with zero steady-state mean) is
g(x) =
δ4
x2
− 3.
2η 4η
Accordingly, the bias constant when the simulation is initiated at 0 is g(0) or −δ 4 /(4η 3 )
which is of order −(1 − ρ)−3 . In the special M/M/c case where δ 2 = 2cµ, the bias
constant is −[cµ(1 − ρ)3 ]−1 . In either case, the order of the bias constant is (1 − ρ)−3 as
reflected in our earlier results.
The variance is
Z ∞
σ2 = 2
f˜(x)g(x)π(dx) = δ 6 /(2η 4 ),
0
which is of order (1 − ρ)−4 in general as pointed out in Whitt [1989], and equal to
4(cµ(1 − ρ)4 )−1 in the M/M/c case where δ 2 = 2cµ.
If we initialize the RBM in the state α = δ 2 /(2η) instead of 0, then the initial bias
when estimating the mean occupancy is −δ 4 /(8η 3 ) which remains of order (1 − ρ)−3 .
ACKNOWLEDGMENTS
It is a privilege to contribute to this issue honoring Don Iglehart and his distinguished career. As a graduate
student, I (Henderson) took classes, including an independent reading course, from Don. I could not have
asked for a better role model. Don is simultaneously a scholar of the highest quality, a superb mentor, and
one of the kindest people you could hope to meet. I will strive to emulate his humble excellence to the best
of my ability.
We are grateful to the editorial team for highly insightful comments that greatly improved the paper.
REFERENCES
S. Asmussen. 1992. Queueing simulation in heavy traffic. Mathematics of Operations Research 17 (1992),
84–111.
H. P. Awad and P. W. Glynn. 2007. On the theoretical comparison of low-bias steady-state estimators. ACM
Transactions on Modeling and Computer Simulation 17, 1 (2007), Article 4.
P. Billingsley. 1968. Convergence of Probability Measures. Wiley, New York.
ACM Transactions on Modeling and Computer Simulation, Vol. V, No. N, Article A, Publication date: January YYYY.
A:20
Ni and Henderson
S. Borst, A. Mandelbaum, and M. I. Reiman. 2004. Dimensioning large call centers. Operations Research 52
(2004), 17–34.
E. Cinlar. 1972. Superposition of point processes. In Stochastic Point Processes: Statistical Analysis, Theory,
and Applications, P. A. W. Lewis (Ed.). Wiley Interscience, New York, 549–606.
M. A. Crane and D. L. Iglehart. 1974a. Simulating Stable Stochastic Systems, I : General Multiserver
Queues. J. ACM 21, 1 (1974), 103–113.
M. A. Crane and D. L. Iglehart. 1974b. Simulating Stable Stochastic Systems, II: Markov Chains. J. ACM
21, 1 (1974), 114–123.
M. A. Crane and D. L. Iglehart. 1975. Simulating Stable Stochastic Systems: III. Regenerative Processes
and Discrete-Event Simulations. Operations Research 23, 1 (1975), 33–45.
J. G. Dai and S. He. 2011. Queues in service systems: customer abandonment and diffusion approximations. In Tutorials in Operations Research: Transforming Research into Action, Joseph Geunes (Ed.).
INFORMS, Hanover MD, Chapter 3, 31–59.
J. G. Dai and S. He. 2012. Many-server queues with customer abandonment: a survey of diffusion and fluid
approximations. Journal of Systems Science and Systems Engineering 21 (2012), 1–36.
D. Down, S. P. Meyn, and R. L. Tweedie. 1995. Exponential and uniform ergodicity of Markov processes.
Annals of Probability 23, 4 (1995), 1671–1691.
O. Garnett, A. Mandelbaum, and M. Reiman. 2002. Designing a call center with impatient customers. Manufacturing & Service Operations Management 4, 3 (2002), 208–227.
P. W. Glynn and S. P. Meyn. 1996. A Liapounov bound for solutions of the Poisson equation. Annals of
Probability 24 (1996), 916–931.
S. Halfin and W. Whitt. 1981. Heavy-traffic limits for queues with many exponential servers. Operations
Research 29, 3 (May - Jun. 1981), 567–588. http://www.jstor.org/stable/170115
J. M. Harrison. 1990. Brownian Motion and Stochastic Flow Systems (2nd ed.). Krieger, Malabar Florida.
D. L. Iglehart and W. Whitt. 1970a. Multichannel queues in heavy traffic I. Advances in Applied Probability
2 (1970), 150–177.
D. L. Iglehart and W. Whitt. 1970b. Multichannel queues in heavy traffic II: sequences, networks, and
batches. Advances in Applied Probability 2 (1970), 355–369.
S. Karlin and H. M. Taylor. 1975. A First Course in Stochastic Processes (2nd ed.). Academic Press, Boston.
S. Karlin and H. M. Taylor. 1981. A Second Course in Stochastic Processes. Academic Press, Boston.
P. Mandl. 1968. Analytical Treatment of One-dimensional Markov Processes. Springer-Verlag, New York.
S. P. Meyn and R. L. Tweedie. 1993a. Markov Chains and Stochastic Stability. Springer-Verlag, London.
S. P. Meyn and R. L. Tweedie. 1993b. Stability of Markovian processes III: Foster-Lyapunov criteria for
continuous-time processes. Advances in Applied Probability 25 (1993), 518–548.
B. L. Nelson. 2013. Foundations and Methods of Stochastic Simulation. International Series in Operations
Research & Management Science, Vol. 187. Springer, New York.
¨
S. I. Resnick. 1992. Adventures in Stochastic Processes. Birkhauser,
Boston.
R. Srikant and W. Whitt. 1996. Simulation run lengths to estimate blocking probabilities. ACM Transactions
on Modeling and Computer Simulation 6, 1 (1996), 7–52.
S. G. Steckley and S. G. Henderson. 2006. The error in steady-state approximations for the time-dependent
waiting time distribution. Stochastic Models 23 (2006), 307–332.
N. M. Steiger, E. K. Lada, J. R. Wilson, J. A. Joines, C. Alexopoulos, and D. Goldsman. 2005. ASAP3: A batch
means procedure for steady-state simulation analysis. ACM Transactions on Modeling and Computer
Simulation 15, 1 (2005), 39–73.
R. J. Wang and P. W. Glynn. 2014. On the Marginal Standard Error Rule and the testing of initial transient
deletion methods. (2014). Submitted for publication.
W. Whitt. 1989. Planning queueing simulations. Management Science 35 (1989), 1341–1366.
W. Whitt. 1992. Asymptotic formulas for Markov processes with applications to simulation. Operations Research 40, 2 (1992), 279–291.
W. Whitt. 2002. Stochastic-Process Limits. Springer, New York.
W. Whitt. 2006. Analysis for design. In Handbook of Simulation, S. G. Henderson and B. L. Nelson (Eds.).
Elsevier, Amsterdam, 381–413.
R. W. Wolff. 1989. Stochastic Modeling and the Theory of Queues. Prentice Hall, Englewood Cliffs NJ.
S. Zeltyn and A. Mandelbaum. 2005. Call centers with impatient customers: many-server asymptotics of the M/M/n + G queue. Queueing Syst. Theory Appl. 51, 3-4 (2005), 361–402.
DOI:http://dx.doi.org/10.1007/s11134-005-3699-8
ACM Transactions on Modeling and Computer Simulation, Vol. V, No. N, Article A, Publication date: January YYYY.
How Hard are Steady-State Queueing Simulations?
Received July 2013; revised July 2013; accepted July 2013
ACM Transactions on Modeling and Computer Simulation, Vol. V, No. N, Article A, Publication date: January YYYY.
A:21