Ch. 11. Numerical methods in queueing theory

D. N. Shanbhag and C. R. Rao, eds., Handbook of Statistics, Vol. 21 © 2003 Published by Elsevier Science B.V.

I | J_ 2[_

Numerical Methods in Queueing Theory

D. Heyman 1. Introduction Queueing theory relies heavily on Markov chains and Laplace transforms to obtain analytical results. Markov chains frequently describe the dependency in queueing processes. Laplace transforms are often convenient for solving differential equations and describing sums of i.i.d, random variables that occur in solving queueing models. Probability generating functions are used for similar reasons. Consequently, this entry describes numerical methods for analyzing Markov chains and for inverting Laplace transforms and probability generating functions. A numerical method that is used extensively is not considered - simulation. Random number generation, experimental design and output analysis raise substantial statistical issues. The numerical methods that are used to address these issues are different from the methods described here. The numerical solution of a queueing model consists of more than obtaining a tractable algorithm. Computer roundoff errors can cause mathematically correct algorithms to produce inaccurate numerical solutions, so numerical accuracy of an algorithm has to be understood. 1.1. Examples The numerical methods we shall consider appear when analyzing the M / G / 1 queue. This queueing model has Poisson arrivals at rate )~. There is a single server that takes a random time with density function g to serve a customer. Customers wait as long as necessary to receive service and are served in order of arrival. Let X (t) be the number of customers present (in queue and in service) at time t, Xn be the number present just after the nth service completion, and W. be the waiting time of the nth customer served. D. Kendall showed that Xn, n = 1, 2 . . . . . is a Markov chain. The transition matrix is

P =

ao ao 0 0 0

al al ao 0 0

a2 a2 al ao 0

a3 a3 a2 al ao

...~ ...| ""l ...| ' ...|

/ 407

(1)

D. Heyman

408 where ak = f0o° e -xt ~ g ( t )

~'~z j ~ j=l

dt. A closed form equation for

P[X(t)

=jJXo =i]e-Stdt

(2)

f~

is available (Cohen, 1969, Eq. 4.49). The transform inversion techniques described in Section 2 allow numerical inversion of this double transform. Let/z -~ be the mean service time, p -- X//z and Y~ij(n) = P [ X n = j l X o = i]. When p < 1 the following limits exist and are independent of i. zrj---- lira /'/--~ OO

rcij(n)---=n lim P[Xn----j] ---> (X3

(3)

and

W(t)= lira e[wn~
t~>0.

(4)

/ 2 - + OO

Moreover, Jr0 = W(0) = 1 - p and W has density for t > 0 (say w) because the service times do. Let ~(s) = f o e-Stw(t) dt; then ~ ( s ) = (l - p )

sg(s) s - x + X~(s)'

~ [ s ] / > 0,

(5)

where ~ is the Laplace transform of g and 91[s] is the real part of s. Let r?(z) = y~-zJrc; with Izl ~< 1" then J=U J ;~(z) = (1 - p)

(z - 1)~()~ - ;~z) z - Z + X~(Z - X z ) '

Izl ~< 1.

(6)

Another important operating characteristic is the busy period, which starts when a customer arrives at an empty system and ends the next time the system is empty. Let b be the density function of a busy period and/~ be its Laplace transform. Then/~ is the smallest positive solution of (7) There are known formal solutions of (5) and (7); these are infinite series and have not been found suitable for numerical computation. Numerical inversion of the Laplace transforms is preferable. The {rrj } can be computed from the balance equations of the Markov chain governed by P, but numerical inversion of the probability generating function in (6) is typically more convenient because it avoids calculating the {ak}. When the customers are partitioned into two classes, one of which is given priority over the other, numerical transform inversion is practical and no potentially computable time-domain expressions are known because the Laplace transform of the "effective service-times" has /~ as part of the argument of ~. The equations were derived by

409

Numerical methods in queueing theory

D. R Gaver and may be found in (Heyman and Sobel, 1982, Section ll-5). Similar issues occur when the customers are served "shortest service-time first" (Tak~ics, 1964). Finite Markov chains are used to model queues with limited (or no) waiting space. When there are K waiting positions, the embedded Markov chain with M / G / 1 assumptions is the (K + 1) x (K + 1) Northwest corner of the matrix P in (1) with the last column augmented so that the row sums are one. The steady-state probabilities can be obtained by numerically solving the balance equations. The special structure of this transition matrix (upper-Hessenberg or skip-free-to-the-right) permits a solution by recursion. The following example is more typical of those models that are solved by the numerical method described in Section 4.2. There are two queues, each with s servers and b waiting positions. A customer that tries to enter queue 1 when all of its waiting positions are full will attempt to join the second queue. Customers that try to join the second queue when all of its waiting positions are full are lost. Assume that the arrival processes are Poisson and the service times at all servers are exponential. Let )vi be the arrival rate (including lost customers but not including overflow customers) and/*i be the service rate at queue i. Let Ni (t) be the number of customers in queue i at time t. The bivariate process {Nl(t), N2(t), t ~> 0} is Markov. Let Q be the generator; Q has the block-partitioned form Bo A1 Q =

Co B1 A2

C1 B2

C2 ,

Ab-i

Bb-1 Ab

(8)

Cb-1 Bb

The blocks are Ak =/*1 min(k, s)I, Ck=)vll,

k = 1, 2 . . . . . b,

k=0,1,2 ..... b-l,

(9) (10)

where I is the (b + 1) x (b + 1) identity matrix. For k = 0, 1 . . . . . b - 1 and i = 1,2 . . . . . b, Bk(i,i -- 1) = b~2 rain(i, s),

Bk(i -- 1, i) =)~2

(11)

Bb(i -- 1,i) = ~ I -t- 5L2.

(12)

and is zero otherwise. Finally, Bb(i,i - 1) = #2min(i, s),

The stationary distribution of Q can be used to calculate performance measures such as the probability that a customer is lost and the probability that a customer arriving at queue 1 is served in queue 2.

410

D. H e y m a n

2. Numerical inversion of Laplace transforms The numerical inversion of Laplace transforms has long been thought to be difficult, but that is not correct for probability distributions. Several programs for inverting Laplace transforms were available at Bell Laboratories in the late 1960's. We describe algorithms for inverting Laplace transforms and probability generating functions (z-transforms) based on Abate and Whitt (1992a, 1992b); the former contains an extensive history and literature review. The expositions by Abate and Whitt contain many probabilistic insights that are lost in the abbreviated treatment given here. Let f be a non-negative function on the positive real line, and let f be the Laplace transform of f :

f(s) =

f0 °

e-St f ( t ) dt,

(13)

where s is a complex variable. The standard inversion formula is the Bromwich contour integral

1 fa+i~

f(t) = ~ i Ja-ioo

e'tf(s) ds,

i = xfZT,

(14)

where a is chosen so that f(s) has no singularities on or to the right of the vertical line s = a in the complex plane. (a = 0 suffices when f is a bounded continuous probability density.) Writing s = a + iu and using the identity from complex variable theory

(15)

e s = e a [cos(u) + i sin(u)]

the right side of (14) yields the following string of equations (,~ [z] denotes the imaginary part of z):

f(t) = ~

lF

at f = e2-7

[cos(ut) + i sin(ut)]f(a + iu) du oo

zc

{~[f(a+iu)]cos(ut)- .~[f(a+iu)sin(ut)]}du

oo

2eat Jr fo ~ 9~[f(a + iu)] cos(ut) du.

(16)

The integral in (16) is evaluated symbolically with the trapezoidal rule. There is a known error bound for the discretization error introduced. Let fh be the value of the integral when a step size of h is used. Then,

f ( t ) ~ fh(t) =

heat Yg

~o

N(f)(a) + 2heat ~ 9 l ( f ) ( a Yg

k=l

+ ikh) cos(kht).

(17)

Numerical methods in queueing theory

411

Taking h = rr/(2t) and a = A/(2t) yields e a/2 °c Z(--1)km(f)

eA/2 ( A )

fh(t) = fit(f) \ 2 t J q- ~

(A+2kzri) 2t

"

(18)

k=l

The Poisson summation formula (see, e.g., Franklin, 1964, p. 514) can be used to show that the discretization error f(t) - fh (t) = ea (say) is given by

ed = ~ _ e - k A f[(2k + 1)t].

(19)

k=l

When f is either a distribution function or a complementary distribution function, f(t) ~< 1 for all t, so e-A

ed <~ 1 -- e -A '

(20)

which is approximately e A when e - A is small. Thus, when A = V l o g 1 0 the discretization error is no larger than 10-×; A = 19.1 achieves a 10 -7 discretization error. This is the only place a probabilistic property of f is used. The first n terms in (18) are summed explicitly and the next m terms are summed via Euler summation (Franklin, 1964). This introduces a truncation error that can be estimated and reduced by increasing n. Abate and Whitt recommend taking m = 11 and n = 15; m and n are increased when the error estimate is larger than desired. Let E(m, n, t) be the Euler sum as described above. Then m

m

E(m,n,t) = ~"~m k 2- Sn+k(t),

(21)

k=0

where

e A/2 ~ /A\ eA/2 n sn(t)--- 2t fit(f)~27) + T ~-~ak(t)

(22)

k=l

and k

~

ak(t)=(--l) ~t~(f)(

A +

2kTci ~ ).

(23)

Eqs. (18), (21), (22) and (23) are embodied in an algorithm called E U L E R (Abate and Whitt, 1992a), which is shown in Figure 1. The Laplace transform variable s is written as s = x + iy and f (s) is written as f (x, y).

412

D. Heyman l. INPUT t 2. Set Parameters: A = 19.1, n = 15, m = 11, U = exp(A/2)/t, c(1) = 1 3. Set Variables: x = A/(2t), h = ~ / t 4. Sum First n Terms of (18): sum = ~ f ( x , 0): For k = 1,2 . . . . . n (a) y -----kh (b) sum ~-- sum 4-(--1)kf(x, y) (c) end for 5. Sum Next m terms of (3) by Euler Summation: Compute Constants: (a) For k = 2 , 3 . . . . . m + l (b) c ( k + 1) = c ( k ) ( m - k + 1)/k (c) end for 6. Euler Summation: s(1) = sum (a) (b) (c) (d) (e)

For k = 1,2 . . . . . m + l j=n+k y = jh s ( k + 1) = s ( k ) + (--1)if(x, y) end for

7. a n s = 0 (a) For k = 1,2 . . . . . m + l (b) ans <-- ans + c ( k ) s ( k + 1) (c) end for 8. ans <-- ans x (U/2 m) 9. OUTPUT: f(t) ---- arts

Fig. 1. Algorithm EULER. The Fourier analysis, trapezoidal rule integration and Poisson summation techniques can also be applied to generating functions (Abate and Whitt, 1992b). Let OO

~' (Z) = Z zkqk'

(24)

k=0 where z is complex. When q is a bound [qkl ~< 1 which is important this bound, ~ is analytic inside the z inside the region of convergence. q0=~(0).

probability mass function, we have the convenient for error analysis but is otherwise inessential. With unit circle. Assume that c) can be computed for any Then (25)

Numerical methods in queueing theory

413

The inversion formula is

qk= ~

G ( r ) + G ( - r ) + 2E(-1)J~t[G(re~Ji/k

k=l,2

.....

j:l

(26) where 0 < r < 1. The error in (26) is at most r2k/(1 -- r2k), which is close to the numerator when the numerator is small. The ideas behind the inversion algorithms described above have been extended to produce algorithms to invert multi-dimensional Laplace transforms, including lattice functions (Choudhury, Lucantoni and Whitt, 1994). In particular, transforms such as

P(z,s) = E z

j

°° f0~

e-Stp(t)dtj

(27)

j=0

can be inverted when each pj is a probability mass function.

3. The ubiquity of Markov chains

Markov chains are ubiquitous in the numerical solution of queueing problems because there is a lot known about numerical solution of Markov chains and because many queueing problems can be described by a Markov chain. The reason for the latter is that even when interarrival or service times are not exponentially distributed, they can be approximated by sums and mixtures of exponential distributions. This idea started with Erlang's method of stages (Heyman and Sobel, 1982, p. 298) and has been extended in recent years. Neuts (1975, 1981 ) introduced the phase-type distributions, which are firstpassage-times of an absorbing continuous-time Markov chain. In addition, some nonMarkovian queueing processes have an embedded Markov chain, such as the M / G / 1 queue. Another way Markov chains arise is that non-Poisson arrival processes are often modeled with a Markov-modulated Poisson Process (MMPP) in which the state of an ergodic continuous-time Markov chain determines the rate of a Poisson stream of events. Each event can be the arrival of a batch of customers. The most general process of this type is the batch Markovian arrival process (BMAP), which is described in Section 6. It can be a difficult task to turn a Markovian model as described above into a computer program to construct the transition matrix. Berson, de Souza de Silva and Muntz (1991) describe a way to generate the transition matrix from a description of the system dynamics.

D. Heyman

414 4. Finite Markov chains

4.1. Transientprobabilities For a continuous-time Markov chain with generator Q, let Hij (t) = P[X (t) =- j [X (0) = i], and 17(0 = (Hij (t)). The balance equations are /:/=/TQ,

t>0,

(28)

where the dot denotes differentiation and H(0) = I. When the structure of Q permits, a transform of the transient probabilities can be obtained and numerical transform inversion applied, as illustrated by the transient probabilities of the M / G / 1 queue described in Section 1. The solution of (28) is / 7 ( t ) = e Qt,

t~>0,

(29)

and for a finite number of states, one can "just compute the matrix exponential" and the solution is obtained. Computing a matrix exponential accurately is not easy. There are at least 19 dubious ways (Moler and Van Loan, 1978); see Golub and Van Loan (1989) for a good way. One could apply numerical differential equation solvers to (28). These methods may display numerical inaccuracies caused by rapidly decaying components of the transient solution (i.e., stiffness). See (Stewart, 1994, Section 8.4) for an introduction to these methods. The most commonly used method is Jensen's algorithm (Grassmann, 1991; Stewart, 1994, Section 8.2) which is based on the uniformization procedure for continuous-time Markov chains. The uniformization procedure (Heyman and Sobel, 1982, Section 87) replaces the generator Q with the transition matrix P = I ÷ Q/qmax, where qmax = maxi {-qii }. Uniformization is used repeatedly in the sequel; it provides both computational and analytic results. Let N(t) be the the number of transitions that occur by time t; it has a Poisson distribution with mean qmaxt and

flij (t) -- ~ Pkp[N(t) -- k].

(30)

k=O Numerical implementations of (30) replace the infinite sum with a finite sum with upper limit K say. The truncation error is

Pkp[N(t)=i] <~ Z k=K+l

P[N(t)=i]=P[N(t) > K]

(31)

k=K+l

which is kept below some number ~ > 0 when K is large enough so that

~ k=0

(qmaxt)k -->/(1--8)e k!

qm~t.

(32)

Numerical methods in queueing theory

415

1. Calculate K from (32) (a) S e t K = 0 , ~ = 1 , o r = l , q = ( 1 - s ) e (b) While o- < r/do i. K < - - K + 1 , ~<--~ex qmaxt K " 2. Calculate

qmaxt.

p(t)

(a) p <-- p(0), y = p(0) (b) F o r k = l , 2 . . . . . K d o qmaxt i. y ~-- yP x k ii. p + - p + y (c) p(t) = e-qmaxtp Fig. 2. Jensen's algorithm.

Since there are no subtractions in (30) it should be resistant to computer round-off errors (Grassmann, 1983), so the computed probabilities should differ from the exact probabilities by no more than E. From (32) we see that for e fixed, the computational effort grows exponentially with the product qmaxt. To mitigate this effect, one can divide the interval (0,t] into subintervals and solve a smaller problem on each subinterval. For 0 = to < t~ < t2 < . . . < t, = t, compute ri = ti - ti 1 and gi 17(ti) = / - / ( t i - 1 ) Z pke-qmaxri (qmax'Ci)k k! '

i = 1, 2 . . . . . n.

(33)

k=O

Figure 2 shows Jensen's algorithm (Stewart, 1994, p. 413) for computing the vector p ( t ) = p ( O ) I 7 ( t ) where p(0) is given vector of initial probabilities. For each time point (ti) that Jensen's algorithm evaluates, the number of multiplications (which is the major computational expense) is K ( s + Sp) where s is the number of states in the Markov chain and Sp is the number of positive entries in P. The vector

p(0), the matrix P and two work vectors of length s need to be stored. 4.2. Steady-state p r o b a b i l i t i e s

Since transient probabilities are often difficult to compute, and since steady-state probabilities, when they exist, describe long-run averages, steady-state probabilities are often computed. The statement of the problem is this. The states of the Markov chain are numbered from 1 to N, and the transition matrix is P = ( P i j ) . We assume P is ergodic and seek a solution of ~=~P,

(34)

with N

Jr~>0

and

ZTri=l' i=1

where Jr is a row vector.

(35)

416

D. Heyman

There are two general classes of algorithms for this problem. Direct algorithms compute an exact solution in exact arithmetic (no roundoff errors) and the number of arithmetic operations can be computed from the description of the problem. An iterative algorithm produces a sequence of approximations that is expected to converge to the exact solution. (Some algorithms appear to converge when applied to problems that don't satisfy the conditions that have been needed to prove convergence.) Gaussian elimination can be applied to (34) and (35); it will require N3/3 (plus terms in N 2 and N) floating point operations. (A floating point operation, or flop, is the arithmetic required to do ai 4-- ai + bj × Ck, which is an addition, a multiplication and some subscripting.) A variant of Gaussian elimination designed for this system of equations is described in Section 4.2. Since rr(0)limk-+~ pk = Jr for any initial probability vector 7r(0), forming successive powers of P is a convergent algorithm which is called the power method. Since computing P x P requires N 3 flops in general, it cannot require fewer flops than Gaussian elimination unless the structure of P can be exploited. The rate of convergence of the power method is controlled by the eigenvalue of P with second largest modulus, and can be very slow. This is particularly noticed with nearly completely decomposable Markov chains (Courtois, 1977). (A nearly completely decomposable Markov chain has a block partitioning such that the blocks along the diagonal have row sums that are close to one. These blocks describe local changes of state which occur much more frequently than transitions to states in other blocks. These Markov chains occur, for example, with MMPP arrivals when many changes in the queue length occur between changes in the arrival rate. Algorithms have been devised that exploit this property (Stewart, 1994, Chapter 6).) Jacobi and Gauss-Seidel iterations are based on computing the iterations Jr (n + 1) +7c(n)P in an effective way. Successive over-relaxation (SOR) methods attempt to speed convergence by averaging successive iterates: rcnew +- x ( n ) P

and (36)

rr(n + 1) +-- cozcnew + (1 - co)Jr(n),

0<0)<2.

The choice of co determines the rate of convergence, and selecting good values is currently more of an art than a science. Seelen (1992) had success in using successive over-relaxation on queueing problems. He developed ways to adapt co as the computations proceed. Since premultiplying a matrix by a vector takes N 2 flops, these iterations should be used when P is sparse or so large that it cannot be stored all at once. Greenberg and Vanderbei (1991) find that these iterations require approximately N iterations to converge, so Gaussian elimination will require less computation unless P is sparse and Gaussian elimination cannot take advantage of this. A direct algorithm The direct method of choice for solving the steady-state balance equations is the GTH algorithm, which was introduced in Grassmann, Taksar and Heyman (1985). It is a

Numerical methods"in queueing theory

417

variant of Gaussian elimination that accurately computes the stationary vector of an irreducible, finite stochastic matrix. The accuracy is achieved by avoiding subtractions, which are a main source of computer round-off errors. Empirical evidence of its accuracy is given in Heyman (1987), and analytic evidence in O'Cinneide (1993). If there are transient states in P, the algorithm will either compute their steady-state probabilities as zero or attempt a division by zero. The latter can be used to discard transient states (Heyman, 1987). The GTH algorithm also finds the steady-state probabilities of a continuous-time Markov chain with an irreducible generator. This follows directly by uniformizing the continuous-time process and treating it as a discrete-time Markov chain, and is described in more detail after the description of the algorithm. The probabilistic reasoning behind the GTH algorithm is this. Suppose state N is invisible, so the sequence of transitions i --+ N --+ j appears to be a transition from i directly to j. With this supposition, the Markov chain with N states would appear to be a Markov chain with N - 1 states, and the transition probabilities would appear to be

Pij + PiN(I + PNN + p 2 N + ' ' ' ) p N j -

PiNPNj

i,j
(37)

the terms in parentheses account for the instantaneous transitions from state N back to itself. This operation is called state reduction. Algebraically, (37) is solving the Nth equation in (34) for rrN and substituting the results in the other equations. The state reduction step can be repeated until there is only one state left in the Markov chain; since there is exactly one redundant equation in (34) when P is irreducible, the result of all the state reduction steps is the equation Jr1 = re1.

(38)

Any positive solution of (34) can be normalized to sum to one to satisfy (35), so we pick 7q = 1. The back-substitution steps reverse the state reduction steps, computing rri from rej, j < i. The normalization step makes the probabilities sum to one. We display the GTH algorithm in Figure 3. The notation a +-- b means b replaces a. For continuous-time Markov chains, the uniformization procedure allows the algorithm to be applied to the transition matrix P = I + Q/qmax. The normalization step cancels the effect of carrying qmax through the state reduction step, and since the diagonal terms are never used, one can apply the algorithm directly to the generator Q. The state reduction step overwrites P so the algorithm has no storage requirements; the matrix P has to be stored if it is to be used in future calculations. The backsubstitution and normalization steps can be combined so overflow of TOT or underflow of zrj are avoided. In Step l(a), Sn should be tested to see if it is zero to machine precision. If it is, then (to machine precision) P is such that states n + 1, n + 2 . . . . . N cannot be reached from state n. The algorithm requires N3/3 additions and an equal number of multiplications, plus terms of order N 2 and order N. The following important special structure can be exploited to reduce the number of computations. P is called banded when either

D. Heyman

418

1. (State Reduction) For n = N, N - 1 . . . . . 2, do the following: (a) Let Sn = Y~,jn11 Pnj. (b) Let Pin ~- Pin/Sn, i < n. (c) Let Pij <--- Pij + PrnPnj, i, j < n. 2. (Initialization) Initialize TOT = 1 and ~1

=

1.

3. (Back-substitution) For j = 2, 3 . . . . . N do the following: j-1 (a) Let rrj = Plj 4- ~ k = 2 rrkPkj • (b) Let TOT ~-- TOT 4- ~j. 4. (Normalization) Let rrj <-- rq/TOT, j = 1,2 . . . . . N.

Fig. 3. GTH algorithm.

Pi, i+k = 0 whenever k > g or Pi,i-k = 0 whenever k > h for some positive numbers g and h (with appropriate modifications at the boundaries - when i ÷ k > N or i - k < 0, the condition is valid). This means that the chain cannot jump too far from the current state. Two examples are the M / M / 1 / k queue and the Engset model with several arrival classes (Cooper, 1981, Section 3.7). The overflow model in Section 1 with generator (8) has a particular banded structure called block tri-diagonal. There is no need to perform the calculations in Steps 1 and 2 that are outside of the band.

5. Infinite Markov chains Markov chains with an infinite number of states often are more analytically tractable than the same model with a finite number of states because a boundary condition that is in the finite state version of the model is absent from the infinite state version. The M / G / 1 queue is an example. Structured infinite-state Markov chains may also exhibit this property, which can be exploited in numerical algorithms. This work started in the mid 1960's (Wallace and Rosenberg, 1966; Evans, 1967) and gained impetus a decade later when computing became easier.

5.1. The main paradigms We assume that the Markov chains under study are irreducible and persistent, so that they possess a unique stationary distribution. For discrete-time chains, we do not require aperiodicity; if this property is present, the stationary distribution is also a limiting distribution. Transition matrices and generators are expressed in block partition form, P (or Q) = (Aij) say, where Aij is a matrix (not necessarily square). This corresponds to a state space of ordered pairs, (i, k) say, where i runs from zero to infinity and is called the level and k has a finite support and is called the phase. The interpretation of Aij is model specific, but the following example shows the general idea. Suppose the arrivals form a 2-state MMPP, the service times are exponential, and there is a single server. Let )~k be the arrival rate in phase k and qk be the rate the MMPP leaves phase k.

Numerical methods in queueing theory

419

Let/~ be the service rate. Then i and j represent the number of customers present and phase respectively, and the generator has Ai i-1 :

Ai i :

'

'

Ai,i+l =

'

Yl q2

Y2

'

( 010)

)~2 '

where Yl and Y2 make the row sums equal zero. The structure that leads to explicit solutions is repeating rows, except possibly for the first row and first column to allow for different boundary conditions. Specifically, the structure is a block version of a Toeplitz matrix,

C1 C2 C3

Ao A-1 A-2

A1 Ao A-1

A2 A1 Ao

(40)

A formal representation of the stationary distribution of (40) as a transition matrix can be obtained (Grassmann and Heyman, 1990), and an algorithm for its computation developed (Grassmann and Heyman, 1993). When the matrices to the left of A_~ are all zero (block upper-Hessenberg or skip-free-to-the-left in blocks) or the matrices to the right of A1 are all zero (block lower-Hessenberg or skip-free-to-the-right in blocks) an explicit solution can be obtained. The former is called the M / G / 1 paradigm (Neuts, 1989) because it's the block form of the Markov chain used to analyze the M / G / 1 queue that is shown in (1). The latter is called the G I / M / 1 paradigm (Neuts, 1981) for an analogous reason. When the chain is skip-free in both directions it is called a quasi-birth-and-death (QBD) process because it's the block form of the generator of a birth-and-death process. We now examine algorithms for these three special cases. 5.2. The G I / M / 1 paradigm We change the subscripting used in (40) to make the subsequent equations neater. The specific matrix we study here is

Co CI C2

m H1 H2

0 H0 H1

0 0 H0

...| ... ] .

? :::)

(41)

Let (41) be a transition matrix, let re = (7r0, Jrl,rr2 . . . . ) be the stationary vector partitioned conformally (7ri is associated with the row that starts with Ci). In the scalar G I / M / 1 queue, the (scalar) stationary probabilities satisfy 7ri+l = rTri where r is the

D. Heyman

420

root of x = Z ~ x i H i inside the unit interval. Thus, rci = (1 - r)r i, which is the geometric distribution. Suppose an analogous condition applies to (41); specifically, suppose there is a matrix R such that 7ri+l =7ciR .¢:>yri =rcoR i,

i >~0.

(42)

This is called a matrix-geometric distribution. The balance equation for zr0 is O<3

7rO= YrOE Rk Ck" 0

(43)

The balance equation for rc~ is oo

oo

(44) 0

0

From (44), (42) can be valid only if

R = ERkHk . 0

(45)

Substituting (45) in the balance equation for re2 yields Oo

G~

:

=

0

nor

R

:

2

(46)

0

Proceeding in this way yields (42). In order for oO

E0 rrkl =

1

(47)

(1 is a column vector of ones), we must have limk-+oc R lc = 0, which implies that ( I - R) -1 exists and all of the eigenvalues of R are inside the unit disk. These results are formally proven in (Neuts, l 981, Chapter 1) and displayed in Theorem 1.2.1 there, along with the additional results that the matrix ~ c RkCk appearing in (43) is stochastic and that R is the minimal non-negative solution of

X = ExkHk. 0

(48)

The matrix R can be found from (48) by successive substitution starting with the initial solution X = 0. This produces an increasing sequence of matrices converging to

Numerical methods in queueing theory

421

R (Neuts, 1981, Lemma 1.2.3). Various other algorithms for computing R have been proposed; Latouche (1993) compares them. Once R is computed, (43) is solved, and then re0 is normalized by (47). The state probabilities then follow from (42). From the matrix-geometric solution for the state probabilities, one can compute that the mean number of customers present in the steady-state is

L = ~--~kzr~l = 7r0 0

k R k 1 = rc0(I - R ) - l l .

(49)

The matrix R has a probabilistic meaning that can be deduced from (45). The skipfree-to-the-left structure implies that when there is a sequence of transitions such that the levels go monotonically from n + k to n (k > 1), when the levels decrease, they decrease by one. Since the rows repeat, the probabilistic laws governing these transitions do not depend on n. Let rij be the mean number of times the chain enters state (n -t- 1, j ) before returning to level n when starting from state (n, i), and R (rij). The repeating rows implies that the (i, j ) t h element of R g is the mean number of visits to state (n + k, j) before returning to level n when starting in state (n, i), and conditioning on the first transition will show that R satisfies (45). In computations, all Hk with k greater than some finite n are zero. The usual stopping criterion for iterating (48) is that the difference between successive iterates is smaller than a given tolerance. This criterion doesn't guarantee that the computed R is within a known, or even small, distance of the true R. These results for a transition matrix formally carry over to the generator of a continuous-time Markov chain by uniformization (Neuts, 1981, Section 1.7). The numerical values are not the same, but (43) with the left side replaced by zero, (47) and (48) are valid. =

5.3. The M / G / I paradigm

The block matrix analog of the embedded Markov chain of the M / G / 1 queue is Bo Ao 0

B~

B2

B3

A1 Ao

A2 A1

A3 A2

...'~ ...1 ...].

o o ? ? :::)

(5o)

Permitting the first row to differ from the second row allows one to model various boundary conditions, such as letting the first customer served in a busy period to have a different service-time distribution from all other customers, and introduces only minor notational changes in the analysis. As above, let (50) be an ergodic transition matrix, and 7r = (7r0, rcl, re2 . . . . ) be the stationary vector partitioned conformally. The skip-free-to-the-left property leads to an analog of (45). Let gij be the probability that when starting in state (n + 1, i), n > 0, the state entered when the Markov chain next

D. Heyman

422

reaches level n is (n, j), and G to

G = ~

=

(gij).

Then G is the minimal non-negative solution

(51)

AkG k o

and G is stochastic (Neuts, 1989, Chapter 2). The solution of (51) can be obtained by successive iterations. Experience has shown that convergence is most rapid when the initial trial solution is a stochastic matrix (Lucantoni, 1993). The matrix G doesn't suffice to compute the steady-state distribution for the M / G / 1 paradigm as the matrix R did for the G I / M / 1 paradigm. The reason is the following. Temporarily use scalar states and let P be a irreducible transition matrix with positiverecurrent states. Then

(52)

I - P = (A - I ) ( B - S),

where A is strictly upper-triangular, B is strictly lower-triangular, and S is diagonal with the upper-left element equal to zero and all other diagonal elements positive; (Grassmann, 1993) for a finite number of states and (Heyman, 1995) for infinitely many states. This is the analog of the L U factorization that arises when solving systems of linear equations. To obtain Jr, one first solves for a row vector ot that satisfies o~(B - S) = 0. Since the top row of B - S is zero, ~ = (d, 0, 0 . . . . ) where d is any number is a solution. Then one solves rc(A - I) = ee

(53)

for 7r. Choosing d = - 1 yields 7r0 = 1 and since A - I is upper-triangular (53) is solved recursively, yielding the unnormalized solution j-1

~j=Z~iaij,

(54)

7to = 1, A = (aij).

i=0

In the the M / G / 1 paradigm, the matrix G determines the matrix B in the factorization. In the G I / M / 1 paradigm, the matrix R determines the matrix A in the factorization and (54) yields the matrix-geometric solution. The requisite formulas for obtaining rc for the M / G / 1 paradigm were derived by Ramaswami (1988). For Jr = (re0, 7rl . . . . ) partitioned conformally with the transition matrix,

7rj+

rroBj(G)+~-~zriAj+l_i(G) i=1

1{

I-AI(G)]

-1,

j>~l,

(55)

Numericalmethodsin queueingtheory

423

where Oo

B j ( G ) = E BiGi-J

and

A j ( G ) = E AiGi-J,

i=j

j ~> 0.

(56)

i=j

One can obtain the generating function of Jr just as in the scalar M / G / 1 queue. Define the generating functions oO

~

oQ

k=O

k=O

k=O

(57) for Lzl ~< 1. Write the balance equations in block form, multiply the kth equation by z k and sum, to obtain ~(Z) = Jr0B'(z) -]- JrlA'(z) + Z.Yr2A'(z) -'b Z27t'3A'(z) q - " "

(58)

whence

(z)[zi - 2(z)] = Jr0[za(z)-

(59)

Obtaining Jr0 is beyond the scope of this presentation; the formula and derivation are given by (Neuts, 1989, Theorem 3.2.1). This generalizes the Pollaczek-Khinchine formula for the scalar case. The difficulty in using (59) in numerical work is the computation of the matrices A'(z) and B(z). The next section describes a special case in which this computation is not difficult. The special case is broad enough to include a wide variety of useful models.

6. The BMAP/G/1 queue The batch Markovian arrival process (BMAP) was introduced by Neuts (1979). The BMAP/G/1 queue is a special case of the M / G / 1 paradigm whose special structure allows the same performance measures to be computed as can be computed for the traditional M / G / 1 queue. The BMAP includes a wide variety of non-Poisson processes that have been found useful in many applied studies, so this queueing model is frequently used. The first detailed solution of the BMAP/G/1 queue was given by Ramaswami (1980). Further results that make the computations tractable were obtained by Lucantoni (1991). This presentation follows Lucantoni (1993). A BMAP is built on an underlying continuous-time Markov chain on a finite state space; these are the phases of the BMAP. Let D be the generator of this Markov process. When the phase process makes a transition to state i, a random number of arrivals occur, possibly zero. Let )~i be the transition rate in phase-state i and Pi (k, j) be the probability that when state i is entered, k arrivals occur and the next phase state will be j. Define the matrices Dk by (Do)ii = --)~i, 1 <~i <~m < ~x), (Dk)ij = L i P i ( k , j ) , i 7£ j and

D. Heyman

424

k/> 0. Then D = } - ~ Dk. Let 3 be the stationary distribution of the generator D and ~. be the arrival rate of the BMAP. Then oo

(60)

= ~ Z kDkl. k=0 The matrix generating function (x)

D(Z) = Z zk Dk' k=0

[zl ~ 1,

(61)

is a key quantity. Let N(t) be the number of arrivals in (0, t] and J(t) be the phase at time t. The process {N(t), J(t); t ) 0} is a continuous-time Markov chain with generator Do

D1 Do 0

0

0

D2 D1 Do

D3 D2 D1

?..o.. ?.

... ... ...

(62)

:::

where the cumulative number of arrivals is the level. Let Pij (n, t) = P [ N (t) = n, J (t) = j I N ( 0 ) = 0, J(0) = i], P(n, t) = (Pij(n, t)), and OO

fi'(Z, t) -= ~

zn P(n, t) = e D(z)t,

[zl ~ 1, t ~> O.

(63)

n=0

It is this exponential form for P which makes the analysis of the B M A P / G / 1 queue tractable; it's analogous to the exponential form of the probability generating function of a Poisson process and permits a useful explicit formula for A(z). Let H be the distribution of the service times, then

A(z) =

fo £

z n P (n, t) dH (t) =

n=0

=

fo

"fi(z, t) d H (t)

e D(z) d H ( t ) = / 4 [ - D ( z ) ] .

(64)

f0 °°

When H has a rational Laplace-Stieltjes transform, say

hi(s) - y~'~ aksk Y2~ bk s k '

m <~n,

(65)

Numerical methods in queueing theory

425

then in

(66)

which requires evaluating two matrix polynomials and one matrix inverse. Among the special cases of a BMAP are renewal processes with Erlang, hyperexponential and phase-type distributions, and MMPPs. An MMPP with underlying phase transition process with generator Q and arrival rates )~i is the BMAP with A = diag(kl . . . . . )~m), Do = Q - A, D1 = A, and D~ = 0 for k > 1. For BMAP arrivals, one can obtain

Bn -- - D -1 ~

Dk+lAn-k

and

k=0 B(Z) = - D o 1 [D(z) - Do]A'(z).

(67)

Substituting (67) into (59) yields

(Z) [ZI -- A'(Z)] ----- 7 r 0 D o 1D(z)A(z).

(68)

To evaluate 7r0, first obtain the matrix G, and the steady-state distribution of G, F say. Let/2 -1 be the mean service time; the traffic intensity of the queue is p = 2/fi, and p < 1 is the ergodicity condition. Then (Lucantoni, 1991, Eq. (54)) fro = (1 - P)F(-Do)/~..

(69)

Taken together, (64), (68) and (69) provide an explicit solution for the probability generating function of the steady-state probabilities. The Laplace transforms of the waiting-time distribution and transient state probabilities can be obtained by generalizing the analysis of the scalar case. Since the arrivals are not a pure Poisson process and may come in batches, care must be taken to distinguish probabilities at customer arrival epochs, at customer departure epochs and at an arbitrary time. Numerical values are then obtained by numerical transform inversion, so both of the general methods described in this article are used together. Special cases such as arrivals occuring one at a time, as a renewal process with a phase-type inter-arrival time distribution, service times having a phase-type distribution, and arrivals forming a 2-state MMPP produce more explicit formulas (Lucantoni, 1991).

7. The quasi birth-and-death process The quasi-birth-and-death process (QBD) is skip-free in both directions, so it enjoys the special properties of the G I / M / 1 and M / G / 1 paradigms, but it's more than just a special

426

D. Heyman

case. By making the blocks large enough, infinite if necessary, the QBD can describe the general case of repeating rows. Latouche and Ramaswami (1999) develop the theory in this fashion. They also show how to exploit the special structure of the blocks to obtain efficient algorithms for the G I / M / 1 and M / G / 1 paradigms. The discussion here is limited to describing how the special structure of the QBD leads to a very fast algorithm for the state probabilities. The QBD is used frequently in applied studies because in prospective studies one often only has estimates of the mean arrival rate and service time, so exponential assumptions are made. The specific matrix we study here is Co Cl 0

H0 H1 //2

0 H0 H1

0 0 HO

...'~ ...1 ...1.

(70)

? 0 ? ? :::)

We interpret this as a transition matrix; continuous-time Markov chains are analyzed by uniformization. From (45) and (51) we have R = Ho + RHI + R2H2

(71)

G = Ho + H1G + H2 + H2G 2.

(72)

and

Let Uij be the probability that starting in level n > 0 and phase i, the next visit to level n starts in phase j and there are no visits to level n - 1 in the interim. The repeating rows property implies that this probability does not depend on n. Then U is substochastic and is the minimal non-negative solution of U = H1 -Jr-H O ( I Jr- U -~- U 2 -Jr-" ")H2 = H1 Jr- HO(I -

U)-IH2.

(73)

These fundamental quantities are related by G ----(I - U ) - I H 2 ,

(74)

R = Ho(I - U) -~

(75)

U = HI + HoG = H1 + RH2.

(76)

and

From the discussion of the G I / M / 1 paradigm we know we can compute the steadystate probabilities from R. Rather than compute R from (71), one uses a rapidly converging recursion for G, then computes U fi'om G using (76), and then computes

N u m e r i c a l m e t h o d s in q u e u e i n g theory

427

1. (Initialize) (a) (b) (c) (d)

Set i = 0. Compute B 0 = ( I - H 1 ) - I H o . Compute B 2 = (I - H 1)-1 H2. Set G = B2 and P = Bo.

2. (Compute G) While mini ~ j Gij < 1 - e do the following:

(a) i <--i+ 1. (b) Compute H~ = BOB2, H~ = B 2 and H i = B22. (c) Compute BO = (I -- H1)*-1 H 0* and B2 = ( I - H1,)-1 H2., (d) G <-- G + PB2 and P +- PBo. 3. (Termination) (a) Compute U = H1 + HoG. (b) Compute R = HO + (I - U) -1 . Fig. 4. Logarithmic reduction algorithm.

R from U using (75). The new recursion for G was developed by Latouche and Ramaswami (1993). They consider censored versions of the Markov chain that visits only even numbered states. This yields a sequence of substochastic matrices Gn that converge monotonically to the stochastic matrix G;

Gn = £ ( I - I B ~ i ) )

"(k)"

G = lim Gn f t --+ o~

k=O \ i=0

(77)

where B}°)

= (I - H1)-IHi,

i = O, 2,

(78)

and B}k+l)

= (I - - ~n(k) n(k) - - ~n(k) n(k)~-I 0 ~2 2 ~0 J (B{k)) 2

(79)

for i = 0, 2 and k = 0, 1,2 . . . . . (An empty product is defined to be I.) Since Gn is substochastic, a convenient convergence criterion is to terminate the iterations when all of the row sums are larger than 1 - e for some specified small number e > 0. The algorithm is displayed in Figure 4. When there are m rows in each block, an iteration of this algorithm requires (25/3)m 3 + O(m) flops. Computational experience has shown that this algorithm converges very rapidly and that all the computations can be done accurately.

428

D. Heyman

References Abate, J. and W. Whitt (1992a). The Fourier-series method for inverting transforms of probability distributions. Queueing Systems 10, 5-88. Abate, J. and W. Whirr (1992b). Numerical inversion of probability generating functions. Oper. Res. Lett. 12, 245-251. Berson, S., E. de Souza de Silva and R. R. Muntz (1991). A methodology for the specification and generation of Markov models. In Stewart (1991), pp. 11-36. Choudhury, G. L., D. M. Lucantoni and W. Whitt (1994). Multidimensional transform inversion with application to the transient M/G/1 queue. Ann. AppI. Probab. 4, 719-740. Cohen, J. W. (1969). The Single Server Queue. North-Holland, Amsterdam. Cooper, R. B. (198i). Introduction to Queueing Theory, 2nd ed. North-Holland, New York. Courtois, R-J. (1977). Decomposability. Academic Press. Evans, R. V. (1967). Geometric distribution in some two dimensional queueing systems. Oper. Res. 5. Franklin, E (1964). A Treatise on Advanced Calculus. Dover, New York. Golub, G. H. and C. F. Van Loan (1989). Matrix Computations, 2rid ed. Johns Hopkins University Press, Baltimore. Grassmann, W. K. (1983). Rounding errors. Technical report, Operations Research Department, Stanford University. Grassmann, W. K. (1991). Finding trmlsient solutions in Markovian event systems. In Stewart (1991), pp. 35737l. Grassmann, W. K. (1993). Means and variances in Markov reward systems. In Meyer and Plemmons (1993). Grassmann, W. K. and D. R Heyman (1990). Equilibrium distribution of block-structured Markov chains with repeating rows. J. Appl. Probab. 27, 557-576. Grassmann, W. K. and D. R Heyman (1993). Computation of steady-state probabilities for infinite-state Markov chains with repeating rows. ORSA Z Comput. 5, 292-303. Grassmann, W. K., M. I. Taksar and D, R Heyman (1985). Regenerative analysis and steady-state distributions for Markov chains. Oper. Res. 33, 1107-1116. Greenberg, A. and R. Vanderbei (1991). Quicker convergence for iterative numerical solutions to stochastic problems: probabilistic interpretation, ordering heuristics, and parallel processing. Probab. Engrg. Inform. Sci. 4, 493-521. Heyman, D. R (1987). Further comparisons of some direct methods for computing stations distributions of Markov chains. SIAM J. Alg. and Disc. Meth. 8, 52-60. Heyman, D. R (1995). A decomposition theorem for infinite stochastic matrices. J. Appl. Probab. 893-901. Heyman, D. R and M. J. Sobel (1982). Stochastic Models in Operations Research, Vol. 1. McGraw-Hill, New York. Latouche, G. (1993). Algorithms for infinite Markov chains with repeating columns. In Meyer and Plemmons (1993). Latouche, G. and V. Ramaswmni (1993). A logarithmic reduction algorithm for the quasi-birth-and-death process. J. Appl. Probab. 30, 650-674. Latouche, G. and V. Ramaswami (1999). Introduction to Matrix Analytic Methods in Stochastic Modeling. SIAM. Lucantoni, D. (1991). New results on the single server queue with a batch Markovian arrival process. Stochastic Models 7, 1-46. Lucantoni, D. (1993). The BMAP/G/1 queue: a tutorial. In Models and Techniques for Performance Evaluation of Computer and Communication Systems, pp. 330-358 (Eds. L. Donatiello and R. Nelson). Springer-Verlag. Meyer, C. D. and R. J. Plemmons (Eds.) (1993). Linear Algebra, Markov Chains and Queueing Models. Springer-Verlag. Moler, C. and C. Van Loan (1978). Nineteen dubious ways to compute the exponential of a matrix. SIAM Rev. 20, 801-836. Neuts, M. E (1975). Probability distributions of phase type. In Liber Amicorium Prof. Emeritus H. Florin, pp. 173-206. University of Louvin.

Numerical methods in queueing theory

429

Neuts, M. E (1979). A versatile Markovian point process. J. Appl. Probab. 16, 764-779. Neuts, M. F. (1981). Matrix-Geometric Solutions in Stochastic Models. The Johns Hopkins University Press. Neuts, M. E (1989). Structured Stochastic Matrices of M / G /1 Type and their Applications. Marcel Dekker. O'Cinneide, C, A. (1993). Error analysis of a variant of Ganssian elimination for steady-state distributions of Markov chains. Numer. Math. 65, 109-120. Ramaswami, V. (1980). The N/G/1 queue and its detailed analysis. Adv. Appl. Probab. 12, 222-261. Ramaswami, V. (1988). A stable recursion for the steady-state vector in Markov chains of M/G/1 type. Stochastic Models 4, 193-188. Seelen, L. E (1992). An algorithm for Ph/Ph/c queues. Eu~ J. Oper. Res. 23, 118-127. Stewart, W. J. (Ed.) (1991). Numerical Solution of Markov Chains. Marcel Dekker. Stewart, W. J. (1994). Introduction to the Numerical Solution of Markov Chains. Princeton University Press, Princeton. Takfics, L. (1964). Priority queues. Oper. Res. 12, 63-74. Wallace, V. L. and R. S. Rosenberg (1966). Markovian models and numerical analysis of computer system behavior. In Proc. AFIPS Spring Joint Computer Conference, Vol. 28. AFIPS Press, New Jersey.